Searching for Physics Beyond the Standard Model in an Off-Axis DUNE Near Detector

Next generation neutrino oscillation experiments like DUNE and T2HK are multi-purpose observatories, with a rich physics program beyond oscillation measurements. A special role is played by their near detector facilities, which are particularly well-suited to search for weakly coupled dark sector particles produced in the primary target. In this paper, we demonstrate this by estimating the sensitivity of the DUNE near detectors to the scattering of sub-GeV DM particles and to the decay of sub-GeV sterile neutrinos ("heavy neutral leptons"). We discuss in particular the importance of the DUNE-PRISM design, which allows some of the near detectors to be moved away from the beam axis. At such off-axis locations, the signal-to-background ratio improves for many new physics searches. We find that this leads to a dramatic boost in the sensitivity to boosted DM particles interacting mainly with hadrons, while for boosted DM interacting with leptons, data taken on-axis leads to marginally stronger exclusion limits. Searches for heavy neutral leptons perform equally well in both configurations.


I Introduction
The near detectors of long-baseline neutrino experiments, once considered an afterthought to reduce systematic uncertainties in oscillation measurements, are nowadays independent experiments in their own right. Besides delivering a wealth of data on neutrino interaction physics, it has been realized that they could also serve to probe the existence of physics beyond the Standard Model (SM) [1] such as heavy neutral leptons or other new particles, possibly connected to the dark matter (DM) puzzle [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. Indeed, signatures of new dark particles at these experiments can be linked in a predictive way to compelling scenarios of light DM. For instance, "invisible" decays of new sub-GeV particles that mediate light DM-SM interactions, can be searched for at neutrino fixed target facilities by looking for scattering of the decay products off nucleons and/or electrons in the near detector. A light DM program at neutrino facilities could complement the next generation light DM direct detection program [21], in particular, the upcoming experiment SENSEI [22]. So far, only the Fermilab-based neutrino experiment MiniBooNE has performed dedicated searches for light DM [23,24] and very recently MicroBooNE has released the first search for heavy neutral leptons (HNL) [18] and for a light Higgs decaying into e + e − [25], but the untapped potential is big, with many past and ongoing experiments having the capability to supersede MiniBooNE's sensitivity [8,9,11,15]. Importantly, these searches can typically be done fully parasitically to the main neutrino program [9,11].
The near detector physics program will be taken to the next level by the DUNE-PRISM detectors, to be installed 574 m downstream from the target [26] at the long-baseline neutrino facility (LBNF) at Fermilab, the neutrino source for the DUNE experiment. These detectors -a liquid argon time projection chamber (TPC) and a magnetized gaseous argon TPC -will be mounted on a movable platform, allowing them to be displaced up to 30.5 m (53 mrad) away from the beam axis. This capability mainly serves the detectors' primary purpose, namely constraining the unoscillated neutrino flux and measuring the neutrino cross sections. In particular, the neutrino spectrum changes as a function of the off-axis angle (for kinematic reasons), while the cross sections obviously do not. Therefore, taking data at different off-axis positions will allow DUNE-PRISM to disentangle the uncertainties in the neutrino spectrum from the uncertainties in the neutrino cross section.
In this paper, we will discuss the impact of the DUNE-PRISM concept on searches for physics beyond the SM. More specifically, we will consider the production of light ( GeV) and very weakly interacting new particles in the target, followed by their interaction or decay inside the DUNE-PRISM detectors.
Among the numerous extensions of the SM that can be probed in DUNE-PRISM and other accelerator neutrino experiments, we will consider in particular: (1) light ( GeV) DM particles produced with a large Lorentz boost and detectable via dark photon-mediated DM-electron scattering; (2) light DM particles detectable via DM-nucleus scattering (leptophobic DM); and (3) heavy neutral leptons (sterile neutrinos) decaying to various combinations of neutrinos, charged leptons, and hadrons. In the following sections, we will introduce these scenarios one by one and discuss the anticipated sensitivity of DUNE-PRISM, both on-axis and off-axis. Specifically, section II will be focused on dark photon-mediated DM, section III will deal with leptophobic DM, and section IV will be about heavy neutral leptons. We will discuss our findings and conclude in section V. Let us comment that similar searches for DUNE have been considered previously in refs. [9,17,27,28] for the case of on-axis detectors, and for additional data taking away from the beam axis in ref. [14].
We will go beyond these studies in two important ways: • We will investigate the usefulness of taking DUNE off-axis data for additional scenarios, for which this has never been done before. In particular, we will consider scattering of light, leptophobic DM and decays of heavy neutral leptons.
• We will also reconsider DM scattering on electrons, previously studied in ref. [14]. As a cross-check, we will reproduce the total rates analysis carried out in this reference, but we will also show that an analysis including the electron recoil spectrum is equally sensitive on-axis and off-axis.

II Light Dark Matter Interacting via a Dark Photon
Most direct searches for DM lose sensitivity at DM masses below a few GeV, motivating a new experimental program that focuses specifically on this mass range [21].
One of the simplest and most generic models for DM in the MeV-GeV mass range augments the SM by a scalar DM particle φ (or a Majorana fermion) and a new U (1) gauge boson, A . The relevant terms in the Lagrangian read with and where ] is the DM current, g is the U (1) gauge coupling, F µν and F µν are the U (1) and U (1) field strength tensors, respectively, and parameterizes the small kinetic mixing between the dark and visible photons. thus ultimately controls the interaction strength between the dark photon and SM particles. For DM lighter than half the dark photon mass (m φ < m A /2), the thermal relic abundance of φ is determined by its annihilation cross section to SM fermions, where we have defined the effective coupling strength As usual, v rel denotes the relative velocity of the two annihilating DM particles, and α D ≡ g 2 /(4π).
In the following, we will present our results in the m φ -Y plane since this choice makes it easiest to highlight those regions of parameter space where the correct DM thermal abundance is obtained [21,29]. One important feature of the annihilation cross section in eq. (4) is its v 2 rel suppression, through which strong constraints based on precise measurements of the temperature anisotropies of the cosmic microwave background radiation [30][31][32][33] are avoided. For models with unsuppressed annihilation, these constraints would rule out thermal freeze-out production of DM candidates with a mass below ∼ 10 GeV. Thanks to the velocity suppression, scalars or Majorana fermions can account for the totality of the DM abundance via thermal freeze-out even at masses below 1 GeV.
The kinetic mixing term in eq. (2) implies that any process that can create a photon can also create a dark photon, provided this is kinematically allowed. In a meson production target like the ones employed in neutrino beam experiments, dark photons can be copiously produced in meson decays such as π 0 → γA or η → γA , with a smaller contribution from bremsstrahlung. The dark photon couples to the dark current with coupling strength g , and to the SM electromagnetic current with coupling strength e. We will consider in particular the case where the dark photon mass, m A , is larger than twice the DM mass, m φ . In this case, any A produced in the target will rapidly decay, almost exclusively to φφ † . Hence, a beam of φ particles will travel alongside the neutrino beam and eventually reach the near detector, where φ particles can scatter on nuclei and electrons. It is in particular the latter channel -φ-electron scattering -that we will focus on because in this channel neutrino-induced backgrounds are smaller [2,11,15].

II.A Dark Matter Production and Detection
In a proton beam dump, dark photons with masses below ∼ 1 GeV are mainly produced in the decays of the lightest neutral mesons, π 0 and η, and in proton bremsstrahlung via the process pp → ppA . Production processes induced by leptonic secondary particles and their bremsstrahlung are usually subdominant, even if not completely negligible as reported in a detailed calculation in ref. [34]. We do not consider the latter type of processes in this work, and therefore our estimates should be considered conservative in this regard. In the mass window considered, production mechanisms that can be described in perturbative quantum chromodynamics (QCD), such as Drell-Yan production, are negligible and are not taken into account here. In the simulation of the dark photon signal, we realistically take into account correlations between the geometric acceptance of the detector and the angular spread of the DM flux. For the detector geometry, we consider for simplicity a cylindrical shape oriented along the beam axis with transverse surface given by a circle of radius 3.5 m. While this simplified detector model does not exactly match the envisioned geometry of the DUNE-PRISM detectors, the error introduced by our approximation should be negligible compared to the intrinsic uncertainties of the flux prediction.
Meson decay. The production of dark photons in meson decays occurs mostly via transitions of the form where X = π 0 , η. Reactions involving higher mass mesons are possible but usually subdominant, as shown for example in refs. [35,36]. The ρ and ω resonances, which lead to a sizeable enhancement of the production rate in a narrow mass region, are effectively taken into account within our formalism for proton bremsstrahlung, as detailed below. The transition in eq. (6) can proceed either via an on-shell or an off-shell A . We assume that the decay is dominated by the on-shell mode and make use of the formula [2,37] BR(X → γA ) BR(X → γγ) for the corresponding branching ratio. This expression has been derived in the narrow width approximation (see ref. [38] for a discussion on the full treatment of off-shell and its implication on the sensitivity). As a first step towards computing the rate of DM-electron scattering events expected in the DUNE near detectors, we need to model the spectra of π 0 s and ηs produced in DUNE's primary target. This is crucial as any systematic bias in these spectra will propagate through the simulation chain and affect the yield of signal events. This is even more relevant if one considers, as done in this work, the possibility offered by the DUNE-PRISM concept to have a movable detector. We consider two samples of mesons: • primary-only (SoftQCD), which includes only the mesons produced in the primary interaction of the proton impinging on the target, assuming that all protons eventually interact there. We employ Pythia v8.230 [39] and adopt the flag SoftQCD as the main mode for the simulation. For the simulation of relatively low-energy fixed target experiments with Pythia, the flag SoftQCD:All should be preferably used, as reported also in ref. [40].  FIG. 1. Expected number of DM-electron scattering events in the DUNE near detectors as a function of the detector position relative to the beam axis for the neutrino mode run. The background (solid gray histogram) corresponds to elastic neutrino and anti-neutrino scattering on electrons. Following ref. [14] we assume that the background from charged current quasi-elastic interactions can be made negligible by applying an energy-dependent cut on the lepton angle, which barely affects the signal and the elastic neutrino-electron scattering events. For the signal + background histograms, we show separately the two signal samples introduced in the main text: the Pythia-based sample of primary mesons, primary-only (SoftQCD) (red dot-dashed), and the sample beam-dump (blue dashed) based on the ancillary material of ref. [34]. See text for a discussion of the differences between these two samples.
• beam-dump, provided as ancillary material of ref. [34], which includes the interaction of secondary particles propagating in the DUNE target using Geant4 [41].
The beam-dump sample is more realistic and is taken as our default choice to study the DUNE sensitivity. We consider the other sample for sanity checks and for comparison with the result of ref. [14], where Pythia was employed for the simulation of mesons as well. We found discrepancies with the latter that we traced back to a different setup in the Pythia generator. 1 Further details and technical aspects on the comparison are reported in appendix A. With the meson fluxes in hand, we use MadDump [42] to generate samples of signal events of DM-electron scattering in the DUNE near detectors for all relevant parameter points. The program, once instrumented with a suitable UFO model file which encodes the coupling structure of the DM and the dark photon [43], takes care of all the steps of the simulation chain, from the decay of the parent mesons into DM particles to the detection process in the DUNE near detectors, including the computation of geometric acceptances and other effects due to the finite size of the detector.
In fig. 1 we compare the results of the two different signal predictions by plotting the rate of signal + background events in each case as a function of the detector location ∆x OA relative to the beam axis. The broad features visible in the figure are consistent with expectations: both the signal and background rates decrease as the detector is moved away from the focus direction of the beam. The signal-to-background ratio, however, generally increases with ∆x OA because the nearly massless neutrinos that are responsible for the background inherit more of the forward boost of their parent mesons than the much heavier A bosons from which the DM originates. Comparing our primary-only (SoftQCD) (red dot-dashed) sample to the beam-dump one (blue dashed), we observe that they are of similar size in the on-axis bins, but diverge with increasing off-axis angle: the event rate predicted by the beam-dump sample drops less rapidly as the detector is moved away from the beam axis. This is indeed consistent with our expectation as there is a strong correlation between the location of the detector and the energy spectrum of the DM particles, which in turn originates from the spectra of the parent mesons. In particular, when the detector is on-axis, the flux it receives is dominated by the decay products of relatively energetic mesons, which are most likely produced in the primary proton-proton interaction. Therefore, this flux is reliably modeled by the primary-only (SoftQCD) method. Off-axis, however, the DM energy spectrum is dominated by much lower-energy particles, making it much more sensitive to the decays of soft mesons produced in secondary interactions in the target. The beam-dump sample includes these mesons, explaining why the discrepancy between the blue and red curves in fig. 1 increases with ∆x OA . We conclude that the modeling of secondary interactions is crucial for the off-axis strategy.
Proton bremsstrahlung. In the dark photon mass range 500 MeV m A 1 GeV, i.e. above the η threshold, proton bremsstrahlung dominates dark photon production. Bremsstrahlung is preferentially emitted in the forward direction (collinear with the incoming proton) and in this limit can be well described by a generalization of the Fermi-Williams-Weizsäcker method [44][45][46] (or "equivalent photon method"). This method is based on the assumption that proton-nucleon scattering is dominated by exchange of vector bosons that are close to on-shell. Again we rely on MadDump for our simulations, which implements bremsstrahlung following refs. [47,48]. Let us parameterize the 4-momentum vector of the emitted A as Here, P is the momentum of the incident proton, z is the fraction of the proton momentum carried by the outgoing A , p T is the momentum perpendicular to the beam momentum, and φ is the azimuthal angle. We generate unweighted A events according to the differential production rate where σ pA (s) denotes the total interaction cross section of the incoming protons with a target nucleus of mass number A, s = 2m p E p is the square of the center-of-mass energy, and s = 2m p (E p − E A ). The ratio of cross sections σ pA (s )/σ pA (s) compares the probability that the incoming proton interacts after having emitted a photon to its total interaction probability. For the proton form factor F 1,p (m 2 A ), we use the parameterization from ref. [49] in the time-like region, so that off-shell mixing with vector mesons such as the ρ and ω is effectively included in our calculation, leading to a resonance peak in the A production rate around m A 770 MeV [50]. Finally, the photon splitting function is with H = p 2 T + (1 − z)m 2 A + z 2 m 2 p . The number of produced A events is estimated by integrating eq. (8) over a region well within the realm of the collinear approximation, i.e. a region where the kinematic conditions hold. In particular, following refs. [47,48,51], we use the integration intervals z ∈ [0.1, 0.9] and p T < 1 GeV.  , we see that, as expected, both the signal and background rates are significantly lower in the off-axis positions, especially at high energy. At low energies, where most of the events are concentrated, the signal-to-background ratio becomes better when going off-axis.

II.B Backgrounds
As the experimental signature of φ-e − scattering is a single, energetic recoil electron, the main backgrounds to this search in a neutrino beam experiment are due to neutrino-electron scattering and, if the final state hadronic system stays unidentified, charged current (CC) ν enucleon interactions. We estimate the backgrounds using GENIE v3.00.06 [52]. In particular, we simulate the relevant processes -ν-e − scattering and CC ν e interactions on nuclei -in argon. The resulting event rates are then weighted bin-by-bin by the various on-and off-axis fluxes. The simulated DUNE neutrino fluxes have been taken from ref. [53], and extrapolated to higher energies, where Monte Carlo statistics in the simulations from ref. [53] is too low. We perform this extrapolation by fitting the tail of the available fluxes linearly in log-space. The expected flux at very large neutrino energies, where the existing fluxes lack statistics, is then obtained by evaluating the fit function.
From the kinematics of elastic scattering, we know that both the signal process φe − → φe − and the background process νe − → νe − obey E e θ 2 e < 2m e , with E e (θ e ) being the final state electron's total energy (scattering angle). Neutrino-nucleon interactions, on the other hand, lead to larger scattering angles. Since the expected angular resolution of DUNE-PRISM is sufficient [54], we impose the above condition as a kinematic cut, which leaves us with a comparatively small number of events involving neutrino-nucleon interactions. We therefore neglect the latter and only consider the elastic neutrino-electron scattering contribution in our analyses.
In fig. 2, we plot the projected electron energy spectra for both the DM-electron scattering signal in the dark photon model and for the neutrino-electron scattering background. We see that both the signal and the background peak at the lowest energies, but feature a long tail towards higher energies. For the signal, the tail is most pronounced for heavier A , corresponding to heavier DM particles, as we impose m A = 3m φ here. Heavy A are produced when heavy mesons decay in the LBNF target and are therefore only kinematically accessible in very hard proton-proton collisions. The high-energy tails of both the signal and the background are strongly suppressed offaxis because the production of very energetic DM particles (for the signal) and neutrinos (for the background) occurs in production events that are strongly boosted in the forward direction. Note that the signal for m A = 0.6 GeV, m φ = 0.2 GeV (green histograms in fig. 2 drops of faster when going off-axis than the rates for lighter dark sectors. This is because at m A = 0.6 GeV production 10 -13 10 -12

II.C Statistical Analysis
To derive sensitivity limits from our predicted signal and background rates, we use standard frequentist techniques. In particular, we test the signal + background hypothesis against simulated background-only data. To do so, we use a Poissonian log-likelihood function, log L(Θ, X), defined in eq. (B1), which depends on the physical model parameters Θ = ( 4 α D , m A , m φ ) and a set of nuisance parameters X. The latter parameterize systematic normalization uncertainties and spectral "tilts" in both the signal and the background spectra. We consider systematic errors that are uncorrelated between different on-/off-axis positions (assuming 1% relative error) in addition to errors which are correlated among all positions (assuming 10% relative error). The sensitivity limits on 4 α D for fixed values of m A and m φ are determined by comparing the log-likelihood ratio, defined in eq. (B4), to the 90% quantile of a χ 2 distribution with one degree of freedom. See appendix B for further details on our statistical procedure.

II.D Sensitivity to Light Dark Matter in the Dark Photon Model
We present our main results for the dark photon model in fig. 3. In this figure, we compare on the one hand different running strategies: all data taken on-axis, all data taken off-axis, and combining data taken at different on-axis and off-axis locations as in DUNE-PRISM. On the other hand, we also compare two different analysis strategies, namely a total rates analysis (n bins = 1 in eq. (B1)) in the left panel and a spectral analysis (n bins = 80 equal-width bins up to 20 GeV) in the right panel. The former type of analysis is similar to the one discussed in ref. [14], and we confirm the main conclusion of these authors, namely that the DUNE-PRISM strategy of combining runs in seven different on-axis and off-axis locations (cyan dot-dashed) benefits the sensitivity to light DM in the dark photon model. It yields better results than both an on-axis-only run (red solid) and an off-axis-only run (purple dashed). Interestingly, though, we reach a different conclusion when including the event spectrum: as shown in the right panel of fig. 3, the sensitivity in this case is about the same for on-axis-only running and for the PRISM strategy. This can be understood by going back to fig. 2, where we see that the signal-to-background ratio at energies 2 GeV is significantly better on-axis than it is off-axis. For a total rates analysis, these high energy events do not contribute because of the steep drop of the event rate compared to the lowest energy bins. A spectral analysis, however, is able to harness also the statistical power of high-energy events and therefore suffers more from the poorer signal-to-background ratio in the off-axis position.
The overall shape of the exclusion curves in fig. 3 can be understood mostly from the φ production rate, which drops at larger masses, where fewer production modes are available. The spectral feature at m φ ∼ 230 MeV, corresponding to m A ∼ 700 MeV is related to the ρ resonance: when m A = m ρ , dark photons and ρ vector mesons exhibit maximal mixing, leading to very efficient A production and thus strong limits.
We compare the DUNE sensitivity to existing limits from BaBar [55] and NA64 [56], to a recast of NuMI off-axis data from NOνA [11] and MiniBooNE (MB) [15], and to the expected sensitivities of ICARUS-NuMI off-axis [15] and SHiP [36]. Note that we have rescaled the ICARUS-NuMI limit to an integrated luminosity of 2.5 × 10 21 protons on target (pot), corresponding to 5 years of NuMI running at the nominal beam power of 700 kW. We choose to present here some recasts that have not been officially approved by the respective experimental collaboration and for this reason have been omitted in some previous studies such as the "Physics Beyond Colliders" study at CERN [57]. We include such unofficial recasts only for neutrino experiments like NOνA [11], but not for other experiments like E137 [34,58] or BEBC [15]. Also we do not present projections for future proposed experiments other than SHiP for the same reason, for a summary see ref. [57]. We see that DUNE-PRISM can probe important new regions of parameter space, with a sensitivity to Y = 2 α D (m φ /m A ) 4 that is up to half an order of magnitude better than existing constraints for some m φ . Moreover, DUNE-PRISM can improve on the projected sensitivity for ICARUS-NuMI both in the small mass region (m A 10 MeV) and at large dark photon mass (m A 100 MeV).

III Leptophobic Dark Matter
The second scenario we are going to consider in this paper is a leptophobic DM model, where the visible and the dark sector interact via a new gauge interaction under which the leptons are neutral. As a concrete example, we will consider the force associated with a gauged baryon number symmetry U (1) B . We will call the leptophobic gauge field Z and the corresponding gauge coupling g Z . The relevant terms in the Lagrangian of the leptophobic model are then where The sum in the last term runs over all quark flavors, and z φ , z q denote the U (1) B charges of DM and of SM quarks, respectively. The leptophobic nature of the Z makes this particle rather elusive to experimental probes, as we will see in the following section. However, strong theoretical constraints arise from the fact that the U (1) B gauge symmetry is anomalous, so that extra fermions need to be added to the model to cancel anomalies [59][60][61][62][63][64]. Since these constraints are somewhat dependent on the ultraviolet completion of the model, we will not include them in our sensitivity plots.

III.A Dark Matter Production and Detection
We focus on Z candidates with masses m Z ≥ 2 GeV. In this range, the dominant production mechanism is given by prompt production, which can be described by standard methods in perturbative QCD.
The number of DM particles produced in the collision of primary protons impinging on the carbon target is given by the formula where N POT is the number of protons on target, A = 12 is the mass number of carbon (the target material used in the LBNF), and σ pA and σ pN stand for the proton-nucleus and proton-nucleon cross sections, respectively. In the above, we have assumed linear scaling with A for the cross section of DM pair production, σ pA→φφ † = Aσ pN →φφ † , and an effective total cross section per nucleon σ pA,tot = A 0.71 σ pN,tot for proton-carbon collisions [65]. We adopt the approximation σ pN,tot = 40 mb [66] and compute σ pA→φφ † numerically in MadDump, using a UFO implementation of the leptophobic model as done in ref. [42]. The computation is carried out at tree-level in QCD according to the standard parton model formula pA→φφ † is the partonic cross section for DM production in the scattering of two partons a and b, and f a/h 1 (f b/h 2 ) are the corresponding parton distribution functions of parton a (b) within hadron h 1 (h 2 ). In our simulation, we employ the leading order PDF set NNPDF2.3LO [67,68], and we fix the factorization scale to µ F = m Z . In the following, we will consider as a benchmark a scalar DM candidate φ of mass m φ = 750 MeV and U (1) B charge z φ = 3, while we vary the mass of the force carrier in the range m Z ∈ [2,7] GeV. In this case, the branching fraction BR(Z → φφ † ) is of O(1) [8] and, throughout this work, we assume it to be equal to 1.
Detection of leptophobic DM occurs via scattering of DM particles with the nuclei in the detector. Our search strategy focuses on the deep-inelastic scattering (DIS) signature, which is the most relevant one at multi-GeV energies. As for the dark photon model discussed in section II, we expect better signal/background discrimination power when going off-axis [8]. The number of signal events is computed according to the formula where z denotes the coordinate along the beam axis, ρ(z) is the number density of nucleons in the detector, S is its surface orthogonal to the z-direction, and σ φ(φ † )N →φ(φ † )N is the deep-inelastic DM-nucleon scattering cross section. The doubly differential flux d 2 N φ (E, S)/(dE dS) of DM particles per unit area dS and per unit energy dE is computed on-the-fly by MadDump.
In fig. 4 we show the geometric acceptance, which indicates what fraction of DM particles crosses the DUNE-PRISM detectors, as function of the off-axis location of the detectors for different masses m Z of the force carrier. We observe a distinctive pattern for higher masses: the acceptance steeply increases as the detectors are moved off-axis and then flattens at the largest off-axis angles, θ OA , attainable at DUNE-PRISM. This result is consistent with the ones shown in fig. 4 of ref. [8] and is a consequence of our DM candidate being a scalar. Indeed, this behavior can be understood by considering the differential production cross section for scalar DM particles, which reads dσ/dθ OA ∝ (1 − cos 2 θ OA ) sin θ OA in the center-of-mass frame. The suppression in the forward direction (θ OA = 0) is a consequence of angular momentum conservation.
Because of the increasing acceptance at large off-axis angles, combined with the expectation of lower neutrino-induced backgrounds far away from the beam axis, we consider not only the off-axis angles attainable in DUNE-PRISM, but also a hypothetical detector location at ∆x OA = 60 m (θ OA = 104.15 mrad). This distance is similar to the one of the ICARUS-NuMI detector relative to the NuMI beam. According to the investigations carried out in ref. [9], we expect that going that far off-axis leads to close-to-optimal sensitivity as it reduces the background to merely a few

III.B Backgrounds
As DM in the leptophobic model scatters only on hadrons, the main irreducible background channel is neutral current (NC) neutrino scattering. Charged current (CC) interactions might also lead to signal-like signatures when the final state charged lepton is misidentified as a charged pion. We expect the latter contribution to be smaller, as the NC and CC cross sections are of similar size, and hence we neglect it in the following.
The expected background distributions are obtained as described in section II.B, i.e. by using GENIE to simulate events and then weighting the events with the appropriate on-and off-axis neutrino fluxes. For studies at ∆x OA = 60 m, we compute the neutrino flux based on DUNE's flux Ntuples published in ref. [53], while for smaller off-axis distances we use directly the flux histograms from the same reference. Considering only data with E vis > 3 GeV allows us to neglect any contributions from resonant scattering, quasi-elastic scattering, and coherent scattering, leaving us with NC deep-inelastic neutrino scattering as the only relevant background process. The predicted signal and background spectra for the leptophobic model are shown in fig. 5. As for the dark photon model (cf. fig. 2), we find that backgrounds are dramatically suppressed when going off-axis. NC interactions with a given visible energy E vis are typically caused by neutrinos with an energy E ν E vis , implying that they are very sensitive to the high-energy tail of the neutrino spectrum. As the latter is strongly suppressed off-axis, so is the rate of NC background events to our DM search. Because of the lower backgrounds as well as the increased geometric acceptance off-axis (see fig. 4), we observe a dramatic increase in the signal-to-background ratio, especially for m Z in the multi-GeV range.

III.C Sensitivity to Leptophobic Dark Matter
To estimate the sensitivity of DUNE-PRISM to leptophobic DM, we follow the same statistical procedure as in section II.D, employing in particular the log-likelihood ratio from eqs. (B1) and (B4) which depends on the model parameters Θ = (g 6 Z , m Z , m χ ) in this case. Our results are shown in fig. 6, comparing once again a total rates (cut & count) analysis (n bins = 1, left panel) to an analysis utilizing spectral information (n bins = 57 equal-width bins between 3 GeV and 60 GeV, right panel). We find that the spectral analysis clearly outperforms the total rates fit, as can be understood from the event spectra in fig. 5. As for leptophilic DM (see section II.D), the DUNE-PRISM approach plays out its strengths especially for the total rates analysis. Unlike for leptophilic DM, however, on-axis only running is never optimal in the leptophobic model, no matter which type of analysis is used. An off-axis-only run would, however, be competitive with the DUNE-PRISM strategy if a spectral analysis is performed. This behavior can be understood from the better signal acceptance and lower backgrounds that we have discussed above in the context of figs. 4 and 5. For the same reason, a hypothetical detector at 60 m off-axis (upper black dotted line in fig. 6) would outperform DUNE-PRISM at any of the available on-axis and off-axis locations. Nevertheless, ICARUS-NuMI will do even better using off-axis neutrinos from  fig. 3), we compare a total rates analysis (left panel) to an analysis harnessing the full energy spectrum of events (right panel). We also compare to existing limits from invisible J/ψ and Υ decays at BaBar [70,71] (gray shaded regions). The black crosses indicate the exemplary model parameter points presented in fig. 5.
the NuMI beam. The excellent performance of ICARUS-NuMI can be attributed to its large mass of 480 t, compared to only 68 t for the combination of the DUNE-PRISM ND-LAr and ND-GAr detector. In addition, the NuMI neutrino flux drops off somewhat faster and becomes softer than the DUNE/LBNF flux when going off-axis. Comparing to existing limits from invisible J/ψ and Υ decays at BaBar [70,71], we find that DUNE-PRISM will compete with these existing constraints only in small regions around m Z ∼ 2 GeV and m Z ∼ 4 GeV.

IV Heavy Neutral Leptons
Heavy neutral leptons (HNLs), often also called sterile neutrinos, are SU (3) c ×SU (2) L ×U (1) Ysinglet fermions whose only (non-gravitational) coupling to the SM is via neutrino mixing. The corresponding operator is where N is the HNL field (a Weyl fermion), L are the left-handed lepton doublets,H = iσ 2 H * is the conjugate of the SM Higgs doublet field, and y is a dimensionless Yukawa coupling. Once the Higgs field acquires a vacuum expectation value, the operator in eq. (15) leads to mass mixing between the HNL and the active neutrinos. The HNL thus acquires the same couplings as the active neutrino and can be produced in any process that produces neutrinos in the SM, unless forbidden by kinematics. Consequently, given a large enough coupling y in eq. (15), meson decays in the DUNE target can copiously produce HNLs. Feynman diagrams involving the mixing operator from eq. (15) also admit HNL decays into various final states involving neutrinos, charged leptons, and/or hadrons. If the HNL is sufficiently long-lived (that is, if y is not too large), some of these decays will occur inside the near detector, leaving unique signatures.
The sensitivity of the DUNE near detectors to HNLs has been studied before in ref. [17,40], albeit for an on-axis configuration only. Current global constraints on HNLs have been compiled for instance in refs. [72][73][74][75].

IV.A Production of Heavy Neutral Leptons
In our simulations of HNL production and decay, we closely follow ref. [17], and we use a modified and expanded version of the NuShock code [76] accompanying ref. [17] 2 . More precisely, we predict the HNL flux in DUNE by starting from the simulated neutrino production events released by the DUNE collaboration [53], assuming a total luminosity of 5.5 × 10 21 pot, all taken in neutrino mode. Each event in these files, which are based on DUNE's full Monte Carlo simulation of the target station and decay volume, corresponds to the production of a single SM neutrino. To obtain HNL production events instead, we extract the kinematics of the neutrinos' parent mesons from the Monte Carlo event files and then use NuShock to re-generate their decays, enforcing the production of an HNL instead of an SM neutrino in the final state. NuShock accounts for the reduced branching ratio to the HNL final state by including an appropriate weight factor for each event.
Unfortunately, the DUNE/LBNF simulation includes only neutrino production in pion and kaon decays, but not in charm decays. This is understandable, given that neutrinos from charm decay are completely irrelevant in the SM. For HNL searches, however, charm decays are very important because at HNL masses larger than the kaon mass, they are the only kinematically allowed HNL production channel. We therefore estimate the flux of charm mesons following once again refs. [17,76]. We estimate the charm production rate based on the cross section σ cc = 2 µb for cc (open charm) production in proton-proton scattering at the center-of-mass energy √ s = 2(1 GeV) · (120 GeV) 15 GeV, see fig. 16 in ref. [77]. Taking the ratio of 12σ cc to the total proton-12 C cross section of σ pC = 331.4 mb [78] gives the rate of charm production per pot. Multiplying further by the fragmentation fraction into D s , f Ds = 0.077 [79], we obtain the rate of D s production. To generate the D s kinematics in the center-of-mass frame, we approximate the momentum dependence of the differential D s production cross section in the center-of-mass frame as [80] d 2 σ dx dp T,Ds where p T,Ds is the D s transverse momentum and x ≡ 2p z,Ds / √ s, with p x,Ds the D s momentum along the beam axis. For the numerical coefficients, we use the values n = 6.1 and b = 1.08 [80], which have been measured in the Fermilab E769 experiment using a 250 GeV proton beam impinging on a fixed target. Our predicted HNL flux as a function of HNL energy is shown in fig. 7 for different HNL masses and different off-axis angles. The left panel (M = 10 MeV) corresponds to a practically massless HNL that is dominantly produced in pion decays; the HNLs in the middle panel (M = 200 MeV) are too heavy to benefit from this production mode, and can only be produced in kaon decays, which explains the much lower flux; at M = 1 GeV (right panel), only D meson decays contribute to HNL production. In fig. 7, we have assumed |U e4 | 2 = |U µ4 | 2 = |U τ 4 | 2 = 1, implying that all meson decays involving neutrinos in the SM are replaced by the corresponding decays into HNLs.
Not surprisingly, the HNL spectrum becomes softer at larger off-axis angles because higher energy parent mesons are more strongly boosted in the forward direction. This softening of the HNL spectrum when going off-axis is fully analogous to the off-axis softening of the SM neutrino spectrum. Both on-axis and off-axis, the spectra of heavier HNLs (from kaon and charm decays) are significantly harder than the spectra of light HNLs because the smaller boost factors of heavy mesons render their decays more isotropic. For illustrative purposes, we have made the unphysical assumption |U e4 | 2 = |U µ4 | 2 = |U τ 4 | 2 = 1. In other words, we have assumed that all neutrino production processes in the DUNE/LBNF beamline are replaced by the corresponding HNL production processes. The fluxes shown here can therefore be linearly rescaled to smaller and more realistic values of |U α4 | 4 . Note the different vertical and horizontal axis scales in the three panels.

IV.B Heavy Neutral Lepton Decay
To determine the rate of HNL decays inside the DUNE near detectors, we follow the trajectories of all simulated HNLs to determine which of them cross the detectors. We consider both the segmented liquid argon time projection chamber referred to as ArgonCube or ND-LAr by the DUNE collaboration and the pressurized gaseous argon TPC ("multi-purpose detector" or ND-GAr). The ND-LAr detector is box-shaped, with a width (perpendicular to the beam axis) of 7 m, a height of 3 m, and a depth (along the beam axis) of 5 m. The ND-GAr detector is cylindrical, with the cylinder axis oriented horizontally and perpendicular to the beam axis. The detector's width is 5 m, and the cylinder radius is 2.6 m [26]. We compute the 3D coordinates at which each trajectory enters and exits the detectors by using some elementary geometry, courtesy of the trimesh Python package [81]. We then weight each HNL with the probability for decaying inside either of the two detectors.
HNL lifetimes and decay branching ratios into different final states as functions of mass are once again computed using NuShock. In setting limits, we will consider in particular the following final states: 1. νe + e − because it has the largest branching ratio among the visible decay channels for HNL masses below the muon threshold. The channel contributes appreciably also at higher masses, and, moreover, backgrounds are small in this channel.
2. νe ± µ ∓ and νµ + µ − due to the large branching ratios for heavy HNLs as well as low background levels 3. νπ 0 , e ± π ∓ and µ ± π ∓ , which dominate at intermediate HNL masses, but are also important at large masses. These channels suffer from large backgrounds that arise from NC or CC neutrino interactions with emission of an extra pion.
For each decay channel, we have used NuShock to simulate a large sample of HNL decays at rest. For each HNL in the DUNE/LBNF beam that decays inside the detector, we randomly pick one of the decay events and boost it from the HNL rest frame to the laboratory frame to obtain the 4-momenta of the observable decay products.  The HNL decay branching ratios are plotted in fig. 8 as a function of HNL mass (see also ref. [17]). We see that, at low HNL mass m π , the dominant decay mode is invisible, N → 3ν, followed by N → νe + e − . Above the pion threshold, two-body final states involving charged or neutral pions begin to dominate. Nevertheless, fully leptonic three-body final states (νe + e − , νµ + µ − , νe ± µ ∓ ) remain relevant, and N → νe ± µ ∓ even becomes the strongest visible decay mode at m N 1.2 GeV. Decay modes involving kaons, ρ mesons, and other hadrons become important as the corresponding kinematic thresholds are crossed.

IV.C Backgrounds and Analysis Cuts
To estimate the backgrounds to an HNL search in the DUNE near detectors, we have used GENIE v3.00.06 [52] to simulate CC and NC neutrino interactions on argon. We have simulated 3 × 10 6 events per flavor, one million each within each of the neutrino energy intervals [0, 20] GeV, [20,40] GeV, and [40,60] GeV. Within each interval, neutrino energies are distributed according to the on-axis DUNE flux. To obtain off-axis event samples, events are appropriately re-weighted with the ratio of the off-axis and on-axis fluxes. The rationale for dividing the simulation into three different energy ranges is to obtain more Monte Carlo statistics at high energies, where much of the sensitivity to HNLs is coming from.
We process the simulated events through the background simulation implemented in NuShock [17,76]. This involves a very simple detector simulation that applies Gaussian smearing to the momenta of the final state particles and uses a kinetic energy threshold to determine which of them are reconstructed. It rejects any events containing hadrons other than pions or atomic nuclei, thus exploiting the fact that the HNL decay modes that we consider do not contain heavy hadrons, while many potential background processes are from deep-inelastic neutrino scattering events that typically involve a lot of hadronic activity. Note that π 0 s, which are by default not decayed in GENIE, are instead decayed in NuShock and the analysis is carried out on the two final state photons. Charged pions, on the other hand, are retained undecayed due to their much longer lifetime.
The next step consists of simple particle misidentification rules. In particular: • Pions are misreconstructed as muons if their randomly chosen track length is sufficiently long (> 2 m).
• e + e − pairs are misreconstructed as (converted) photons if their angular separation is below a threshold (3 • ).
• Photons that convert after less than 2 cm are reconstructed as electrons or positrons.
The list of final state particles after threshold cuts and misidentification is then compared to the list of final state particles expected for the signal channel under consideration. Only if the number and type of each particle match between the signal and background, the event is retained (exclusive analysis). We do, however, conservatively assume no charge identification capabilities even though the ND-GAr detector will have a magnetic field and should therefore be able to efficiently distinguish positively and negatively charged final state particles.
To further suppress backgrounds, we exploit the fact that HNL decay products are typically strongly boosted in the forward direction, while background events have a more isotropic topology. We implement a cut on the angle θ between the mean direction of the two visible would-be HNL decay products and the beam axis. For given HNL mass M , we set the threshold at θ < M/(E 1 + E 2 ), where E 1 and E 2 are the energies of the two would-be HNL decay products. This threshold value is an estimate for the forward boost of the decay products in a real HNL decay.
We find the dominant background contributions for the different HNL decay channels to be as follows: 1. For N → νe + e − , by far the most important background is due to misidentified photons.
2. N → νµ + µ − is most easily mimicked by CC ν µ interactions with one real muon and one misidentified charged pion.
3. Similarly, N → νe ± µ ∓ suffers from a background due to CC ν µ interactions with a real muon and a photon that is misidentified as an electron.
4. Backgrounds to N → e ± π ∓ arise mostly from CC ν e interactions with pion production. A smaller contribution comes from NC neutrino interactions with a real pion and a photon misidentified as an electron.

5.
Similarly, N → µ ± π ∓ is affected by a large background due to CC ν µ interactions with pion production.
6. Finally, in N → νπ 0 , the source of the large background is NC neutrino interactions with pion production.
We have also considered neutrino trident production (which is not simulated by GENIE) as a possible source of background. However, based on the calculation in ref. [82], we conclude that this background will be negligible. In fig. 9, we compare the signal and background predictions for several HNL masses and decay channels, both on-axis (blue) and off-axis (purple). We have chosen the same representative HNL masses as in fig. 7, and the decay channel shown for each of them is the most sensitive one at this particular mass (in the absence of backgrounds). The values of the mixing matrix elements |U α4 | 2 (assumed to be the same for all flavors α = e, µ, τ ) are chosen the 95% CL limits at the given masses. All channels shown here involve two visible final state particles, and we show N → e ± π ∓ in the middle panel, N → ν + (π 0 → γγ) in the right panel), we show the spectra of the two visible final state particles separately. If the particles are identical (up to their charge, which we conservatively assume not to be measured), the upper histogram corresponds to the more energetic one, while the lower histogram is for the softer particle. Note that the sensitivity estimates discussed below are based on the full two-dimensional event distributions. We observe that the signal-to-background ratio improves significantly when going off-axis.
the spectra for both of them separately. (In the case of N → νπ 0 , the two final state particles are the two photons from π 0 decay.) If the two final state particles are identical (or identical up to their charge, which we assume is not measured), the upper plot shows the spectrum of the harder of the two, while the lower plot is for the softer one. We observe that, for all channels, background levels at high energies are lowered by several orders of magnitude when going off-axis. The signal spectra, on the other hand, drop somewhat more slowly. This is true especially at low energies, where the HNLs' forward boost is relatively small and so their angular distribution is more isotropic than the one of SM neutrinos, which are responsible for the background.

IV.D Sensitivity to Heavy Neutral Leptons
We are now ready to estimate the sensitivity of the DUNE near detectors to HNL decays. We use again the log-likelihood function from eq. (B1), but here we do so separately for ND-GAr and ND-LAr. This way, our analysis benefits from the much better signal-to-background ratio in ND-GAr which comes from the fact that the HNL signal scales with the detector volume, while the background scales with its mass. The two independent χ 2 values are only added up in the very end. As the HNL decay final states that we consider here contain two visible particles, we bin the events in two dimensions, corresponding to the energies of the two particles. This turns out to be very important for optimizing the sensitivity. One problem with working with two-dimensional histograms is that it is difficult to generate sufficient background Monte Carlo statistics in all bins. This can be problematic because bins that contain zero background events simply due to limited simulation statistics can lead to sensitivity estimates that are too optimistic. To mitigate this problem, we choose larger bins at higher energies: our bin width is 1 GeV for particle energies below 10 GeV, 2 GeV for particle energies up to 20 GeV, and 4 GeV for particle energies up to 40 GeV. To be on the safe side, we also identify bins in which the background prediction is exactly zero. We then consider averages of the background rate over increasingly larger neighborhoods of each such bin until we obtain a non-zero rate. If that average rate is above 0.1 events, we exclude the problematic bin from our analysis.
As systematic uncertainties, we include a 10% normalization error, which we assume to be uncorrelated between different HNL decay channels and between the signal and background. These uncertainties are described in terms of nuisance parameters with Gaussian priors. To obtain the final sensitivity, our χ 2 function is minimized over all nuisance parameters.
Our results are shown in figs. 10 and 11. The four panels in fig. 10 correspond to different assumption on the HNL couplings (coupling to ν e only, to ν µ only, to ν τ only, and to all three active neutrino flavors universally). We observe rich structure in these plots: first, we see that the projected exclusion regions have the wedge shape that is typical for long-lived particle searches. At too small mixing, HNL production is suppressed beyond the detectable level and HNLs are so longlived that they mostly decay far beyond the detector. At too large mixing, HNLs are abundantly produced, but most of them decay before reaching the detector. Away from kinematic thresholds, the upper limits on the mixing matrix elements |U e4 | 2 , |U µ4 | 2 , |U τ 4 | 2 (bottom edges of the wedgeshaped regions in fig. 10) scale roughly as |U α4 | 2 limit ∝ M −3 . This can be understood as follows: the HNL production rate scales as |U α4 | 2 M 2 over wide mass ranges, with the proportionality to M 2 reflecting the chiral suppression that otherwise affects many leptonic meson decays. The HNL decay rate in the HNL rest frame scales as |U α4 | 2 M 5 , as can be seen from dimensional analysis. In the laboratory frame, this scaling changes to |U α4 | 2 M 6 due to relativistic time dilation. Simultaneously, the opening angle of the HNL beam coming from its Lorentz boost grows with M , implying that the fraction of HNLs crossing the detector drops as M −2 . Overall, these arguments show that the experimental count rate scales as |U α4 | 4 M 6 , with deviations being observed close to kinematic thresholds. Moreover, the geometric scaling factor is not always exactly M −2 , depending on how exactly the flux in a given channel drops with angle.
Let us now discuss the kinematic thresholds visible in fig. 10. For instance, for N -ν e coupling only (top left panel of fig. 10), we see that, at M ∼ m π ∼ 140 MeV, HNL production in pion decays becomes kinematically forbidden, leading to a dent in the sensitivity in the N → νe + e − channel. Simultaneously, however, new HNL decay modes (N → e ± π ∓ and N → νπ 0 ) become allowed, compensating to some extent for this loss of sensitivity due to their large branching ratios. The next kinematic threshold occurs at M ∼ m K ∼ 490 MeV, when also HNL production in kaon decays becomes forbidden. Beyond this threshold, only HNLs from charm decay contribute, but since charm production in DUNE/LBNF is relative inefficient due to the low primary proton energy of 120 GeV, this leads to a significant drop in sensitivity. Nevertheless, the DUNE near detectors will be able to probe a large chunk of parameter space up to M ∼ m Ds ∼ 2 GeV, where also HNL production in charm decays becomes forbidden.
For HNL couplings to ν µ (top right panel of fig. 10), all thresholds are shifted by about m µ ∼ 100 GeV because each HNL needs to be produced together with a muon. For HNL couplings exclusively to ν τ (bottom left panel of fig. 10), production thresholds do not play as important a role because charm decays are the only production mode at all HNL masses due to the requirement of producing a τ alongside the HNL. The small kink visible at M ∼ 190 MeV can be understood from the D s -τ mass difference. Beyond the kink, HNL production in D s → N + τ decays is kinematically forbidden, leaving only the off-shell decays D s → ν τ + (τ * → N eν e , N µν µ , N + hadrons) as viable HNL production modes.
Comparing on-axis and off-axis sensitivities, we notice that going off-axis does not lead to significant benefits, in spite of the much better signal-to-background ratio. The reason is that, with the cuts discussed above, background suppression is fairly effective even on-axis. It is important to note, though, that going off-axis also does not significantly harm the sensitivity for any channel, in spite of the lower off-axis fluxes. This means that the search for HNLs can be carried out truly parasitically to DUNE's main oscillation program. No matter where the various near detectors are placed at any given moment, they will contribute in a useful way to the HNL sensitivity. Notably, the SAND (System for On-Axis Neutrino Detection) beam monitor, which will always remain on-axis and has not been included in our estimates, can be harnessed for the HNL search  as well.
In fig. 11, we compare DUNE's sensitivity to an array of current and future HNL constraints, combining all six HNL decay channels analyzed above. We do so for the same three scenarios as in fig. 10: on-axis running only, the realistic DUNE-PRISM running strategy with equal amounts of data collected at seven different locations relative to the beam axis (0 m, 6 m, 12 m, 18 m, 24 m, 30 m, 36 m [53]), and a hypothetical background-free DUNE-PRISM search. The comparison reveals that DUNE will be able to somewhat improve on existing limits for HNL couplings to ν e and ν µ , and will go far beyond what is currently possible for HNL couplings to ν τ . (While this conclusion hinges to some extent on our treatment of charm production in LBNF, we are confident that, given the very conservative choice we have made for the charm production cross section, the sensitivity of the real experiment will likely be better than our estimate.) Note that the strongest existing limits on HNL mixing with ν τ shown in fig. 11 are based on CHARM data [83], which has recently been reanalyzed [84]. We also remind the reader that the limits from PS191 [85] shown here have recently been called into question [86] and may be weaker than indicated in fig. 11. This would make future DUNE limits on HNL mixing with ν µ even more important.
Comparing to planned experiments that may happen on a timescale similar to DUNE, we note that the ones that are most competitive with DUNE at low HNL masses (below the kaon threshold) are the Fermilab short-baseline experiments ("SBN"), that is MicroBooNE, SBND, and ICARUS. Other planned experiments (FASER-2, CODEX-b, NA62++, MATHUSLA, SHiP, FCC-ee) probe higher HNL masses than DUNE because they operate at higher beam energies. They are, however, not able to reach the luminosities achievable in DUNE.
We should also keep in mind that the limits shown here are conservative in several ways, in particular they are based on only 5 years of data taking, and we have assumed no charge identification capabilities.

V Summary and Conclusions
In this work, we have explored the benefits of the DUNE-PRISM detectors for exploring physics beyond the SM, namely two models of light DM and a scenario featuring heavy neutral leptons (or sterile neutrinos). We have focused in particular on the capability of these detectors to move in and out of the beam axis and have found that exploiting this capability is highly beneficial in some scenarios, and never disadvantageous.
For MeV-GeV-scale scalar DM coupled via a dark photon with kinetic mixing, we have found that a detector placed off-axis will see a significant reduction in backgrounds from SM neutrino scattering. While also the signal due to DM-electron scattering is suppressed, the signal-tobackground ratio is improved significantly. This implies that a realistic DUNE-PRISM running strategy that combines on-axis and off-axis measurements will be ideal for constraining such a scenario. This is especially true if a simple cut & count analysis is performed; for a full spectral analysis, an on-axis-only measurement would perform equally well as DUNE-PRISM. Overall, a 5-year run of DUNE-PRISM (5.5 × 10 21 pot) will be able to improve existing limits on dark photon-mediated light DM by up to a factor of a few, as shown in fig. 3.
If light scalar DM couples through a leptophobic Z gauge boson rather than a dark photon, we are in the interesting situation that, for some parameter ranges (namely those with relatively heavy Z ), the DM flux increases away from the beam axis, while the background from SM neutral current neutrino interactions falls off. Consequently, the DUNE-PRISM strategy is always advantageous compared to an on-axis-only run, see fig. 6. After 5 years, DUNE-PRISM limits will be competitive with existing ones, and better in some parameter regions.
Turning finally to heavy neutral leptons ( fig. 11), we find once again that taking data both on-axis and off-axis as in DUNE-PRISM never hurts the sensitivity. For some decay channels, especially those with large backgrounds like N → νπ 0 , off-axis running is highly beneficial. For channels where backgrounds are lower or can be well distinguished from the signal by using spectral information, the sensitivity is similar off-axis and on-axis. Compared to existing limits, DUNE-PRISM after 5 years will be competitive and in some parameter regions better than existing limits for heavy neutral lepton mixing with ν e and ν µ . Significant new territory will be covered for couplings to ν τ thanks to charm production in the DUNE/LBNF target.
We conclude that the DUNE-PRISM detectors are a versatile new tool for probing numerous extensions of the SM. Their unique capability of moving off-axis will only improve the sensitivity to such scenarios, implying that a rich program of new physics searches can be carried out concurrently with DUNE's neutrino oscillation program without requiring an adaptation of the running strategy.  11. The DUNE near detectors compared to other experiments sensitive to HNLs. The top, middle, and bottom panels show the sensitivity to the squared mixing matrix elements |U e4 | 2 , |U µ4 | 2 , and |U τ 4 | 2 , respectively, assuming that only one of them is non-zero at a time. In this plot, the decay channels νe + e − , νµ + µ − , νe ± µ ∓ , νπ 0 , π ± e ∓ , π ± µ ∓ are combined. As in fig. 10, solid blue curves show the sensitivity of DUNE-PRISM (5.5 × 10 21 pot in neutrino mode, equally split between the on-axis position and six different off-axis locations). Dotted blue contours show results for on-axis running only, and thin dashed blue curves represent a hypothetical background-free analysis. The contours shown in the background correspond to existing limits (filled) and sensitivities of planned experiments (unfilled dashed and dotted). These limits are taken from the compilation in ref. [  FIG. 12. Same as fig. 1, but including the result of ref [14] in solid green, and our own primary-only (HardQCD) sample in dotted orange.
Machado for helping us in the comparison of our results to the results of ref. [14].

A Simulation of Dark Matter Production in the Dark Photon Model
In the context of the dark photon model introduced in section II, we have compared our signal and background predictions with those of ref. [14]. We have found that our estimate of the background due to neutrino-electron scattering agrees well with the results of ref. [14] in both the neutrino-and anti-neutrino mode. In the following, we focus on the estimate of the number of signal events.
As discussed in sec. II.A a crucial aspect in the simulation of the signal events is the modeling of the light meson spectra from the DUNE/LBNF target. In ref. [14], only the production of mesons in the primary proton interaction was considered, and Pythia was used as the main simulation tool. In fig. 12 we compare our results with those of fig. 3 of ref. [14] using the same benchmark point m A = 90 MeV = 3m φ , 4 α D = 10 −15 . Contrary to our expectation, we observe that our primary-only (SoftQCD) curve fails to accurately reproduce the results of that paper (solid green line) at any value of ∆x OA . As anticipated in section II.A, the source of the discrepancy is that the authors of ref. [14] have used a different configuration for their Pythia simulation. More precisely, they activated the HardQCD:All flag, and if we do the same (dotted orange curve in fig. 12), we are indeed able to reproduce their results.
We remark that using the SoftQCD:All flag is preferable for modeling the proton-proton primary interactions for the case of relatively low energy fixed target experiments such as DUNE/LBNF. (This was shown in ref. [40], which appeared after ref. [14].) In fact, primary-only (HardQCD) events tend to be much softer and have a larger angular spread than what is ex- pected from primary proton-proton interactions, as also show in fig. 13. Hence, the primary-only (HardQCD) sample underestimates the number of signal events on-axis, while being coincidentally similar to the beam-dump sample when going off-axis.

B Statistical Analysis
In the following, we describe the methodology used to derive sensitivity limits from our predicted signal and background rates in the analyses of sections II to IV in more detail. As stated in these sections we use standard frequentist techniques to test the signal + background hypothesis against simulated background-only data. We employ a Poissonian log-likelihood function, i.e. The signal and background rates in the i-th energy bin at the j-th off-axis position including systematic biases, S ij (Θ, X) and B ij (X), are defined in terms of the rates without biases, S ij (Θ, 0) and B ij (0), according to S ij (Θ, X) = 1 + X norm S,correl 1 + X norm We collectively denote the vector of physical model parameters Θ and the vector of nuisance parameters X. The nuisance parameters X norm S,corr and X norm S,j describe systematic normalization uncertainties in the signal, while the parameters X tilt S,correl and X tilt S,j parameterize spectral "tilts": their effect is to pivot the spectrum about its midpoint. The meaning of the corresponding parameters affecting the background (X norm B,corr , X norm B,j , X tilt B,correl and X tilt B,j ) is analogous. (Note that tilt errors are not considered in the HNL analysis of section IV.) All systematic errors are treated as Gaussian and are constrained by the pull terms in the second row of eq. (B1). As systematic normalization and tilt uncertainties we assume σ pos = 1% (uncorrelated between different on-/off-axis positions) and σ correl = 10% (correlated among all positions). The sums in the first line of eq. (B1) run over n bins = 80 (n bins = 57) energy bins, equally spaced in the interval [0, 20] GeV ( [3,60] GeV) in case of the dark photon model (leptophobic model), and over n pos = 7 on-and off-axis positions (0 m, 6 m/10. 45  . For comparison we also present results for an on-axis-only run, or for data taken at a fixed off-axis location. In this case, of course, we set n pos = 1. For HNL searches, we use the two-dimensional binning described in section IV.D. The statistical significance at which a given parameter point Θ is excluded can be estimated from the log-likelihood ratio [89] Z(Θ) ≡ −2 log L(Θ,X) where (Θ,X) is the combination of model parameters and nuisance parameters that maximizes the likelihood (minimizes χ 2 ), andX are the nuisance parameters that maximize the likelihood for fixed Θ. Z(Θ) follows a χ 2 distribution, with the number of degrees of freedom equal to the dimension of Θ. Consequently, the nσ exclusion region can be estimated by comparing Z(Θ) to the nσ quantile of that χ 2 distribution. However, experimental collaboration often present their results as constraints on the "signal strength" µ, which is defined as an overall, energy-independent rescaling factor of the event rate relative to some reference point. In the dark photon model, the signal strength can be taken as the product of coupling constants µ ≡ 4 α D , while in the leptophobic model it is µ ≡ g 6 Z . In this way of analyzing the data, Z(Θ) is replaced by a one-parameter function whose values are then compared to the χ 2 distribution with one degree of freedom. In eq. (B4), (μ,X) is once again the parameter point in (µ, X) at which the likelihood is maximal. In a signal strength analysis, however, only µ is varied and all other model parameters are kept fixed. The statistical question answered by a signal strength analysis is "What is the constraint on the signal strength, assuming we already know all other model parameters." A multi-parameter fit, on the other hand, asks "What are the preferred parameter regions, assuming the chosen model is the correct one, but we do not have any prior information on any of its parameters." To make our results comparable to those shown in the literature, we use signal strength analyses throughout.