Neutrino Portals, Terrestrial Upscattering, and Atmospheric Neutrinos

We consider the upscattering of atmospheric neutrinos in the interior of the Earth producing heavy neutral leptons (HNLs) which subsequently decay inside large volume detectors (e.g. Super-Kamiokande or DUNE). We compute the flux of upscattered HNLs arriving at a detector, and the resultant event rate of visible decay products. Using Super-Kamiokande's atmospheric neutrino dataset we find new leading constraints for dipole couplings to any flavor with HNL masses between roughly 10 MeV and 100 MeV. For mass mixing with tau neutrinos, we probe new parameter space near HNL masses of $\sim 20$ MeV with prospects for substantial future improvements. We also discuss prospects at future experiments such as DUNE, JUNO, and Hyper-Kamiokande.


I. INTRODUCTION
Heavy neutral leptons (HNLs) are well motivated extensions of the Standard Model (SM).They appear ubiquitously in dark sector models, and are especially important because their coupling to neutrinos represents one of three unique renormalizable "portals" between a generic dark sector and the Standard Model [1][2][3][4].They are natural partners to left-handed neutrinos, reflecting a matching of chiral degrees of freedom observed among all other SM fermions, and their non-observation is easily explained due to their being a SM gauge singlet.They are further motivated by anomaly-cancellation arguments that become essential in models with new gauge groups e.g. a gauged B − L [5,6].HNLs appear as necessary ingredients in certain grand unified theories (GUTs) e.g.SO (10) [7], and they are intimately connected to neutrino masses. 1  Despite their strong theoretical motivation, there is no model-independent prediction for the HNL mass scale.HNLs may be O(eV) in mass (i.e.sterile neutrinos) and connected to neutrino masses via a naive Dirac mass mechanism (like all other SM fermions), in which case they are most easily searched for using short-baseline oscillation experiments or cosmological observables.Alternatively, neutrino masses may be generated by a type-I [11][12][13] seesaw mechanism each of which leads to different expected mass scales for different Yukawa couplings.
The ubiquity of HNLs in generic models of a dark sector, and their unconstrained mass range therefore motivates a broad search strategy that targets many decades of HNL mass parameter space, ranging from the eV to the GeV (or even TeV) scale [14][15][16][17][18][19][20][21][22][23][24][25][26][27].In this work we will demonstrate that large volume detectors can efficiently search for HNLs via their decays by leveraging atmospheric neutrino upscattering inside the Earth.This search strategy is ideally suited to HNL decay lengths, λ, satisfying 10 m ≲ λ ≲ 6000 km with the upper limit set by the Earth's radius.This complements fixed target "beam dump" experiments and missing energy searches which typically lose sensitivity as the HNL decay length becomes much longer than the experimental apparatus, which tend to range from 10s to 100s of meters (much shorter than the 1000s of km that characterize the Earth's radius).We derive new and leading constraints for HNL masses between ∼ 10 MeV and 200 MeV, a mass range which has interesting implications for models that address the Hubble tension [28].
Much of the literature on HNLs focuses on the aforementioned renormalizable coupling between HNLs and active neutrinos available in the SM.This is the so-called mass-mixing portal, which results in a small admixture of HNL contamination among the active neutrinos where U is the mixing angle, and where α ∈ {e, µ, τ } labels the active neutrino species in the flavor basis.This then induces a transition matrix elements within the weak current, e.g.L int ⊃ U αN Ni γ µ P L ν α J µ .Above the weak scale the mass-mixing portal is relevant in the Wilsonian sense, however below the weak scale the mixing angle accompanies an irrelevant dimension-6 operator in the 4-Fermi effective theory that governs low-energy neutrino phenomenology.Despite being Wilsonian-irrelevant, the mass-mixing portal can still be efficiently probed at low energies because of its strength relative to SM neutrino interactions which proceed through the same dimension-6 contact operators.
There is, however, one unique portal that is dimension-5 and so can come to dominate over SM weak currents at low energies even if it is sub-dominant at high energies.This is the so-called "dipole portal" which first received substantial attention in the context of the Mini-BooNE and LSND anomalies [29][30][31][32][33].The authors of [34] further pointed out interesting "double-bang" phenomenology that could be probed in experiments such as IceCube.Ref. [35] initiated a broad study of the relevant parameter space for a dipole portal, and this was recently complemented by a thorough analysis of low-energy and cosmological phenomena [36].Since these early studies, the viable parameter space for a neutrino dipole portal has received considerable attention [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54], and has persisted as a potential explanation of the MiniBooNE excess [40,55].
The interaction Lagrangian for the dipole portal is conventionally taken to be ( where d α is the (flavor dependent) transition dipole moment between ν α and the singlet fermion N .In complete generality one could consider a linear combination of magnetic and electric transition dipole portals (see [56] for a recent discussion and [57,58] for related work in the context of angular distributions in HNL decays).It suffices, however, to consider only the magnetic dipole portal as a simplified model in the majority of parameter space, and we restrict our attention to this case here.
Constructing UV-completions that yield sizeable dipole operators is a non-trivial model building task.One constraint stems from neutrino masses, since loop diagrams involving a photon insertion on the incoming and outgoing neutrino can alter neutrino textures.This can be avoided if N is a Dirac or psuedo-Dirac fermion [35].
Ref. [36] discusses possible UV completions connected to leptoquarks and recent B-anomalies, and other models have been discussed in [59].In this paper we work purely at the level of low-energy Lagrangian Eq. ( 2) and remain agnostic to the UV-origin of the dipole portal.For the results show below, we consider N to be a Dirac fermion.From a phenomenological standpoint, the primary difference between Dirac and Majorana HNLs is that the Majorana decay length is half of the Dirac decay length [60].We will show later that the lower bounds we obtain on the model are independent of decay length, so we expect similar bounds for Majorana HNLs.
Finally, we note that HNLs appear generically in dark sector models and that it is natural to consider models in which there are additional light degrees of freedom.While observational evidence demands that lowmass dark sector particles are weakly coupled to the SM, they needn't be weakly coupled to one another and it is consistent and arguably generic for there to exist complicated dark sector dynamics (e.g.[2]).A simple example is a model with three HNLs and one massive Z ′ that interacts with the HNLs via O(1) couplings, but is secluded from the SM except for small kinetic mixing terms with e.g. the SM photon (e.g.[45,[61][62][63]).New dark sector dynamics can modify e.g.decay lengths and couplings to nuclei, however, while details can change, the same basic phenomenology proceeds: neutrinos scatter on nuclei and produce HNLs, and those HNLs decay inside detectors producing a visible decay signature.This is illustrative of the fact that many neutrino-portal dark sectors may be efficiently probed via terrestrial upscattering provided the dark sector has a long-lived visibly decaying particle in its spectrum.
In this work we derive new constraints on neutrino portals using existing data by leveraging atmospheric neutrino upscattering (see Fig. 1 for a pictorial description), with our results summarized in Fig. 2. Recent work [60,64] has identified cosmic ray showers as a potentially useful source of HNLs, however our search strategy differs in that the HNLs are produced via volumetric upscattering within the Earth, rather than being produced directly via meson decays in a cosmic rays shower.Because of the large volume of the Earth, this search strategy is ideally suited for regions of parameter space in which the HNL decay length is smaller than, or comparable to, the radius of the Earth.The signature of interest is the through-going decay of an HNL into some visible SM degrees of freedom and the production mechanism is νA → N X with A some SM particle (typically a nucleus) that is naturally abundant within the Earth's mantle and/or core.The treatment of atmospheric neutrino upscattering is considerably more complicated than solar neutrino upscattering which was pursued previously by one of us in [60,64].The largest technical challenge is that atmospheric neutrinos oscillate over O(km) length scales.This demands a detailed treatment that includes electron number density profiles along arbitrary line segments within the Earth.In this paper we develop a Monte Carlo routine that is capable of computing the expected event yield inside a detector taking into account all relevant physical details.
The rest of this paper is organized as follows.In Section II we discuss the upscattering of atmospheric neutrinos.This includes a discussion of neutrino oscillations, atmospheric neutrino intensities, and relevant formulae for upscattering cross sections (both coherent and incoherent).Details of the numerical implementation are deferred to Appendix A. In Section II C we discuss the visible signatures of through-going HNLs in large volume detectors.For the dipole portal the signal is always a broad spectrum of photons, whereas for the mass-mixing portal branching ratios vary depending on the HNL mass.Next in Section III we derive new constraints on neutrino-portal couplings to HNLs using Super-Kamiokande (SK) data and discuss potential improvements with both Hyper-Kamiokande (HK) and DUNE in Section IV.Finally in Section V we summarize our findings and discuss potential future directions.

II. UPSCATTERING IN THE EARTH
Atmospheric neutrinos supply a broadband spectrum of electron and muon flavor neutrinos ranging from ∼ 100 MeV up to ∼100s of GeV.At these energies the neutrinos pass through the Earth without scattering, but do undergo substantial flavor oscillations that depend nontrivially on the matter profile encountered by the neutrinos in transit.At a typical point inside the Earth, this results in a quasi-isotropic intensity of neutrinos with a O(1) contributions from ν e , ν µ , and ν τ with a broad range of energies as described above.In what follows we outline how to formalize the problem of atmospheric neutrino upscattering νA → N X, with A a SM nucleus and X some SM final state particles.
For our upscattering formalism, we begin with the incoming flux of atmospheric neutrinos.The flux of these neutrinos is sensitive to the neutrino energy and the zenith angle relative to the neutrino entry point. 2 For flavor dependent couplings we include neutrino oscillations, which are affected by the matter profile between the entry position, W, and the interaction point X.The result is an angle and energy dependent neutrino intensity I να (E ν , Ω ν , X) that depends on the neutrino flavor α, and the position inside the Earth, X.This intensity can be related to the standard atmospheric neutrino intensity of flavor β, I ν β , (at the surface) via where a sum over β is implied, and cos ϕ zen is chosen so that the neutrino points from W to X (see Fig. 1).Neutrino oscillation probabilities, denoted by P αβ for ν β → ν α , depend on the position inside the Earth and the angle of incidence since these two parameters determine the neutrino's path through the Earth.The zenith angle at which the neutrino is produced also depends on both X and Ω ν .The neutrino oscillations must be computed separately for each angle, energy, and point inside the Earth.In what follows we take recent best fit values from the NuFit collaboration [71]  .The effect of varying neutrino oscillation parameters within their allowed range of values produces only a small effect.This is because the flux of HNLs arriving at the detector depends on the volume-averaged oscillation probabilities weighted by the broadband atmospheric flux.The result therefore samples a wide range of L/E ν and is relatively insensitive to e.g.δ CP and the mass mixing hierarchy.
While most of the atmospheric neutrino flux has energies of 100s of MeV, the flux extends to 10s of TeV.Momentum transfers can then be much larger than the scale of nuclear coherence Q coh ∼ 100 MeV such that scattering will not be entirely coherent.Instead, the cross section be composed of a coherent and incoherent piece dσ = dσ coh + dσ in.where dσ coh ∼ O(Z 2 , Q 2 w ) and dσ in.∼ O(Z, Q w ) with Z an Q w the electric and weak charge of the nucleus.The former is relevant for the dipole portal and the latter is relevant for the mass mixing portal.The coherent contribution can be reliably treated working in the infinite mass limit because nuclear form factors ensure that |Q| ∼ 0.3 GeV whereas nuclei are generically heavy M A ∼ 30 GeV such that nuclear recoil energies are small In general the upscattering cross section will depend on both the scattering angle and HNL energy d 2 σ/d cos ΘdE N .We are interested in angles, Θ, such that the resulting HNL is directed toward the detector.When considering coherent scattering, the nucleus' recoil can be neglected such that For elastic scattering on free nucleons, recoil effects must be included and the delta function instead relates E N as a function of both E ν and cos Θ.In this work we include the contribution from coherent scattering and incoherent scattering on the constituent nucleons.We model the incoherent contribution as if the scattering took place on free nucleons and neglect detailed nuclear effects.In models where the upscattering is dominated by coherent scattering (e.g. the dipole portal) the nuclear uncertainties are drastically reduced.For the mass-mixing portal .existing limits derived in the literature (for a helpful compilation of projected sensitivities see [42]).We have included two bounds from Super-Kamiokande for the dipole portal; one where σsys = 0 and another where σsys is 5% of the background events.Notice that due to atmospheric neutrino oscillations averaged over the interior of the Earth our constraints are flavor independent up to an O(1) factor.Constraints from NOMAD (µ-only) are taken from [34,65], MiniBooNE (µ-only), supernova bounds (assumed flavor democratic), and LSND constraints are taken from [35], and solar constraints are taken from [60].The τ -only mass-mixing constraints are taken from [47,64,[66][67][68][69].The dashed lines in (d) are meant to illustrate the theoretical uncertainty in the production rate due to incoherent scattering on nuclei.The lower (upper) dashed line is where we double (halve) the contribution due to incoherent scattering.The coherent contribution (which is nearly free of nuclear uncertainties) guarantees an irreducible flux atmospheric upscattered HNLs.
we find that incoherent scattering provides an O(1) contribution to the total rate, and treat the theory uncertainty from nuclear effects conservatively (see Fig. 2d).
The HNL created in this interaction is unstable, with a decay length of where τ is the characteristic decay time in the rest frame of the HNL.Given the decay length, λ, the probability for an HNL produced at a location X directed toward our detector located at Y to decay visibly within the fiducial volume, of length ℓ and area A ⊥ is where B vis = Γ N →vis /Γ N is the branching ratio to visible SM decay products (an experiment and search strategy dependent quantity).The probability of an HNL being directed towards the detector is proportional to the solid angle of the detector as seen from the point of emission, Ω det ∼ A ⊥ /|X − Y| 2 .We note that in the limit λ ≫ ℓ we may approximate 1 − e −ℓ/λ ≈ ℓ/λ such that the over-all rate scales as the volume of the detector V = A ⊥ ℓ.
For λ ∼ ℓ the rate will depend somewhat on the geometry of the detector and so the ceiling of our constraints will be modified by an O(1) number; because this region is already ruled out by complementary search strategies this is not important, except perhaps for electron-only coupled dipole portals (see Fig. 2a).Putting everything together the event rate for HNLs produced by a ν α neutrino portal (with α ∈ {e, µ, τ }) to decay visibly inside the fiducial volume of a detector is given by The term in square brackets is the differential flux of HNLs per unit volume produced at location X and the factor P vis weights the spectrum by the probability of decay within the detector.The event spectrum for a given visible decay product can be found by folding the differential rate dR/dE N computed using Eq. ( 7) with the spectrum of daughter particles produced in the lab frame by an HNL with energy E N .Equation ( 7) cannot be calculated analytically for simple matter profiles due to the complex dependence of the oscillated neutrino intensity as a function of X.Even without oscillations, a realistic density and composition profile of the Earth demands a numerical treatment.We have developed a purpose built Monte Carlo program capable of solving Eq. ( 7) efficiently using conditional importance sampling.The details of our implementation are discussed in Appendix A, however we briefly sketch the procedure here.First, we generate an ensemble of neutrino energies by importance sampling an approximate atmospheric neutrino flux curve.For each neutrino energy we calculate the maximum HNL decay length, which corresponds to the case when We then sample a position inside the Earth, X, from an exponential distribution defined relative to the detector.At each point X the density and composition of the Earth is computed.For each production mechanism (e.g.coherent v.s.incoherent scattering off 56 Fe and 16 O), we generate a random initial neutrino angle (defined relative to Y − X) using a nonuniform sampling that accounts for correlations induced by the differential cross section 3 dσ i /d cos Θ.Given the incident neutrino angle, we then propagate backwards to the point on the Earth's surface where the neutrino would have originated, W, and we calculate the zenith angle relative to the Earth's tangent at that point.The neutrino intensity is then calculated using NuFlux [74] with E ν taken from the first step, and at the required zenith angle to propagate from W to X.All events are saved in an event record, with appropriate weights ac- 3 For highly forward scattering these correlations can make certain numerical methods (e.g.VEGAS) highly inefficient [72,73].By "working backwards" from the detector to the source of neutrinos we efficiently account for these correlations.
counting for the various terms in Eq. ( 7).Finally, in a post-processing stage we calculate the relative weights for the various neutrino flavors at the location X by numerically solving the Schrödinger equation along the line segment connecting the upscattering location X to the neutrino point of origin at the Earth's surface.This is done self-consistently using the same density and composition profile as was used to generate the upscattering events.
We now specialize our discussion to the relevant neutrino portals discussed herein.

A. Dipole Portal
For the dipole portal, the visible decay signal is N → γν with a branching ratio of BR = 1.The analysis is consequently straightforward.The decay length of the HNL is given by which is comparable to the size of the Earth for much of the parameter space of interest.For a neutrino dipole portal, coherent scattering on nuclei is the dominant upscattering process for all energies.The upscattering cross section assumes a simple form when one works in the large mass limit and drops all mass-suppressed effects including nuclear recoil corrections, and contributions from a nuclear magnetic dipole moment.This is a valid approximation because the nucleus' charge form factor ensures that Q ≪ 300 MeV, such that recoil corrections are always small.The resulting differential cross section is of the form The angular dependence dσ/d cos Θ can be obtained by a simple change of variables.The charge form factor of each nucleus is modelled as a Helm form factor with parameters fitted to the tabulated two parameter Fermi distributions from [75] (see Appendix B 3 for more discussion).
Although it is subdominant to the coherent contribution, we also include an incoherent sum over nucleons in our model of upscattering.This cross section is given by dσ A = Zdσ p + (A − Z)dσ n with the proton and neutron cross sections parameterized in terms of standard Dirac and Pauli form factors.
We can see in Fig. 2 that our dipole coupling bounds have a relatively flat region in d when m N is between 0.01 and 0.08 GeV.In Fig. 3, we see that this flat region corresponds to decay lengths satisfying the hierarchy ℓ ≪ λ ≪ R ⊕ .This bound can be estimated through a relatively simple approximation: treat the Earth as being of a constant density composed of a single element, consider the neutrino flux as isotropic, ignore angular dependence on the cross section, consider elastic scattering such that E ν = E N , and set all terms of O(ℓ/λ) and O(λ/R ⊕ ) to zero.Within this approximation, we find where f vis is the fraction of HNL decays in our detector that are in the visible energy range.Here, the only depenedence on the transition dipole moment appears in the cross section as a d 2 (the rate is independent of the decay length for ℓ ≪ λ ≪ R ⊕ ).We can define σ = σ/d 2 , and then our estimate for the floor of our constraint is where d(m N ) is an estimate for the "floor" of our constraint as a function of m N and R exp N →vis is the rate of FIG. 4. Flavor independent curves corresponding to 197 HNL events at Super-Kamiokande with 5326 days of data using the approximation from Eq. ( 11) and the Monte Carlo simulation.
In the approximation, we consider the Earth as composed entirely of silicon with a density of 5 g/cm 3 .
visible energy deposition that can be excluded by the experiment under consideration.In Fig. 4, we see that this approximation closely matches the true bounds that we get for the full Monte Carlo, meaning that this approximation can be used to see how the lower bound will scale with exposure time, volume, decreased background, etc.Since this approximation ignores decay length, it should also be valid for setting approximate bounds on Majorana HNLs, which have half the decay length of Dirac HNLs

B. Mass Mixing Portal
The mass mixing portal is more complicated phenomenologically because production cross sections rise with energy, and for m N ≳ m π many new hadronic decay channels open.We have included many of these details in our simulation, however a posteriori it is clear that searches relying on terrestrial upscattering are only competitive with existing constraints for masses below the pion threshold.We therefore focus our discussion on the case of the decay channel N → e + e − ν which is the only visible decay mode for m N ≤ 135 MeV.
The decay of an HNL to an e + e − pair depends on the flavor structure of the mass mixing portal and the flavor of the invisible SM neutrino.These effects introduce an O(1) prefactor that depends on the final state which can be found in [22,24,76], however the dominant effect is the muon decay like formula for the partial width The result is that HNL decay lengths are extremely long for low masses, and can easily exceed the radius of the Earth by orders of magnitude.In this regime terrestrial upscattering offers substantial benefits over traditionally laboratory based searches, and can offer leading sensitivity on |U τ N | 2 .HNLs can always decay to three neutrinos, N → ννν for any non-zero mixing angle.The result is a branching ratio that is O(10%) for N → νe + e − for all HNL masses below the pion threshold.We include this effect in our simulations computing the full decay length and taking B vis = Γ e + e − ν /Γ.
HNL upscattering proceeds via the weak neutral current and for the relatively low HNL masses that we focus on here all of coherent (i.e.CEvNS), quasi-elastic, and deep inelastic scattering contribute to the upscattering yield.We find that for regions of parameter space where atmospheric upscattering is competitive that the scattering mechanisms are dominantly coherent and incoherent scattering on nucleons, with deep inelastic events contributing only a few percent to the total flux.
The coherent contribution is relatively insensitive to the nuclear species and can be calculated from first principles.We model the weak nuclear form factor by setting it equal to the charge nuclear form factor. Incoherent scattering on nuclei is modelled as described above for the dipole portal.This neglects all effects of nuclear structure and we therefore expect a sizable theoretical error from our modelling.Unlike the dipole portal case, we find that incoherent scattering makes up roughly twothirds the total upscattered flux.Owing to its relative importance, we have included "error bands" in Fig. 2 in which the incoherent scattering cross section has been doubled, and halved respectively; we believe this to be a conservative overestimate of the theoretical uncertainty.

C. Decays inside the Detector
In the presentation above we have outlined how to calculate the flux of unstable particles arriving at a given large volume detector.This flux is not directly visible, and the bona fide observable is the energy and angular distribution of an HNL's visible daughter particles.For illustration, we discuss the case of a dipole portal decay N → νγ in detail below.The case of a three-body decay, as in N → νe + e − is qualitatively similar, but slightly more involved due to the three body final state.The details of the decay distribution do not substantially impact our rate-only estimate, although their details may be relevant for future searches that we outline in Section IV.
In the case of the dipole portal, when the HNL decays, it decays into a photon and neutrino.The angular distribution of a dipole-mediated decay in the HNL rest frame depends on the level of CP violation [56,57,[77][78][79], with dΓ/d cos ζ ′ ∼ 1 + α cos ζ ′ and α ∈ [−1, 1] and ζ ′ the angle between the photon and the HNL polarization.For simplicity, we take α = 0 such that the decays are isotropic.Our sensitivity is only mildly sensitive to this choice; α > 0 leads to a somewhat harder photon spectrum in the lab frame, while α < 0 leads to a somewhat softer spectrum (see e.g.related discussion in [60]).In the rest frame, E ν,rest = E γ,rest = E N /2, and this leads to the following lab frame kinematic variables A flat (i.e.isotropic) distribution in cos ζ ′ results in a "box distribution" for E γ,lab ranging between [E Knowing the initial momentum of the HNL, we can sample cos ζ ′ uniformly between [−1, 1] and generate a random sample of angles of the detected photon relative to the horizon at the detector ϕ det .

III. SUPER-KAMIOKANDE CONSTRAINTS
We now turn to our analysis of public data from Super-Kamiokande, which when coupled with our Monte Carlo simulation, allows us to set new limits on neutrino portal couplings.SK is a large volume (22.5 × 10 3 m 3 fiducial volume) Cherenkov detector whose primary background is the scattering of atmospheric neutrinos passing through the detector.It is well suited to search for through-going HNL decays and has a large statistical sample of atmospheric neutrino events which can be used to set limits on the rate of visible HNL decay [80,81].
The SK collaboration classifies events as sub-GeV (30 MeV < E vis < 1.33 GeV) and multi-GeV (E vis > 1.33 GeV) with sub-classifications for each event type.In the sub-GeV sample events are classified as e-like, µ-like or π 0 -like, single-ring or two-ring, and 0 decay-e, 1 decaye, or 2 decay-e.The decay-e classification is meant to capture Michel-electrons from muon decay, while the particle identification (PID) is based on characteristic Cherenkov ring patterns of each particle.The multi-GeV sample is split into partially contained and fully contained, the former applying exclusively to muon events.In the fully contained sample events are classified as single-ring or multi-ring and are then further sub-divided as ν e -like, νe -like, or µ-like.The ν e -like vs νe -like samples are defined by a cut on the number of decay-e events which ultimately stem from a ν e n → e − π + n interaction with subsequent pion decay at rest, followed by muon decay [82].Not all ν e interactions produce a π + and so there is substantial cross-contamination between the two samples; by way of contrast the µ-like sample is relatively pure.
In what follows we describe a simple rate-only analysis based on the published results in [82].For each model we focus on a the relevant experimental signature and use the experimental collaboration's Monte Carlo prediction as the expected Poisson mean of the event sample.Given their observed data, we then set limits at the 95% confidence level on the number of allowed events in the energy range as defined by the experiment.We consider both systematics and statistically limited searches with a conservative estimate of a 5% systematic uncertainty on the collaboration's Monte Carlo prediction for their sub-GeV sample of 0-decay-e events.

A. Dipole Portal
The dipole portal's only signature is a single photon which will be classified as an e-like 0 decay-e signature in the sub-GeV analysis and as a fully-contained νe -like event in the multi-GeV sample.The energy distribution of photons is broad for all HNL masses but its precise shape depends on both m N and d.The multi-GeV and sub-GeV samples therefore provide complementary tools with which to probe the HNL parameter space.
We set limits by taking the union of the excluded regions from the multi-GeV and sub-GeV analyses separately and these are shown in Fig. 5.We note that the constraints cross around d ≃ 5 × 10 −9 MeV −1 , with the multi-GeV search dominating for larger d and the sub-GeV dominating below.We consider d e = d µ = d τ = d (flavor independent), and flavor dependent couplings accounting for neutrino oscillations in each case.Based on Table II of [82], we assume a Poisson mean from the collaboration's Monte Carlo simulation of µ MC = 10266 sub-GeV 0 decay-e events while the observed count is N obs = 10294.
It is tempting, given the close agreement between Monte Carlo and observation to infer that systematic uncertainties on the Monte Carlo prediction are fully under control, however the Super-Kamiokande Monte Carlo is tuned to their data to self-consistently determine, e.g. the atmospheric flux normalization.In the presence of new physics this tuning could be compromised and so it is important to estimate a systematic uncertainty on an experiment such as Super-Kamiokande.First, note that while the overall normalization of the atmospheric flux is poorly constrained, the ν e : ν µ ratio is known to within a few-percent [70,82].Therefore the flux normalization at Super-Kamiokande can be fixed using muonexclusive sub-samples, and the electron flux can be subsequently inferred.Second, it is worth noting that the sub-GeV 0-decay e bin of the Super-Kamiokande data set is relatively insensitive to neutrino oscillations, and its background modelling is therefore reasonably robust.Quantitatively, one can compare the predicted flux with and without oscillations from Fig. 14.4 of [83]; the flux changes by only 2.7%.We therefore conclude that a 5% systematic uncertainty can be conservatively applied to the Super-Kamiokande Monte Carlo prediction of the 0decay e sub-GeV background from atmospheric ν e scattering.
For finding constraints, we take the statistical uncertainty as σ stat = √ µ M C + µ HN L .We take a conservative upper bound on the systematic uncertainty at σ sys = 0.05µ M C .We then solve P (x ≤ N obs |µ, σ) = 0.05 where our probability distribution function is (2πσ) −1/2 exp (x−µ) 2 2σ and µ = µ MC + µ HNL and σ = σ 2 stat + σ 2 sys .For σ sys = 0, we find that µ HNL = 197 is excluded at 95%-CL; this corresponds to the number of events per 328 kt-yr (corresponding to 5326 live-days at Super-K).For σ sys = 0.05µ MC , we find that our 95%-CL bound now corresponds to µ HNL = 893.For the multi-GeV analysis we take the νe -like sample which has µ MC = 2194 and N obs = 2142 for a 328 ktyr exposure.Following the same procedure as above, the σ sys = 0 95%-CL bound is µ HNL = 26 and the σ sys = 0.05µ MC 95%-CL bound is µ HNL = 145.It is worth noting that the excluded number of excess events is determined solely by uncertainties in the rate of Standard Model events at Super-K.Uncertainties in our BSM theory only affect how we translate these exclusion bounds into parameter space.We also expect that any uncertainties in our theory are accounted for by our conservative estimate on the Super-K uncertainty.
Using our Monte Carlo integrator, we compute the rate of HNLs passing through and decaying within the detector.The photon spectrum is generated using the lab frame decay distribution of the HNL.For flavor dependent dipole couplings we re-weight the ensemble of Monte Carlo events by adjusting the intensity of neutrinos at the upscattering location according to the oscillation probabilities computed along the line segment connecting the upscatter location to the position on the Earth's surface above which the atmospheric neutrino is produced.Photon detection efficiencies are taken to be unity, ϵ γ = 1 which we believe to be reasonable as Super-K can reach 100% trigger efficiency on events with 4.49 MeV of energy [84], and our photon energies are well above this.The photon spectrum is integrated from E γ = 30 MeV to E γ = 1.33 GeV for the sub-GeV sample and from E γ = 1.33 GeV to the highest energy in the Monte Carlo sample for the multi-GeV sample.

B. Mass Mixing Portal
For the mass mixing portal we focus our analysis on N → e + e − ν which is the only visible decay mode for m N ≲ 130 MeV, and contributes for all HNL masses.As we will show the only region in which terrestrial upscattering can compete with fixed target experiments is in the low mass regime and so this suffices for our purposes.
An e + e − pair will appear as highly collimated and result in an electromagnetic shower that is difficult (or impossible) to distinguish from a single electron or single photon.The HNLs are sufficiently boosted such that wide-angle e + e − pairs are a non-issue and the decay signature maps onto the same search channels as the single FIG. 5. Exclusion contours assuming a statistically limited search (i.e.σsys = 0) from Fig. 2 (a), (b), and (c) and the equivalent constraint for a flavor independent dipole portal.Our constraints are only moderately sensitive to neutrino flavor due to substantial oscillations within the interior of the Earth.Point A is representative of parameter space that is dominated by the sub-GeV sample at Super-K.Point B is representative of a point in parameter space that is dominated by the multi-GeV sample.We discuss this further in Section IV and show distributions for Point A in Figs. 6 and 8 and Point B in Figs.7 and 9.
photon analysis.We can therefore take the rate-only exclusions from above an apply them directly to the mass mixing portal.The sub-GeV sample provides the best sensitivity to HNL mass-mixing over the full range of parameter space and we find that new regions of parameter space for τ -coupled HNLs can be probed with existing Super-Kamiokande data.
For our upscattering simulation we include coherent, quasi-elastic, and deep inelastic scattering channels.We do not include resonance production, nor do we account for nuclear structure (e.g.Pauli blocking, giant dipole resonances etc. ).We note that in our region of sensitivity, incoherent scattering off of nucleons is the dominant contribution (contributing to around 2/3 of the rate).Coherent scattering contributes to around 25 percent of the total rate, and DIS contributes to less than 10 percent of the rate.

C. Scaling with increased sensitivity
Before moving on to considering future experiments, let us discuss how the constraints described above scale with increased sensitivity.Importantly, the right-most boundary of our exclusions corresponding to an upper bound on m N is set by our sensitivity rather than by any kinematic thresholds.This is because the flux of atmospheric neutrinos is broad and there is no fundamental limitation on the mass of HNLs which can be produced.As sensitivity improves, either by collecting more data, improving background discrimination, or by leveraging new detector technologies, heavier HNLs can be probed.Furthermore, in the regions of parameter space highlighted in Fig. 3, limits on the dipole coupling scale as d ∼ (sensitivity) 1/2 which is extremely advantageous relative to the naive scaling of d ∼ (sensitivity) 1/4 that one would expect in the long-lifetime limit.Taken together, this suggests that improved sensitivity using atmospheric upscattering as a source of HNLs has a high return on investment.

IV. SEARCH STRATEGIES AT FUTURE EXPERIMENTS
In this section we discuss how to improve future searches for HNLs.All of the exclusion contours in this work are based exclusively on the simple rate-only estimates of the previous section and the reader may view the following discussion as an outlook towards future improvements.Our Monte Carlo routine can generate kinematic distributions such as the energy and zenith angle of the HNL's decay products, and these can be leveraged to improve signal to background ratios.We also discuss potential improvements in background rejection using different detector technology (e.g.liquid scintillator and/or liquid argon time projection chambers instead of a Cherenkov detector), and the impact of a larger fiducial volume in a detector such as Hyper-Kamiokande (HK).
Let us begin by discussing improvements that can be had by taking into account the energy and angular distributions of the observed photons.For illustration, let us examine two dipole model parameter points which are at the Super-Kamiokande exclusion boundary: the two points marked by triangles in Fig. 5. Recall that these searches exclude a window of couplings, originating from the requirement that decay lengths satisfy: 10 m ≲ λ ≲ R ⊕ .Consequently we will refer to the "floor" and "ceiling" of the coupling exclusion region.
Let us first examine the angular and energy distributions of the detected photons from a parameter point on the floor of the exclusion region, m N = 0.03 GeV and d e = 2 × 10 −10 MeV −1 .We include an example of these distributions on our lower bound for electron-flavor coupling Fig. 8.We see the angular distribution highly favors angles less than π/2, corresponding to upward going photons in our case.Since the photon direction is highly correlated with the HNL direction, this means most of our signal comes from HNLs produced below the detector.Our decay length is large for these parameters, so there is far more volume available for scattering below the detector than above it.The energy distribution is peaked at lower energies both as a consequence of the atmospheric neutrino flux, which falls off quickly with energy, and because of the skewed distribution of photons from the decay of relativistic HNLs [see discussion near Eq. ( 15)].The most probable photon energy in the lab frame is given by one half the HNL's energy.At low FIG.6. Angular distribution of the detected photons for Point A in Fig. 5 (mN = 0.03 GeV and d = 2 × 10 −10 MeV −1 ).At low HNL masses and small dipole couplings, the angular distribution of photons is primarily up-going (cos ϕ det < 0).This can be used in future searches to cut backgrounds.
HNL masses, almost all of the atmospheric neutrino flux is capable of producing HNLs.For coherent scattering E N = E ν , and since the atmospheric neutrino flux falls like a power-law the HNL inherits this feature resulting in lower energy photons.
The angular distribution of observed events at Super-Kamiokande is nearly uniform [80].Therefore, we can choose to only look at upward-going events and cut our background in half while keeping our signal virtually unchanged.We expect that this will extend our lower bound of d by a factor of 2 1/4 , as we will need 1/2 as many events to reach the same level of uncertainty, and the number of upscattering events goes as d 2 .Now let us turn our attention to a parameter point on the ceiling of the excluded region.In this case, the energy and angular distributions are qualitatively different.
We show examples of the electron flavor case for m N = 0.1 GeV and d e = 2 × 10 −9 MeV −1 (Fig. 7 and Fig. 9).We see that the photons are now more uniformly distributed in angle.This is because the decay length is now much shorter than the depth of the detector, so the volume available for upscattering is approximately spherical, and we expect our flux of HNLs at the detector to be roughly uniform in angle.We see that there is a slight peak near cos ϕ det = 0, since more neutrinos come from the horizontal direction than from the vertical.
The energy distribution for HNLs with shorter decay lengths (corresponding to large masses and strong couplings) is nearly uniform in the sub-GeV sample at Super-Kamiokande.This can be understood as follows: the flux of HNLs roughly mimics the flux of atmospheric neutrinos and so in the sub-GeV regime is relatively flat.The flux of photons is given roughly by dN/dE γ ∼ A ⊥ λP vis (λ, ℓ) we find λP vis ≈ ℓ and the energy distribution is given by dN The box distribution of photons is flat, with a height that scales as ∼ 1/E N such that the difference between the height of two bins of the histogram is given by For λ ≫ ℓ, as in Fig. 8, P vis ≈ ℓ/λ and the overall rate is independent of λ as discussed above.For λ ∼ ℓ, however, the probability of decaying inside the detector becomes some O(1) number P vis ≈ 1 − exp[−ℓ/λ] and the photon spectrum becomes proportional to λ.Since λ ∝ E N ∝ E γ , this cancels against the 1/E γ denominator of Eq. ( 15), and the spectrum for shorter decay lengths inherits the shape of Φ N which is relatively flat for sub-GeV energies; this explains the shape of Fig. 9.While our code is able to estimate the distribution of photons within Super-K, we only provide results from a rate-only analysis.A full analysis including the angular and spectral distributions of the decay photons would require considerations such as energetic and angular reconstruction of photon events in the detector, which are beyond the scope of this paper.We encourage those who have strong familiarity with Super-K to perform a more detailed analysis to improve constraints on new physics couplings.
We now consider future prospects and specifically upcoming experiments with larger fiducial volumes and stronger background rejection methods.We consider Hyper-K, DUNE, and JUNO, all with 10 years of live time.We assume that the Atm Rate Fiducial volume for all three of these experiments is the same as Super-K.We also assume that the energy range of interest will be the same FIG.8. Energy distribution of the detected photons for Point A in Fig. 5 (mN = 0.03 GeV and d = 2 × 10 −10 MeV −1 ).At low masses and small dipole couplings, the energy spectrum is IR peaked.The binning in Eγ is chosen to correspond to the binning in Super-Kamiokande's analysis of the sub-GeV atmospheric neutrino event sample [82].
as Super-Kamiokande (30 MeV − 1.33 GeV).When possible, we will consider an angular cut, meaning that we will only look at upward going events.For Super-K, the Monte Carlo predicts 51% of the e-like 0-decay-e events to be down-going, while 49% are up-going [81].We assume this holds for all experiments.
For ten years at Hyper-K, we expect roughly 70,000 sub-GeV e-like 0 decay-e atmospheric events.Super-Kamiokande run-IV is already doped with gadolinium and data from this exposure will have lower backgrounds from atmospheric neutrinos due to high-efficiency neutron tagging [85,86].In our projections, we assume that Hyper-Kamiokande will be doped with gadolinium, which will cut the background from atmospheric neutrinos roughly in half.An angular cut will let us cut another ∼ 50% of the background, leaving us with roughly 17,000 background events.Assuming no systematic uncertainty, a 95%-CL bound can be set with 218 HNL events.This is, however, not likely to be realistic.The background uncertainties at Hyper-Kamiokande suffer from the same issues as at Super-Kamiokande where statistics are already high enough that a 5% systematic uncertainty makes the search entirely systematics limited.The increased statistical sample will not be helpful unless the systematic uncertainty on the Monte Carlo prediction can be brought down to sub-percent levels even after accounting for reduced background rates from neutron tagging.
We expect roughly 9000 atmospheric events in a "sub-GeV sample" at DUNE over a ten year run-time.Since DUNE is a liquid-argon time projection chamber, it will be easier to distinguish the HNL decay products from a neutrino interaction.In particular, LArTPC technology offers: i) the ability to statistically discriminate ).At larger masses and stronger dipole couplings, the energy spectrum is relatively flat.The binning in Eγ is chosen to correspond to the binning in Super-Kamiokande's analysis of the sub-GeV atmospheric neutrino event sample [82].
between electrons and photons using measurements of dE/dx at the beginning of a track, ii) MeV-scale reconstruction capabilities that can tag gamma rays from nuclear de-excitations [87], and iii) the ability to measure final state charged hadrons including protons, pions with kinetic energies above ∼ 10 MeV [88].Finally, recent work has demonstrated that neutron tagging, using "sprays", may also be possible [89].Importantly, for our background tagging purposes we only need to veto nuclear scattering events and/or single electron showers which is a much easier task than the energy reconstruction considered in [89].While the ultimate capabilities of DUNE will require detailed simulation, we estimate that requiring upward going events will cut 50% of the background, that proton tagging will catch 80% of the remaining events, and that searches for neutron sprays and associated gamma rays from nuclear de-excitation can cut out 80% of the remaining events for which no final state proton is produced.A naive combination of these estimates then suggests that 98% of the background could be rejected at DUNE, however as we have already mentioned above, the precise value will require dedicated simulations.Performing an angular cut as well (requiring an up-going shower), we would expect a background of 88 events.With zero systematic uncertainty, 18 HNL event would set a 95%-CL bound, and this conclusion remains unchanged even if we allow for a ∼ 5% systematic on the background uncertainty.Note that unlike SK, DUNE will be statistically limited provided it can achieve background rejections that are better than ∼ 90%.
Finally, for JUNO, we estimate roughly 6,500 atmospheric events over 10 years.We expect that JUNO will be able to effectively cut out background from events that  [51], and solar DUNE beam upscattering events were derived in [42].We have included on Hyper-Kamiokande and JUNO for when σsys = 0 and when σsys is 5% of the projected background.For our DUNE projection, the background was so low that the two bounds were nearly identical, so we only included the σsys = 0 curve.For the mass-mixing constraints, we use current bounds from [47,64,[66][67][68][69], along with a projected bound for a multi-purpose DUNE near detector from [26].
eject a proton due its low ∼ 1 MeV detection threshold, and will have a 50% efficiency at cutting background events that eject a neutron by leveraging the 2.2 MeV gamma ray from np → dγ [90].Using the relative distributions of neutrinos and anti-neutrinos, we estimate that 20% of the atmospheric events will not produce a free proton and will instead produce a free neutron [91].Assuming a ≳ 95% efficiency at tagging protons, we therefore estimate that the background from atmospheric neutrino CCQE will be ∼ 10% of our re-scaled background from SK i.e. ∼ 650 events.At these levels of statistics JUNO's statistical and assumed systematic uncertainties, taken again as 5%, are comparable.
For our projected limits on |U τ N | 2 , we consider a ten year period at DUNE.As mentioned before, we expect to have significant background reduction by rejecting events that have evidence of scattering off of a nucleon.As a benchmark, we take 30 BSM events in ten years at DUNE to be statistically significant, and we show this contour in Fig. 10d.We see that the projected constraint from a multi-purpose near-detector at DUNE [26] is quite close to constraints stemming from BBN (see Fig. 10d).These two constraints contain much of the parameter space that would be probed by our projected atmospheric upscattering method, and so while we expect improved sensitivity it is likely that the DUNE near detector can cover most of the relevant parameter space.It would be interesting for future studies to identify near-term experiments that may be able to supply new constraints using existing data sets by leveraging atmospheric upscattering.

V. CONCLUSIONS
Atmospheric neutrinos have already given us a wealth of information regarding the nature of neutrino oscillation physics.Here we have seen that atmospheric neutrinos also provide a powerful and qualitatively distinct window into the nature of neutrino interactions and heavy sterile neutrinos.In particular, we have examined the upscattering of atmospheric neutrinos to HNLs and their subsequent decay inside of terrestrial detectors, finding that current data from Super-Kamiokande already yields leading constraints on both the dipole and mass-mixing portals.For the dipole portal, these bounds eat into new parameter space for HNL masses around 10 MeV ≲ m N ≲ 100 MeV, with precise constraints depending on which active neutrino flavor coupling dominates.Similarly, for the mass-mixing portal our Super-Kamiokande constraints provide leading constraints on the tau-sterile mixing angle for HNLs around ∼ 20 MeV.
In the near future, experiments such as DUNE, Hyper-Kamiokande, and JUNO, will be able to take advantage of improved particle identification, background rejection, and employ dedicated search strategies incorporating angular distribution of the events to improve the bounds on the dipole couplings by around a factor of ∼ 2.5 at low HNL masses.As such, this search strategy nicely complements the DUNE beam upscattering [42] and double bang searches [47] which provide better sensitivity to higher HNLs masses, m N ≳ 100 MeV on the dipole portal coupling, and HNLs produced from meson decays at DUNE for the mass-mixing portal [76].As an illustration of the strength of the future bounds, it is striking to observe that JUNO, Hyper-K, and DUNE appear poised to overlap in coupling reach with SN1987A [35], and will therefore close off an allowed gap in couplings.
while for ϵ = 0 we have For sampling the scattering angle, we sample a uniform number χ ∈ [0, 1], such that the angles are sampled as for ϵ = 1 and ϵ = 0 respectively.Finally, we want to sample the neutrino entry position W. We define v in = X − W. Using our scattering angle Θ that we sampled and another angle ψ uniformly sampled from 0 to 2π, then Where Y − X, v 1⊥ , v 2⊥ are all mutually orthogonal.To find the length of the path travelled, we use Having sampled our energy, scattering angle, interaction position, and neutrino entry position, we can calculate other necessary values.From W and X, we calculate the zenith angle ϕ zen of the incoming neutrinos (Fig. 1).We use NuFlux [74] to calculate I incoming (E ν , ϕ zen ).When working with flavor dependent couplings, we calculate a 1D density profile from W to X using the Preliminary Earth Reference Model (PREM) [92].Oscillations are calculated by integrating along this density profile to obtain I να (E ν , ϕ zen , X).
At X, we use the PREM [92] to calculate the density, and call a saved dictionary of the number density for each of the elements [93,94].We calculate the cross sections for the scattering using methods described in Appendix B to get We calculate the decay length of the HNL using Eq. ( 5), and probability for producing a visible decay from Eq. ( 6)

Weighting
Since we preferentially sample our values, we must include a weighting when calculating the rate of decays from the Monte Carlo.These weights are calculated by taking the ratio of our sampling distribution to the true integrand.Explicitly, the weights are given by, where In a Monte Carlo without preferential sampling, we would have standard "Reimann weights", and N the number of samples in our Monte Carlo.Notice that upon combination (multiplying all weights together) the denominators in our new sampling weights, Eq. (A13), cancel against the uniform weights such that we obtain the weight for the i th sample as We now have everything needed to compute the rate This is the numerical cousin of Eq. (7).In this routine, between 20,000 and 300,000 events were generated for each simulation.The larger simulations were necessary for the mass-mixing model.For dipole upscattering, the cross section can be decomposed into where A and Z are the atomic mass and atomic number of the nucleus respectively and F (−t) is the form-factor for the transferred momentum.We work in the infinite mass limit, so our transferred momentum goes as The coherent scattering is given in Eq. ( 9) with Q 2 = −t.Meanwhile, incoherent scattering off of protons or neutrons will go as Here, m p is the mass of a nucleon, and s is the centerof-mass energy given by m 2 p + 2m p E ν .The value of t for incoherent scattering goes as If the neutrino scatters incoherently, then the energy of the HNL is where and µ p = 2.793, µ n = -1.913and Q 2 = t.
In our routine, while we use Eq.(B2) and Eq.(B4) for finding the transferred momentum, we do not know if the true scattering is coherent or incoherent (Eq.(B1) says that we have components of both).Therefore, we let the energy of the propagating HNL be E ν .
When implementing the full cross section, we find that the coherent part still dominates, so most of our phenomenology can be explained by considering the coherent case.

Mass Mixing Routine
In the mass-mixing portal, we must consider coherent elastic scattering, incoherent scattering on nucleons, and deep-inelastic scattering (DIS).Unlike the dipole portal, where nearly all scattering is coherent, the mass-mixing model has significant contributions from incoherent scattering.Since different forms of scattering leads to different HNL energies (and therefore different observed energies), we run the simulations independently for each type of scattering, and then sum together the rates to get the total contribution.
For coherent scattering, we have In Eq. (B8), G F is the Fermi coupling constant and Q w is the weak charge of the nucleus.
Incoherent scattering is modelled by treating nuclei as collections of free nucleons and using standard hadronic form factors for the nuclei.We take dipole parameterizations of the vector, magnetic, and axial form factors and rely on a partially conserved axial current ansatz for the pseudoscalar form factor; explicit expressions can be found in [97,98].
For DIS, we consider scattering off of individual quarks.To find the cross section, we first find cross section for scattering off of quarks σ f as a function of the momentum carried by the quark.This is parameterized by x, the fraction of the total longitudinal nucleon momentum carried by the quark.Finally, we have, where f f (x) is the parton distribution function (PDF) for the particular quark.We numerically performed the integral in Eq. (B9), using PDFs from [99] and treating the HNL as massless, since the masses we are sensitive to are far below the GeV energy scale where DIS becomes important.Although the cross section for DIS is scales linearly with neutrino energy the DIS scattering only contributes on the order of a few percent in the region of parameter space that is not covered by existing searches.Therefore, in our code, we use a simple form for DIS with the leading coefficient determined by Eq. (B9),

(B10)
From these scattering channels, we can determine the number of visible decays (in this case, N → νe + e − ) expected for Super-K.We assume that Super-Kamiokande will be unable to resolve both the e + and the e − , so the decay will appear as a sub-GeV 0 decay-e event.We calculate the energy of the e + e − pair by using the invariant mass distribution in [100], and require that the energy be between 30 MeV and 1.33 GeV.We see the resulting bounds in Fig. 2d.

Form Factor Fitting
We can see that for Eq.(B1), we need a way to calculate the nuclear form factor.The Helm form factor allows us to accomplish this × e −(sQ) 2 /2 , (B11) where j 1 is a spherical Bessel function of the first kind and Rather than using default parameters as a global description of all nuclei, we fit the values of R A , r 0 , and s independently for each nucleus, and then store these values for later calculation.We begin by taking the 2-parameter Fermi distribution for the radial nuclear charge distribution from [75].Taking the 3D Fourier transform of this charge distribution gives us the charge form factor.Using initial values of R A , r 0 , and s, we define the difference between the Fermi and Helm form factors as We then use gradient ascent to iteratively improve our fit (i.e.R 1 → R 1 − (dS/dR 1 )δR 1 where δR 1 is some predetermined constant).

FIG. 1 .
FIG. 1. Schematic of upscattering within the Earth and decaying in the detector.A neutrino enters at W with angle ϕz relative to the vertical.The neutrino scatters into an HNL at X with scattering angle Θ.The HNL reaches the detector at Y, decays within the detector into a neutrino and a visible photon.The photon is emitted with angle ζ relative to the HNL, and detected with angle ϕ det relative to the detector.

( a )FIG. 2 .
FIG.2.Comparison of our dipole-portal limits (a)-(c) and mass-mixing limit for |UτN | 2 (d) vs. existing limits derived in the literature (for a helpful compilation of projected sensitivities see[42]).We have included two bounds from Super-Kamiokande for the dipole portal; one where σsys = 0 and another where σsys is 5% of the background events.Notice that due to atmospheric neutrino oscillations averaged over the interior of the Earth our constraints are flavor independent up to an O(1) factor.Constraints from NOMAD (µ-only) are taken from[34,65], MiniBooNE (µ-only), supernova bounds (assumed flavor democratic), and LSND constraints are taken from[35], and solar constraints are taken from[60].The τ -only mass-mixing constraints are taken from[47,64,[66][67][68][69].The dashed lines in (d) are meant to illustrate the theoretical uncertainty in the production rate due to incoherent scattering on nuclei.The lower (upper) dashed line is where we double (halve) the contribution due to incoherent scattering.The coherent contribution (which is nearly free of nuclear uncertainties) guarantees an irreducible flux atmospheric upscattered HNLs.
7 n F I y V u p m c J 2 y 4 B z n s D U B K w J 7 P D O A z 7 J I T e x f I b e I p i e W k m 8 m S P L H H t A T e n b u n U f n x X m d R 5 e c x c w B + g H n 7 R P K + 5 s h < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " x A B H G c P f M 2 b V 7 5 B P a e m E T u 1 a N A I = " > A A A C E n i c d V D N S g M x G M z 6 b / 2 r e v Q S L E J 7 W X a 3 t d W D U P T i s Y J V o V v K t 9 m 0 h i a 7 S 5 I V y l J f w Y u v 4 s W D I l 4 9 e f N t z G o F F R 0 I G W b m I / k m S D h T 2 n H e r K n p m d m 5 + Y X F w t L y y u p a c X 3 j T M W p J L R N Y h 7 L i w A U 5 S y i b c 0 0 p x e J p C A FIG. 7. Angular distribution of the detected photons for Point B in Fig. 5 (mN = 0.1 GeV and d = 2 × 10 −9 MeV −1 ).The angular spectrum is relatively flat, because at shorter decay lengths, the mountain above Super-Kamiokande contributes an O(1) fraction of the upscattering events.

FIG. 9 .
FIG.9.Energy distribution of the detected photons for Point B in Fig.5(mN = 0.1 GeV and d = 2 × 10 −9 MeV −1 ).At larger masses and stronger dipole couplings, the energy spectrum is relatively flat.The binning in Eγ is chosen to correspond to the binning in Super-Kamiokande's analysis of the sub-GeV atmospheric neutrino event sample[82].

FIG. 10 .
FIG.10.Comparison of our projections for future dipole-portal and mass-mixing-portal limits vs. other projections derived in the literature.Constraints from FLArE 100 and FASERν2 were derived in[51], and solar DUNE beam upscattering events were derived in[42].We have included on Hyper-Kamiokande and JUNO for when σsys = 0 and when σsys is 5% of the projected background.For our DUNE projection, the background was so low that the two bounds were nearly identical, so we only included the σsys = 0 curve.For the mass-mixing constraints, we use current bounds from[47,64,[66][67][68][69], along with a projected bound for a multi-purpose DUNE near detector from[26].