Searching for Physics Beyond the Standard Model in an Off-Axis DUNE Near Detector
Moritz Breitbach, Luca Buonocore, Claudia Frugiuele, Joachim Kopp, Lukas Mittnacht
MMITP-21-004, ZU-TH-5/21, CERN-TH-2021-018
Searching for Physics Beyond the Standard Modelin an Off-Axis DUNE Near Detector
Moritz Breitbach,
1, a
Luca Buonocore,
2, b
Claudia Frugiuele,
3, c
Joachim Kopp,
4, 1, d and Lukas Mittnacht
1, e PRISMA Cluster of Excellence and Mainz Institute for Theoretical Physics,Johannes Gutenberg-Universität Mainz, Germany Physik Institut, Universität Zürich, Switzerland INFN, Sezione di Milano, Via Celoria 16, I-20133 Milano, Italy. Theoretical Physics Department, CERN, Geneva, Switzerland (Dated: February 9, 2021)Next generation neutrino oscillation experiments like DUNE and T2HK are multi-purposeobservatories, with a rich physics program beyond oscillation measurements. A special roleis played by their near detector facilities, which are particularly well-suited to search forweakly coupled dark sector particles produced in the primary target. In this paper, wedemonstrate this by estimating the sensitivity of the DUNE near detectors to the scatteringof sub-GeV DM particles and to the decay of sub-GeV sterile neutrinos (“heavy neutralleptons”). We discuss in particular the importance of the DUNE-PRISM design, whichallows some of the near detectors to be moved away from the beam axis. At such off-axislocations, the signal-to-background ratio improves for many new physics searches. We findthat this leads to a dramatic boost in the sensitivity to boosted DM particles interactingmainly with hadrons, while for boosted DM interacting with leptons, data taken on-axisleads to marginally stronger exclusion limits. Searches for heavy neutral leptons performequally well in both configurations.
Contents
I Introduction 2II Light Dark Matter Interacting via a Dark Photon 3
II.A Dark Matter Production and Detection . . . . . . . . . . . . . . . . . . . . . . . . . 4II.B Backgrounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6II.C Statistical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7II.D Sensitivity to Light Dark Matter in the Dark Photon Model . . . . . . . . . . . . . 8
III Leptophobic Dark Matter 9
III.A Dark Matter Production and Detection . . . . . . . . . . . . . . . . . . . . . . . . . 9III.B Backgrounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11III.C Sensitivity to Leptophobic Dark Matter . . . . . . . . . . . . . . . . . . . . . . . . 12
IV Heavy Neutral Leptons 12
IV.A Production of Heavy Neutral Leptons . . . . . . . . . . . . . . . . . . . . . . . . . 13IV.B Heavy Neutral Lepton Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14IV.C Backgrounds and Analysis Cuts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15IV.D Sensitivity to Heavy Neutral Leptons . . . . . . . . . . . . . . . . . . . . . . . . . . 17
V Summary and Conclusions 20A Simulation of Dark Matter Production in the Dark Photon Model 22B Statistical Analysis 23References 25 a [email protected] b [email protected] c [email protected] d [email protected] e [email protected] a r X i v : . [ h e p - ph ] F e b I Introduction
The near detectors of long-baseline neutrino experiments, once considered an afterthought toreduce systematic uncertainties in oscillation measurements, are nowadays independent experi-ments in their own right. Besides delivering a wealth of data on neutrino interaction physics, ithas been realized that they could also serve to probe the existence of physics beyond the StandardModel (SM) [1] such as heavy neutral leptons or other new particles, possibly connected to thedark matter (DM) puzzle [2–16]. Indeed, signatures of new dark particles at these experimentscan 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 atneutrino fixed target facilities by looking for scattering of the decay products off nucleons and/orelectrons in the near detector. A light DM program at neutrino facilities could complement thenext generation light DM direct detection program [17], in particular, the upcoming experimentSENSEI [18]. So far, only the Fermilab-based neutrino experiment MiniBooNE has performeddedicated searches for light DM [19, 20], but the untapped potential is big, with many past andongoing experiments having the capability to supersede MiniBooNE’s sensitivity [7, 8, 10, 14].Importantly, these searches can typically be done fully parasitically to the main neutrino program[8, 10].The near detector physics program will be taken to the next level by the DUNE-PRISM detec-tors, to be installed
574 m downstream from the target [21] at the long-baseline neutrino facility(LBNF) at Fermilab, the neutrino source for the DUNE experiment. These detectors – a liquidargon time projection chamber (TPC) and a magnetized gaseous argon TPC – will be mountedon a movable platform, allowing them to be displaced up to . (
53 mrad ) away from thebeam axis. This capability mainly serves the detectors’ primary purpose, namely constrainingthe unoscillated neutrino flux and measuring the neutrino cross sections. In particular, the neu-trino spectrum changes as a function of the off-axis angle (for kinematic reasons), while the crosssections 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 theneutrino cross section.In this paper, we will discuss the impact of the DUNE-PRISM concept on searches for physicsbeyond the SM. More specifically, we will consider the production of light ( (cid:46)
GeV ) and veryweakly interacting new particles in the target, followed by their interaction or decay inside theDUNE-PRISM detectors.Among the numerous extensions of the SM that can be probed in DUNE-PRISM and otheraccelerator neutrino experiments, we will consider in particular: (1) light ( (cid:46)
GeV ) DM parti-cles χ produced with a large Lorentz boost and detectable via dark photon-mediated χ –electronscattering; (2) light DM particles detectable via χ –nucleus scattering (leptophobic DM); and (3)heavy neutral leptons (sterile neutrinos) decaying to various combinations of neutrinos, chargedleptons, and hadrons. In the following sections, we will introduce these scenarios one by one anddiscuss the anticipated sensitivity of DUNE-PRISM, both on-axis and off-axis. Specifically, sec-tion 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 insection V. Let us comment that similar searches for DUNE have been considered previously inrefs. [8, 16] for the case of on-axis detectors, and for additional data taking away from the beamaxis in ref. [13].We will go beyond these studies in two important ways:• We will investigate the usefulness of taking DUNE off-axis data for additional scenarios, forwhich 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. [13]. As across-check, we will reproduce the total rates analysis carried out in this reference, but wewill also show that an analysis including the electron recoil spectrum is equally sensitiveon-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 newexperimental program that focuses specifically on this mass range [17].One of the simplest and most generic models for DM in the MeV–GeV mass range augmentsthe SM by a scalar DM particle φ (or a Majorana fermion) and a new U (1) (cid:48) gauge boson, A (cid:48) . Therelevant terms in the Lagrangian read L DM = L A (cid:48) + L φ , (1)with L A (cid:48) = − F (cid:48) µν F (cid:48) µν + m A (cid:48) A (cid:48) µ A (cid:48) µ − (cid:15) F (cid:48) µν F µν , (2)and L φ = ig (cid:48) A (cid:48) µ J φµ + ( ∂ µ φ † )( ∂ µ φ ) − m φ φ † φ , (3)where J φµ = (cid:2) ( ∂ µ φ † ) φ − φ † ( ∂ µ φ (cid:1) ] is the DM current, g (cid:48) is the U (1) (cid:48) gauge coupling, F µν and F (cid:48) µν are the U (1) and U (1) (cid:48) field strength tensors, respectively, and (cid:15) parameterizes the small kineticmixing between the dark and visible photons. (cid:15) thus ultimately controls the interaction strengthbetween the dark photon and SM particles. For DM lighter than half the dark photon mass( m φ < m A (cid:48) / ), the thermal relic abundance of φ is determined by its annihilation cross sectionto SM fermions, σ ( φφ † → f ¯ f ) v rel ∼ πv rel Ym φ , (4)where we have defined the effective coupling strength Y ≡ (cid:15) α D (cid:18) m φ m A (cid:48) (cid:19) . (5)As usual, v rel denotes the relative velocity of the two annihilating DM particles, and α D ≡ g (cid:48) / (4 π ) .In the following, we will present our results in the m φ – Y plane since this choice makes it easiest tohighlight those regions of parameter space where the correct DM thermal abundance is obtained[17, 22].One important feature of the annihilation cross section in eq. (4) is its v rel suppression, throughwhich strong constraints based on precise measurements of the temperature anisotropies of thecosmic microwave background radiation [23–26] are avoided. For models with unsuppressed an-nihilation, these constraints would rule out thermal freeze-out production of DM candidates witha mass below ∼
10 GeV . Thanks to the velocity suppression, scalars or Majorana fermions canaccount for the totality of the DM abundance via thermal freeze-out even at masses below .The kinetic mixing term in eq. (2) implies that any process that can create a photon can alsocreate a dark photon, provided this is kinematically allowed. In a meson production target like theones employed in neutrino beam experiments, dark photons can be copiously produced in mesondecays such as π → γA (cid:48) or η → γA (cid:48) , with a smaller contribution from bremsstrahlung. Thedark photon couples to the dark current with coupling strength g (cid:48) , and to the SM electromagneticcurrent with coupling strength (cid:15) e . We will consider in particular the case where the dark photonmass, m A (cid:48) , is larger than twice the DM mass, m φ . In this case, any A (cid:48) produced in the targetwill rapidly decay, almost exclusively to φφ † . Hence, a beam of φ particles will travel alongsidethe neutrino beam and eventually reach the near detector, where φ particles can scatter on nucleiand electrons. It is in particular the latter channel – φ –electron scattering – that we will focus onbecause in this channel neutrino-induced backgrounds are smaller [2, 27, 28]. II.A Dark Matter Production and Detection
In a proton beam dump, dark photons with masses below ∼ are mainly producedin the decays of the lightest neutral mesons, π and η , and in proton bremsstrahlung via theprocess pp → ppA (cid:48) . Production processes induced by leptonic secondary particles and theirbremsstrahlung are usually subdominant, even if not completely negligible as reported in a detailedcalculation in ref. [29]. We do not consider the latter type of processes in this work, and thereforeour 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 simulationof the dark photon signal, we realistically take into account correlations between the geometricacceptance of the detector and the angular spread of the DM flux. For the detector geometry, weconsider for simplicity a cylindrical shape oriented along the beam axis with transverse surfacegiven by a circle of radius . . While this simplified detector model does not exactly match theenvisioned geometry of the DUNE-PRISM detectors, the error introduced by our approximationshould be negligible compared to the intrinsic uncertainties of the flux prediction. Meson decay.
The production of dark photons in meson decays occurs mostly via transitionsof the form X → γA (cid:48) → γφφ † , (6)where X = π , η . Reactions involving higher mass mesons are possible but usually subdominant,as shown for example in refs. [30, 31]. The ρ and ω resonances, which lead to a sizeable enhance-ment of the production rate in a narrow mass region, are effectively taken into account within ourformalism for proton bremsstrahlung, as detailed below. The transition in eq. (6) can proceedeither via an on-shell or an off-shell A (cid:48) . We assume that the decay is dominated by the on-shellmode and make use of the formula [2, 32]BR ( X → γA (cid:48) ) BR ( X → γγ ) (cid:39) (cid:15) (cid:18) − m A (cid:48) m X (cid:19) , X = π , η (7)for the corresponding branching ratio. This expression has been derived in the narrow widthapproximation (see ref. [33] for a discussion on the full treatment of off-shell and its implicationon the sensitivity).As a first step towards computing the rate of DM–electron scattering events expected in theDUNE near detectors, we need to model the spectra of π s and η s produced in DUNE’s primarytarget. This is crucial as any systematic bias in these spectra will propagate through the simulationchain and affect the yield of signal events. This is even more relevant if one considers, as done inthis work, the possibility offered by the DUNE-PRISM concept to have a movable detector. Weconsider two samples of mesons:• primary-only (SoftQCD) , which includes only the mesons produced in the primary interac-tion of the proton impinging on the target, assuming that all protons eventually interactthere. We employ Pythia v8.230 [34] and adopt the flag
SoftQCD as the main mode forthe 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. [35].• beam-dump , provided as ancillary material of ref. [29], which includes the interaction ofsecondary particles propagating in the DUNE target using
Geant4 [36].The beam-dump sample is more realistic and is taken as our default choice to study the DUNEsensitivity. We consider the other sample for sanity checks and for comparison with the result ofref. [13], where
Pythia was employed for the simulation of mesons as well. We found discrepancies e v e n t s / po s iti on / y ea r bkgsig+bkg primary-only (SoftQCD)sig+bkg beam-dump m A ' =90 MeV = m φ , α D ε =10 -15 neutrino-mode ( s i g + bkg ) / bkg Δ x OA [m] FIG. 1. Expected number of DM–electron scattering events in the DUNE near detectors as a function ofthe detector position relative to the beam axis for the neutrino mode run. The background (solid grayhistogram) corresponds to elastic neutrino and anti-neutrino scattering on electrons. Following ref. [13]we assume that the background from charged current quasi-elastic interactions can be made negligibleby applying an energy-dependent cut on the lepton angle, which barely affects the signal and the elasticneutrino–electron scattering events. For the signal + background histograms, we show separately the twosignal 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 ofref. [29]. See text for a discussion of the differences between these two samples. with the latter that we traced back to a different setup in the
Pythia generator. Further detailsand technical aspects on the comparison are reported in appendix A.With the meson fluxes in hand, we use
MadDump [37] to generate samples of signal eventsof DM–electron scattering in the DUNE near detectors for all relevant parameter points. Theprogram, once instrumented with a suitable UFO model file which encodes the coupling structureof the DM and the dark photon [38], takes care of all the steps of the simulation chain, fromthe decay of the parent mesons into DM particles to the detection process in the DUNE neardetectors, including the computation of geometric acceptances and other effects due to the finitesize of the detector.In fig. 1 we compare the results of the two different signal predictions by plotting the rate ofsignal + background events in each case as a function of the detector location ∆ x OA relative tothe beam axis. The broad features visible in the figure are consistent with expectations: both thesignal and background rates decrease as the detector is moved away from the focus direction of thebeam. The signal-to-background ratio, however, generally increases with ∆ x OA because the nearlymassless neutrinos that are responsible for the background inherit more of the forward boost oftheir parent mesons than the much heavier A (cid:48) bosons from which the DM originates. Comparingour primary-only (SoftQCD) (red dot-dashed) sample to the beam-dump one (blue dashed), weobserve 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 movedaway from the beam axis. This is indeed consistent with our expectation as there is a strongcorrelation 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 We thank the authors of ref. [13] for helping us understand the
Pythia flags that were used in their simulation. 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 fluxis reliably modeled by the primary-only (SoftQCD) method. Off-axis, however, the DM energyspectrum is dominated by much lower-energy particles, making it much more sensitive to thedecays of soft mesons produced in secondary interactions in the target. The beam-dump sampleincludes these mesons, explaining why the discrepancy between the blue and red curves in fig. 1increases with ∆ x OA . We conclude that the modeling of secondary interactions is crucial for theoff-axis strategy. Proton bremsstrahlung.
In the dark photon mass range
500 MeV (cid:46) m A (cid:48) (cid:46) , i.e. abovethe η threshold, proton bremsstrahlung dominates dark photon production. Bremsstrahlung ispreferentially emitted in the forward direction (collinear with the incoming proton) and in thislimit can be well described by a generalization of the Fermi–Williams–Weizsäcker method [39–41](or “equivalent photon method”). This method is based on the assumption that proton–nucleonscattering 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. [42, 43]. Let usparameterize the 4-momentum vector of the emitted A (cid:48) as p A (cid:48) = ( E A (cid:48) , p T cos( φ ) , p T sin( φ ) , zP ) ,with E A (cid:48) (cid:39) zP + ( p T + m A (cid:48) ) / (2 zP ) . Here, P is the momentum of the incident proton, z is thefraction of the proton momentum carried by the outgoing A (cid:48) , p T is the momentum perpendicular tothe beam momentum, and φ is the azimuthal angle. We generate unweighted A (cid:48) events accordingto the differential production rate d N A (cid:48) d z d p T = σ pA ( s (cid:48) ) σ pA ( s ) F ,p ( m A (cid:48) ) w ba ( z, p ) , (8)where σ pA ( s ) denotes the total interaction cross section of the incoming protons with a targetnucleus of mass number A , s = 2 m p E p is the square of the center-of-mass energy, and s (cid:48) =2 m p ( E p − E A (cid:48) ) . The ratio of cross sections σ pA ( s (cid:48) ) /σ pA ( s ) compares the probability that theincoming proton interacts after having emitted a photon to its total interaction probability. Forthe proton form factor F ,p ( m A (cid:48) ) , we use the parameterization from ref. [44] in the time-likeregion, so that off-shell mixing with vector mesons such as the ρ and ω is effectively included inour calculation, leading to a resonance peak in the A (cid:48) production rate around m A (cid:48) (cid:39)
770 MeV [45].Finally, the photon splitting function is w ba ( z, p ) = (cid:15) α πH (cid:20) − z ) z − z (1 − z ) (cid:18) m p + m A (cid:48) H − z m p H (cid:19) + 2 z (1 − z )[1 + (1 − z ) ] m p m A (cid:48) H + 2 z (1 − z ) m A (cid:48) H (cid:21) , (9)with H = p T + (1 − z ) m A (cid:48) + z m p . The number of produced A (cid:48) events is estimated by integratingeq. (8) over a region well within the realm of the collinear approximation, i.e. a region where thekinematic conditions E p , E A (cid:48) , E p − E A (cid:48) (cid:29) m p , m A (cid:48) , p T (10)hold. In particular, following refs. [42, 43, 46], we use the integration intervals z ∈ [0 . , . and p T < . II.B Backgrounds
As the experimental signature of φ – e − scattering is a single, energetic recoil electron, themain backgrounds to this search in a neutrino beam experiment are due to neutrino–electron electron energy [GeV] e v e n t s / G e V /5 . × p o t /68 t e l a s t i c ν – e − neutrino-mode, on-axis electron energy [GeV] off-axis (18 m) electron energy [GeV] m A = 3 m φ off-axis (30 m) m A = 6 MeV (cid:15) α D = 10 -19 m A = 60 MeV (cid:15) α D = 10 -16 m A = 0 . (cid:15) α D = 10 -11 FIG. 2. Signal and background spectra in DUNE-PRISM for the dark photon model after 5 years of datataking ( . × pot ) in neutrino mode. Colored histograms show the DM–electron scattering signal fordifferent exemplary sets of model parameter points, imposing m A (cid:48) = 3 m φ . The gray shaded backgroundhistogram is based on a simulation of elastic neutrino–electron scattering. Comparing data taken all on-axis (left panel) to data taken at
18 m ( .
36 mrad ) off-axis (center panel) or
30 m ( .
26 mrad ) off-axis(right panel), we see that, as expected, both the signal and background rates are significantly lower in theoff-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. scattering and, if the final state hadronic system stays unidentified, charged current (CC) ν e –nucleon interactions. We estimate the backgrounds using GENIE v3.00.06 [47]. In particular, wesimulate 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. [48], and extrapolated to higherenergies, where Monte Carlo statistics in the simulations from ref. [48] is too low. We performthis extrapolation by fitting the tail of the available fluxes linearly in log-space. The expectedflux at very large neutrino energies, where the existing fluxes lack statistics, is then obtained byevaluating the fit function.From the kinematics of elastic scattering, we know that both the signal process φe − → φe − andthe background process νe − → νe − obey E e θ e < m e , with E e ( θ e ) being the final state electron’stotal energy (scattering angle). Neutrino–nucleon interactions, on the other hand, lead to largerscattering angles. Since the expected angular resolution of DUNE-PRISM is sufficient [49], weimpose the above condition as a kinematic cut, which leaves us with a comparatively small numberof events involving neutrino–nucleon interactions. We therefore neglect the latter and only considerthe elastic neutrino–electron scattering contribution in our analyses.In fig. 2, we plot the projected electron energy spectra for both the DM–electron scatteringsignal in the dark photon model and for the neutrino–electron scattering background. We see thatboth the signal and the background peak at the lowest energies, but feature a long tail towardshigher energies. For the signal, the tail is most pronounced for heavier A (cid:48) , corresponding to heavierDM particles, as we impose m A (cid:48) = 3 m φ here. Heavy A (cid:48) are produced when heavy mesons decayin the LBNF target and are therefore only kinematically accessible in very hard proton–protoncollisions. The high-energy tails of both the signal and the background are strongly suppressedoff-axis because the production of very energetic DM particles (for the signal) and neutrinos (forthe background) occurs in production events that are strongly boosted in the forward direction. II.C Statistical Analysis
To derive sensitivity limits from our predicted signal and background rates, we use standardfrequentist techniques. In particular, we test the signal + background hypothesis against simulated -3 -2 -1 m φ [GeV] -13 -12 -11 -10 -9 -8 Y = (cid:15) α D ( m χ / m A ) M B N O ν A I C A R U S - N u M I BaBar
N A S H i P α D = 0 . m A = 3 m φ total rates analysis DUNEon-axisoff-axis (30 m)PRISM relic density (complex scalar)
90% CL -3 -2 -1 m φ [GeV] -13 -12 -11 -10 -9 -8 Y = (cid:15) α D ( m χ / m A ) M B N O ν A I C A R U S - N u M I BaBar
N A S H i P α D = 0 . m A = 3 m φ spectral analysis (∆ E = 250 MeV) DUNEon-axisoff-axis (30 m)PRISM relic density (complex scalar)
90% CL
FIG. 3. Expected upper exclusion limits for the dark photon model, assuming a total DUNE-PRISMrunning time of 5 years ( . × pot) in neutrino mode. We compare results for on-axis-only running (redsolid), off-axis-only running (purple dashed), and a realistic DUNE-PRISM strategy with equal amountsof data taken at seven different locations (cyan dot-dashed). The left panel corresponds to a total rates(cut & count) analysis as in ref. [13], while the analysis in the right panel includes spectral information.Gray shaded regions indicate existing limits on the model, while dotted black lines show projections forother future experiments. Black crosses indicate the exemplary model parameter points presented in fig. 2. background-only data. To do so, we use a Poissonian log-likelihood function, log L ( Θ , X ) , definedin eq. (B1), which depends on the physical model parameters Θ = ( (cid:15) α D , m A (cid:48) , m φ ) and a setof nuisance parameters X . The latter parameterize systematic normalization uncertainties andspectral “tilts” in both the signal and the background spectra. We consider systematic errors thatare uncorrelated between different on-/off-axis positions (assuming 1% relative error) in additionto errors which are correlated among all positions (assuming 10% relative error). The sensitivitylimits on (cid:15) α D for fixed values of m A (cid:48) and m φ are determined by comparing the log-likelihoodratio, defined in eq. (B4), to the 90% quantile of a χ 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 compareon the one hand different running strategies: all data taken on-axis, all data taken off-axis, andcombining data taken at different on-axis and off-axis locations as in DUNE-PRISM. On the otherhand, we also compare two different analysis strategies, namely a total rates analysis ( n bins = 1 ineq. (B1)) in the left panel and a spectral analysis ( n bins = 80 equal-width bins up to
20 GeV ) in theright panel. The former type of analysis is similar to the one discussed in ref. [13], and we confirmthe main conclusion of these authors, namely that the DUNE-PRISM strategy of combining runsin seven different on-axis and off-axis locations (cyan dot-dashed) benefits the sensitivity to lightDM 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 conclusionwhen including the event spectrum: as shown in the right panel of fig. 3, the sensitivity in this caseis about the same for on-axis-only running and for the PRISM strategy. This can be understoodby going back to fig. 2, where we see that the signal-to-background ratio at energies (cid:38) issignificantly better on-axis than it is off-axis. For a total rates analysis, these high energy eventsdo 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 eventsand 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. Thespectral feature at m φ ∼
230 MeV , corresponding to m A (cid:48) ∼
700 MeV is related to the ρ resonance:when m A (cid:48) = m ρ , dark photons and ρ vector mesons exhibit maximal mixing, leading to veryefficient A (cid:48) production and thus strong limits.We compare the DUNE sensitivity to existing limits from BaBar [50] and NA64 [51], to a recastof NuMI off-axis data from NO ν A [27] and MiniBooNE (MB) [28], and to the expected sensitivitiesof ICARUS-NuMI off-axis [28] and SHiP [31]. Note that we have rescaled the ICARUS-NuMI limitto an integrated luminosity of . × protons on target (pot), corresponding to 5 years of NuMIrunning at the nominal beam power of
700 kW . We choose to present here some recasts that havenot been officially approved by the respective experimental collaboration and for this reason havebeen omitted in some previous studies such as the “Physics Beyond Colliders” study at CERN[52]. We include such unofficial recasts only for neutrino experiments like NO ν A [10], but not forother experiments like E137 [29, 53] or BEBC [28]. Also we do not present projections for futureproposed experiments other than SHiP for the same reason, for a summary see ref. [52]. We seethat DUNE-PRISM can probe important new regions of parameter space, with a sensitivity to Y = (cid:15) α D ( m χ /m A (cid:48) ) that is up to half an order of magnitude better than existing constraints forsome m φ . Moreover, DUNE-PRISM can improve on the projected sensitivity for ICARUS-NuMIboth in the small mass region ( m A (cid:48) (cid:46)
10 MeV ) and at large dark photon mass ( m A (cid:48) (cid:38)
100 MeV ). III Leptophobic Dark Matter
The second scenario we are going to consider in this paper is a leptophobic DM model, wherethe visible and the dark sector interact via a new gauge interaction under which the leptonsare neutral. As a concrete example, we will consider the force associated with a gauged baryonnumber symmetry U (1) B . We will call the leptophobic gauge field Z (cid:48) and the corresponding gaugecoupling g Z . The relevant terms in the Lagrangian of the leptophobic model are then L leptophobic ⊃ ig Z z φ Z (cid:48) µ J φµ + ∂ µ φ † ∂ µ φ − m φ φ † φ + g Z z q (cid:88) q ¯ qγ µ q Z (cid:48) µ , (11)where J φµ = (cid:2) ( ∂ µ φ † ) φ − φ † ( ∂ µ φ (cid:1) ] . 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 (cid:48) 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 isanomalous, so that extra fermions need to be added to the model to cancel anomalies [54–56].Since these constraints are somewhat dependent on the ultraviolet completion of the model, wewill not include them in our sensitivity plots. III.A Dark Matter Production and Detection
We focus on Z (cid:48) candidates with masses m Z (cid:48) ≥ . In this range, the dominant produc-tion mechanism is given by prompt production, which can be described by standard methods inperturbative QCD.The number of DM particles produced in the collision of primary protons impinging on thecarbon target is given by the formula N φ = 2 N POT σ pA → φφ † σ pA, tot = A . N POT σ pN → φφ † σ pN, tot (12)where N POT is the number of protons on target, A = 12 is the mass number of carbon (the targetmaterial used in the LBNF), and σ pA and σ pN stand for the proton–nucleus and proton–nucleon0 -3 -2 -1 ε [ % ] Δ x OA [m] m Z ' =2 GeVm Z ' =4 GeVm Z ' =7 GeV geometric acceptance
FIG. 4. Geometric acceptance for leptophobic scalar DM particles as a function of the detector locationrelative to the beam axis. The different colored lines correspond to different masses of the Z (cid:48) mediator.The drop of the acceptance when going towards the on-axis location is a consequence of the scalar natureof the DM candidate. cross sections, respectively. In the above, we have assumed linear scaling with A for the crosssection of DM pair production, σ pA → φφ † = Aσ pN → φφ † , and an effective total cross section pernucleon σ pA, tot = A . σ pN, tot for proton–carbon collisions [57]. We adopt the approximation σ pN, tot = 40 mb [58] and compute σ pA → φφ † numerically in MadDump , using a UFO implementa-tion of the leptophobic model as done in ref. [37]. The computation is carried out at tree-level inQCD according to the standard parton model formula σ pA → φφ † = (cid:88) a,b (cid:90) d x d x f a/h ( x ) f b/h ( x ) ˆ σ ( a,b ) pA → φφ † , (13)where ˆ σ ( a,b ) pA → φφ † is the partonic cross section for DM production in the scattering of two partons a and b , and f a/h ( f b/h ) are the corresponding parton distribution functions of parton a ( b ) withinhadron h ( h ). In our simulation, we employ the leading order PDF set NNPDF2.3LO [59, 60], andwe fix the factorization scale to µ F = m Z (cid:48) . In the following, we will consider as a benchmark ascalar DM candidate φ of mass m φ = 750 MeV and U (1) B charge z φ = 3 , while we vary the mass ofthe force carrier in the range m Z (cid:48) ∈ [2 ,
7] GeV . In this case, the branching fraction BR ( Z (cid:48) → φφ † ) is of O (1) [7] 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 thedetector. Our search strategy focuses on the deep-inelastic scattering (DIS) signature, which isthe 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 [7]. The number ofsignal events is computed according to the formula N sig = (cid:90) d z ρ ( z ) (cid:90) S d S (cid:90) d E d N φ ( E, S )d E d S σ φ ( φ † ) N → φ ( φ † ) N , (14)where z denotes the coordinate along the beam axis, ρ ( z ) is the number density of nucleons in thedetector, S is its surface orthogonal to the z -direction, and σ φ ( φ † ) N → φ ( φ † ) N is the deep-inelasticDM–nucleon scattering cross section. The doubly differential flux d N φ ( E, S ) / (d E d S ) of DMparticles per unit area d S and per unit energy d E is computed on-the-fly by MadDump .In fig. 4 we show the geometric acceptance, which indicates what fraction of DM particlescrosses the DUNE-PRISM detectors, as function of the off-axis location of the detectors for dif-ferent masses m Z (cid:48) of the force carrier. We observe a distinctive pattern for higher masses: the1 Experiment Baseline Mass Off-axis Exposure n bkg ε m Z (cid:48) =2 GeV ε m Z (cid:48) =7 GeV [m] [t] [mrad] [pot]DUNE on-axis 754 68 0 . × . × . × − . × − DUNE @
30 m
754 68 52.26 . × . × . × − . × − DUNE @
60 m
754 68 104.15 . × . × − . × − ICARUS-NuMI 789 480 ∼
100 2 . × . × − . × − TABLE I. Summary of experimental setups considered for leptophobic DM. We have used a primaryproton energy of
120 GeV and the expected exposure after 5 years for both DUNE and ICARUS-NuMI. Acut on the visible energy ( E vis > ) is applied to all signal and background events, with n bkg beingthe number of background events that pass the cut. The two rightmost columns indicate the geometricacceptance at two different Z (cid:48) masses. acceptance steeply increases as the detectors are moved off-axis and then flattens at the largestoff-axis angles, θ OA , attainable at DUNE-PRISM. This result is consistent with the ones shownin fig. 4 of ref. [7] and is a consequence of our DM candidate being a scalar. Indeed, this behaviorcan be understood by considering the differential production cross section for scalar DM particles,which reads d σ/ d θ OA ∝ (1 − cos θ OA ) sin θ OA in the center-of-mass frame. The suppression inthe 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 oflower neutrino-induced backgrounds far away from the beam axis, we consider not only the off-axisangles 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 relativeto the NuMI beam. According to the investigations carried out in ref. [8], we expect that goingthat far off-axis leads to close-to-optimal sensitivity as it reduces the background to merely a fewthousand events for an exposure of about . × pot. While exploiting off-axis distances aslarge as
60 m in DUNE would require significant civil construction, doing so might be interestingif hints for leptophobic DM interactions should be found in DUNE-PRISM. We use
MadDump to compute the simulated deep-inelastic neutrino–nucleon scattering events within the detector.For reference, in table I we collect the main parameters of different DUNE-PRISM configurationsas well as ICARUS-NuMI. In simulating ICARUS-NuMI, we have used the NuMI flux from [61].
III.B Backgrounds
As DM in the leptophobic model scatters only on hadrons, the main background channel isneutral current (NC) neutrino scattering. The expected background distributions are obtained asdescribed in section II.B, i.e. by using GENIE to simulate events and then weighting the eventswith the appropriate on- and off-axis neutrino fluxes. For studies at ∆ x OA = 60 m , we computethe neutrino flux based on DUNE’s flux Ntuples published in ref. [48], while for smaller off-axisdistances we use directly the flux histograms from the same reference. Considering only datawith E vis > allows us to neglect any contributions from resonant scattering, quasi-elasticscattering, and coherent scattering, leaving us with NC deep-inelastic neutrino scattering as theonly relevant background process.The predicted signal and background spectra for the leptophobic model are shown in fig. 5. Asfor the dark photon model (cf. fig. 2), we find that backgrounds are dramatically suppressed whengoing off-axis. NC interactions with a given visible energy E vis are typically caused by neutrinoswith an energy E ν (cid:29) E vis , implying that they are very sensitive to the high-energy tail of theneutrino spectrum. As the latter is strongly suppressed off-axis, so is the rate of NC backgroundevents to our DM search. Because of the lower backgrounds as well as the increased geometricacceptance off-axis (see fig. 4), we observe a dramatic increase in the signal-to-background ratio,especially for m Z (cid:48) in the multi-GeV range.2 visible energy [GeV] e v e n t s / G e V /5 . × p o t /68 t N C D I S neutrino-mode, on-axis visible energy [GeV] off-axis (30 m) visible energy [GeV] m χ = 750 MeV off-axis (60 m) m Z = 2 GeV g Z = 10 -5 m Z = 4 GeV g Z = 10 -4 m Z = 7 GeV g Z = 10 -3 FIG. 5. Predicted DM scattering signal in the leptophobic model and corresponding backgrounds forDUNE-PRISM with . × pot of neutrino-mode data, corresponding to 5 years of running. Coloredhistograms show the signal rate for different exemplary model parameter points. The gray shaded back-ground histogram is based on a simulation of NC deep-inelastic neutrino–nucleus scattering. The left panelis for the on-axis location, while the center and right panels correspond to data taken at
30 m ( .
26 mrad )off-axis and at a hypothetical detector location
60 m ( .
15 mrad ) off-axis, respectively. The red dashedline indicates the energy threshold E vis > imposed in our analysis. We observe that backgroundsare significantly suppressed off-axis, especially at high energies. The signal, on the other hand, is evenenhanced there for large Z (cid:48) masses, as explained in section III.A. III.C Sensitivity to Leptophobic Dark Matter
To estimate the sensitivity of DUNE-PRISM to leptophobic DM, we follow the same statisti-cal 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 Z , m Z (cid:48) , m χ ) in this case. Our resultsare shown in fig. 6, comparing once again a total rates (cut & count) analysis ( n bins = 1 , leftpanel) to an analysis utilizing spectral information ( n bins = 57 equal-width bins between and
60 GeV , right panel). We find that the spectral analysis clearly outperforms the total ratesfit, 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. Unlikefor leptophilic DM, however, on-axis only running is never optimal in the leptophobic model, nomatter which type of analysis is used. An off-axis-only run would, however, be competitive withthe DUNE-PRISM strategy if a spectral analysis is performed. This behavior can be understoodfrom the better signal acceptance and lower backgrounds that we have discussed above in thecontext of figs. 4 and 5. For the same reason, a hypothetical detector at
60 m off-axis (upperblack dotted line in fig. 6) would outperform DUNE-PRISM at any of the available on-axis andoff-axis locations. Nevertheless, ICARUS-NuMI will do even better using off-axis neutrinos fromthe NuMI beam. The excellent performance of ICARUS-NuMI can be attributed to its large massof
480 t , compared to only
68 t for the combination of the DUNE-PRISM ND-LAr and ND-GArdetector. In addition, the NuMI neutrino flux drops off somewhat faster and becomes softer thanthe DUNE/LBNF flux when going off-axis. Comparing to existing limits from invisible
J/ψ and Υ decays at BaBar [62, 63], we find that DUNE-PRISM will compete with these existing constraintsonly in small regions around m Z (cid:48) ∼ and m Z (cid:48) ∼ . IV Heavy Neutral Leptons
Heavy neutral leptons (HNLs), often also called sterile neutrinos, are SU (3) c × SU (2) L × U (1) Y -singlet fermions whose only (non-gravitational) coupling to the SM is via neutrino mixing. The3 m Z [GeV] -1 g Z Υ → χχ J/ Ψ → χχ D U N E @ m I C A R U S - N u M I m χ = 750 MeV total rates analysis DUNEon-axisoff-axis (30 m)PRISM
90% CL m Z [GeV] -1 g Z Υ → χχ J/ Ψ → χχ D U N E @ m I C A R U S - N u M I m χ = 750 MeV spectral analysis (∆ E = 1 GeV) DUNEon-axisoff-axis (30 m)PRISM
90% CL
FIG. 6. Projected upper exclusion limits for the leptophobic DM model, assuming the DM particle φ carries a U (1) B charge of z φ = 3 . Colored exclusion curves correspond to different DUNE-PRISM runningstrategies namely on-axis-only running (red solid), off-axis-only running (purple dashed), and equal runningtimes at seven different locations (cyan dot-dashed). Black dotted exclusion curves show the constraintsthat could be obtained at ICARUS-NuMI using neutrinos from the NuMI beam and at a hypotheticaldetector placed
60 m off-axis in the DUNE/LBNF beam. Like for the dark photon model (cf. fig. 3), wecompare 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 [62, 63] (grayshaded regions). The black crosses indicate the exemplary model parameter points presented in fig. 5. corresponding operator is L HNL ⊃ y ¯ L ˜ HN , (15)where N is the HNL field (a Weyl fermion), L are the left-handed lepton doublets, ˜ H = iσ H ∗ isthe conjugate of the SM Higgs doublet field, and y is a dimensionless Yukawa coupling. Once theHiggs field acquires a vacuum expectation value, the operator in eq. (15) leads to mass mixingbetween the HNL and the active neutrinos. The HNL thus acquires the same couplings as theactive neutrino and can be produced in any process that produces neutrinos in the SM, unlessforbidden by kinematics. Consequently, given a large enough coupling y in eq. (15), meson decaysin the DUNE target can copiously produce HNLs.Feynman diagrams involving the mixing operator from eq. (15) also admit HNL decays intovarious final states involving neutrinos, charged leptons, and/or hadrons. If the HNL is sufficientlylong-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. [16], albeitfor an on-axis configuration only. Current global constraints on HNLs have been compiled forinstance in refs. [64–67]. IV.A Production of Heavy Neutral Leptons
In our simulations of HNL production and decay, we closely follow ref. [16], and we use amodified and expanded version of the
NuShock code [68] accompanying ref. [16]. More precisely,we predict the HNL flux in DUNE by starting from the simulated neutrino production eventsreleased by the DUNE collaboration [48], assuming a total luminosity of . × pot , all taken inneutrino mode. Each event in these files, which are based on DUNE’s full Monte Carlo simulationof the target station and decay volume, corresponds to the production of a single SM neutrino. To4obtain HNL production events instead, we extract the kinematics of the neutrinos’ parent mesonsfrom the Monte Carlo event files and then use NuShock to re-generate their decays, enforcingthe production of an HNL instead of an SM neutrino in the final state.
NuShock accounts forthe reduced branching ratio to the HNL final state by including an appropriate weight factor foreach event.Unfortunately, the DUNE/LBNF simulation includes only neutrino production in pion andkaon decays, but not in charm decays. This is understandable, given that neutrinos from charmdecay are completely irrelevant in the SM. For HNL searches, however, charm decays are veryimportant because at HNL masses larger than the kaon mass, they are the only kinematicallyallowed HNL production channel. We therefore estimate the flux of charm mesons followingonce again refs. [16, 68]. We estimate the charm production rate based on the cross section σ c ¯ c = 2 µ b for c ¯ c (open charm) production in proton–proton scattering at the center-of-massenergy √ s = (cid:112) · (120 GeV) (cid:39)
15 GeV , see fig. 16 in ref. [69]. Taking the ratio of σ c ¯ c to the total proton– C cross section of σ pC = 331 . [70] gives the rate of charm production perpot. Multiplying further by the fragmentation fraction into D s , f D s = 0 . [71], we obtain therate of D s production. To generate the D s kinematics in the center-of-mass frame, we approximatethe momentum dependence of the differential D s production cross section in the center-of-massframe as [72] d σ d x d p T,D s ∝ (1 − | x | ) n exp( − b p T,D s ) , (16)where p T,D s is the D s transverse momentum and x ≡ p z,D s / √ s , with p x,D s the D s momentumalong the beam axis. For the numerical coefficients, we use the values n = 6 . and b = 1 . [72],which have been measured in the Fermilab E769 experiment using a
250 GeV proton beam im-pinging on a fixed target.Our predicted HNL flux as a function of HNL energy is shown in fig. 7 for different HNL massesand different off-axis angles. The left panel ( M = 10 MeV ) corresponds to a practically masslessHNL 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 contributeto HNL production. In fig. 7, we have assumed | U e | = | U µ | = | U τ | = 1 , implying that allmeson 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 higherenergy parent mesons are more strongly boosted in the forward direction. This softening of theHNL spectrum when going off-axis is fully analogous to the off-axis softening of the SM neutrinospectrum. 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 heavymesons render their decays more isotropic. IV.B Heavy Neutral Lepton Decay
To determine the rate of HNL decays inside the DUNE near detectors, we follow the trajectoriesof all simulated HNLs to determine which of them cross the detectors. We consider both thesegmented liquid argon time projection chamber referred to as ArgonCube or ND-LAr by theDUNE 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 ,a height of , and a depth (along the beam axis) of . The ND-GAr detector is cylindrical,with the cylinder axis oriented horizontally and perpendicular to the beam axis. The detector’swidth is , and the cylinder radius is . [21]. We compute the 3D coordinates at whicheach trajectory enters and exits the detectors by using some elementary geometry, courtesy of the trimesh Python package [73]. We then weight each HNL with the probability for decaying insideeither of the two detectors.5 HN L s /5 . × p o t / c m / G e V ] M = 10 MeV HN L s /5 . × p o t / c m / G e V ] M = 200 MeV HN L s /5 . × p o t / c m / G e V ] M = 1000 MeV
FIG. 7. The differential flux of HNLs at the DUNE near detector site (
574 m baseline) for different HNLmasses (left, middle, and right panels) and for different off-axis angles (different colors and line styles).For illustrative purposes, we have made the unphysical assumption | U e | = | U µ | = | U τ | = 1 . In otherwords, we have assumed that all neutrino production processes in the DUNE/LBNF beamline are replacedby the corresponding HNL production processes. The fluxes shown here can therefore be linearly rescaledto smaller and more realistic values of | U α | . Note the different vertical and horizontal axis scales in thethree panels. HNL lifetimes and decay branching ratios into different final states as functions of mass are onceagain computed using
NuShock . In setting limits, we will consider in particular the followingfinal states:1. νe + e − because it has the largest branching ratio among the visible decay channels forHNL masses below the muon threshold. The channel contributes appreciably also at highermasses, and, moreover, backgrounds are small in this channel.2. νe ± µ ∓ and νµ + µ − due to the large branching ratios for heavy HNLs as well as low back-ground levels3. νπ , e ± π ∓ and µ ± π ∓ , which dominate at intermediate HNL masses, but are also importantat large masses. These channels suffer from large backgrounds that arise from NC or CCneutrino interactions with emission of an extra pion.For each decay channel, we have used NuShock to simulate a large sample of HNL decays atrest. For each HNL in the DUNE/LBNF beam that decays inside the detector, we randomly pickone of the decay events and boost it from the HNL rest frame to the laboratory frame to obtainthe 4-momenta of the observable decay products.The HNL decay branching ratios are plotted in fig. 8 as a function of HNL mass (see alsoref. [16]). We see that, at low HNL mass (cid:46) m π , the dominant decay mode is invisible, N → ν ,followed by N → νe + e − . Above the pion threshold, two-body final states involving chargedor 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 modeat m N (cid:38) . . Decay modes involving kaons, ρ mesons, and other hadrons become importantas 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 usedGENIE v3.00.06 [47] to simulate CC and NC neutrino interactions on argon. We have simulated6 − − − − b r a n c h i n g r a t i o ν η e ± ρ ∓ µ ± K ∓ ν ω µ ± K ∗ ∓ ν e + e − µ ± ρ ∓ ν e ± µ ∓ e ± K ∗ ∓ ν µ + µ − µ ± π ∓ νη ν γ e ± K ∓ e ± π ∓ ν φ ν ν ρ ν π FIG. 8. Branching ratios of HNL decay modes as a function of mass. We have assumed flavor-universalmixing, i.e. | U e | = | U µ | = | U τ | . × events per flavor, one million each within each of the neutrino energy intervals [0 , GeV, [20 , GeV, and [40 , GeV. Within each interval, neutrino energies are distributed accordingto the on-axis DUNE flux. To obtain off-axis event samples, events are appropriately re-weightedwith the ratio of the off-axis and on-axis fluxes. The rationale for dividing the simulation intothree different energy ranges is to obtain more Monte Carlo statistics at high energies, where muchof the sensitivity to HNLs is coming from.We process the simulated events through the background simulation implemented in
NuShock [16, 68]. This involves a very simple detector simulation that applies Gaussian smearing to themomenta of the final state particles and uses a kinetic energy threshold to determine which ofthem are reconstructed. It rejects any events containing hadrons other than pions or atomicnuclei, thus exploiting the fact that the HNL decay modes that we consider do not contain heavyhadrons, while many potential background processes are from deep-inelastic neutrino scatteringevents that typically involve a lot of hadronic activity. Note that π s, which are by default notdecayed in GENIE, are instead decayed in NuShock and the analysis is carried out on the twofinal state photons. Charged pions, on the other hand, are retained undecayed due to their muchlonger lifetime.The next step consists of simple particle misidentification rules. In particular:• Pions are misreconstructed as muons if their randomly chosen track length is sufficientlylong ( > ).• e + e − pairs are misreconstructed as (converted) photons if their angular separation is belowa threshold ( ◦ ).• Photons that convert after less than are reconstructed as electrons or positrons.The list of final state particles after threshold cuts and misidentification is then compared tothe list of final state particles expected for the signal channel under consideration. Only if the7number 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 capabilitieseven though the ND-GAr detector will have a magnetic field and should therefore be able toefficiently distinguish positively and negatively charged final state particles.To further suppress backgrounds, we exploit the fact that HNL decay products are typicallystrongly 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-beHNL decay products and the beam axis. For given HNL mass M , we set the threshold at θ We are now ready to estimate the sensitivity of the DUNE near detectors to HNL decays. Weuse again the log-likelihood function from eq. (B1), but here we do so separately for ND-GArand ND-LAr. This way, our analysis benefits from the much better signal-to-background ratio in8 − e ± M = 10 MeV, | U α | = 10 − e ± M = 10 MeV, | U α | = 10 − N → νe + e − − e v e n t s /5 . × p o t / G e V e ± e v e n t s /5 . × p o t / G e V e ± − e ± M = 200 MeV, | U α | = 10 − e ± M = 200 MeV, | U α | = 10 − N → e ± π ∓ − e v e n t s /5 . × p o t / G e V π ± e v e n t s /5 . × p o t / G e V π ± − γ M = 1000 MeV, | U α | = 10 − γ M = 1000 MeV, | U α | = 10 − N → νπ − e v e n t s /5 . × p o t / G e V γ e v e n t s /5 . × p o t / G e V γ FIG. 9. Signal and background event spectra at the DUNE near detector complex for specific HNL decaymodes. Blue and purple unshaded histograms correspond to the number of signal events per GeV at theon-axis location and at 18 m off-axis, respectively. Shaded histograms represent the corresponding back-grounds. The three panels correspond to the most sensitive decay channels in the absence of backgroundsat M = 10 MeV (left panel), M = 200 MeV (middle panel), and M = 1000 MeV (right panel), respectively,as indicated in the plots. For all decay channels included here ( N → νe + e − in the left panel, N → e ± π ∓ in the middle panel, N → ν + ( π → γγ ) in the right panel), we show the spectra of the two visiblefinal state particles separately. If the particles are identical (up to their charge, which we conservativelyassume not to be measured), the upper histogram corresponds to the more energetic one, while the lowerhistogram is for the softer particle. Note that the sensitivity estimates discussed below are based on the fulltwo-dimensional event distributions. We observe that the signal-to-background ratio improves significantlywhen going off-axis. ND-GAr which comes from the fact that the HNL signal scales with the detector volume, whilethe background scales with its mass. The two independent χ values are only added up in the veryend. As the HNL decay final states that we consider here contain two visible particles, we bin theevents in two dimensions, corresponding to the energies of the two particles. This turns out to bevery important for optimizing the sensitivity. One problem with working with two-dimensionalhistograms is that it is difficult to generate sufficient background Monte Carlo statistics in allbins. This can be problematic because bins that contain zero background events simply due tolimited simulation statistics can lead to sensitivity estimates that are too optimistic. To mitigatethis problem, we choose larger bins at higher energies: our bin width is for particle energiesbelow 10 GeV , for particle energies up to 20 GeV , and for particle energies up to 40 GeV . To be on the safe side, we also identify bins in which the background prediction is exactlyzero. We then consider averages of the background rate over increasingly larger neighborhoods ofeach such bin until we obtain a non-zero rate. If that average rate is above 0.1 events, we excludethe problematic bin from our analysis.As systematic uncertainties, we include a 10% normalization error, which we assume to beuncorrelated between different HNL decay channels and between the signal and background. Theseuncertainties are described in terms of nuisance parameters with Gaussian priors. To obtain thefinal sensitivity, our χ function is minimized over all nuisance parameters.Our results are shown in figs. 10 and 11. The four panels in fig. 10 correspond to differentassumption on the HNL couplings (coupling to ν e only, to ν µ only, to ν τ only, and to all threeactive neutrino flavors universally). We observe rich structure in these plots: first, we see that theprojected exclusion regions have the wedge shape that is typical for long-lived particle searches. Attoo small mixing, HNL production is suppressed beyond the detectable level and HNLs are so long-lived that they mostly decay far beyond the detector. At too large mixing, HNLs are abundantlyproduced, but most of them decay before reaching the detector. Away from kinematic thresholds,the upper limits on the mixing matrix elements | U e | , | U µ | , | U τ | (bottom edges of the wedge-shaped regions in fig. 10) scale roughly as | U α | limit ∝ M − . This can be understood as follows:9the HNL production rate scales as | U α | M over wide mass ranges, with the proportionalityto M reflecting the chiral suppression that otherwise affects many leptonic meson decays. TheHNL decay rate in the HNL rest frame scales as | U α | M , as can be seen from dimensionalanalysis. In the laboratory frame, this scaling changes to | U α | M 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 − . Overall, these argumentsshow that the experimental count rate scales as | U α | M , with deviations being observed close tokinematic thresholds. Moreover, the geometric scaling factor is not always exactly M − , dependingon 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 couplingonly (top left panel of fig. 10), we see that, at M ∼ m π ∼ 140 MeV , HNL production in piondecays 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 → νπ ) becomeallowed, 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 kaondecays becomes forbidden. Beyond this threshold, only HNLs from charm decay contribute, butsince charm production in DUNE/LBNF is relative inefficient due to the low primary protonenergy of 120 GeV , this leads to a significant drop in sensitivity. Nevertheless, the DUNE neardetectors will be able to probe a large chunk of parameter space up to M ∼ m D s ∼ , wherealso 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 couplingsexclusively to ν τ (bottom left panel of fig. 10), production thresholds do not play as important arole because charm decays are the only production mode at all HNL masses due to the requirementof producing a τ alongside the HNL. The small kink visible at M ∼ 190 MeV can be understoodfrom the D s – τ mass difference. Beyond the kink, HNL production in D s → N + τ decays iskinematically 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 only leads to signif-icant benefits in the N → νπ decay channel, which, however, never dominates the sensitivity.The reason for the larger benefit in this channel is the very large SM background. Here, the bettersignal-to-noise ratio in the off-axis locations, caused by the SM neutrino flux being more stronglyboosted in the forward direction than the flux of heavier HNLs, has a large impact. In otherdecay channels, backgrounds are not as big of a problem, but nevertheless going off-axis does notsignificantly harm the sensitivity for any of them. This means that the search for HNLs can becarried out truly parasitically to DUNE’s main oscillation program. No matter where the variousnear detectors are placed at any given moment, they will contribute in a useful way to the HNLsensitivity. Notably, the SAND (System for On-Axis Neutrino Detection) beam monitor, whichwill always remain on-axis and has not been included in our estimates, can be harnessed for theHNL 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 asin fig. 10: on-axis running only, the realistic DUNE-PRISM running strategy with equal amountsof data collected at seven different locations relative to the beam axis ( , , 12 m , 18 m , 24 m , 30 m , 36 m [48]), and a hypothetical background-free DUNE-PRISM search. The comparisonreveals 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 thisconclusion hinges to some extent on our treatment of charm production in LBNF, we are confidentthat, given the very conservative choice we have made for the charm production cross section, thesensitivity of the real experiment will likely be better than our estimate.)Comparing to planned experiments that may happen on a timescale similar to DUNE, we notethat, if any of these experiments is build, it will likely probe higher HNL masses than DUNE.0 HNL mass [MeV]10 − − − − − m x i n g | U e | . × pot, ν -mode 90% CL, 1 dofN– ν e mixing 10 HNL mass [MeV]10 − − − − − m x i n g | U µ | . × pot, ν -mode 90% CL, 1 dofN– ν µ mixing10 HNL mass [MeV]10 − − − − − m x i n g | U τ | 90% CL, 1 dofN– ν τ mixing PRISM on-axis no bg νπ e ± π ∓ µ ± π ∓ νe + e − νµ + µ − νe ± µ ∓ HNL mass [MeV]10 − − − − − m x i n g | U e | = | U µ | = | U τ | . × pot, ν -mode 90% CL, 1 dofflavor-universal mixing FIG. 10. Sensitivity to HNLs in the DUNE near detectors ND-LAr and ND-GAr as function of HNLmass. Solid lines correspond to the DUNE-PRISM running strategy with data taking both on-axis and atsix different off-axis locations ( , 12 m , 18 m , 24 m , 30 m , 36 m ), with equal exposure for each location.The total number of protons on target assumed here is . × , corresponding to 5 years of running inneutrino mode. Dotted lines show the sensitivity for on-axis data taking only, and thin dashed lines havebeen computed from a DUNE-PRISM analysis that neglects backgrounds. The four panels correspond todifferent HNL coupling structures, namely mixing with ν e only (top left panel), with ν µ only (top rightpanel), with ν τ only (bottom left panel), and flavor-universal mixing | U e | = | U µ | = | U τ | ≡ | U | (bottom right panel). This is not surprising because the experiments in questions (FASER-2, CODEX-b, NA62++,MATHUSLA, SHiP, FCC-ee) operate at higher beam energies. On the other hand, DUNE iswithout competition in terms of luminosity. Therefore, at HNL masses below the kaon threshold,DUNE’s sensitivity is unchallenged and will remain so for the foreseeable future.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 chargeidentification capabilities. V Summary and Conclusions In this work, we have explored the benefits of the DUNE-PRISM detectors for exploring physicsbeyond 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 movein and out of the beam axis and have found that exploiting this capability is highly beneficial in1 − HNL mass [GeV]10 − − − − − − m i x i n g | U e | BB N S ee s a w F CC - ee 90% CL, 1 dof5 . × pot, ν -mode CMSDELPHIEWPD B e ll e S H i P ( n o B c ) M A T HU S L A P S C HA R M F A S E R - C O D E X - b NA ++ N– ν e mixing DUNE-PRISMon-axis onlyno background10 − HNL mass [GeV]10 − − − − − − m i x i n g | U e | BB N S ee s a w F CC - ee 90% CL, 1 dof5 . × pot, ν -mode CMSDELPHIEWPD B e ll e S H i P ( n o B c ) M A T HU S L A P S N u T e V L H C b E C HA R M F A S E R - C O D E X - b NA ++ N– ν µ mixing DUNE-PRISMon-axis onlyno background10 − HNL mass [GeV]10 − − − − − − m i x i n g | U e | BB N S ee s a w F CC - ee 90% CL, 1 dof5 . × pot, ν -mode DELPHI S H i P ( n o B c ) M A T HU S L A C HA R M F A S E R - C O D E X - b NA ++ N– ν τ mixing DUNE-PRISMon-axis onlyno background FIG. 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 e | , | U µ | , and | U τ | ,respectively, assuming that only one of them is non-zero at a time. In this plot, the decay channels νe + e − , νµ + µ − , νe ± µ ∓ , νπ , π ± e ∓ , π ± µ ∓ are combined. As in fig. 10, solid blue curves show the sensitivityof DUNE-PRISM ( . × pot in neutrino mode, equally split between the on-axis position and sixdifferent off-axis locations). Dotted blue contours show results for on-axis running only, and thin dashedblue curves represent a hypothetical background-free analysis. The contours shown in the backgroundcorrespond to existing limits (filled) and sensitivities of planned experiments (unfilled dashed and dotted).These limits are taken from the compilation in ref. [75]. . × pot ) will be able to improve existing limits on darkphoton-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 (cid:48) gauge boson rather than a dark photon,we are in the interesting situation that, for some parameter ranges (namely those with relativelyheavy Z (cid:48) ), the DM flux increases away from the beam axis, while the background from SMneutral current neutrino interactions falls off. Consequently, the DUNE-PRISM strategy is alwaysadvantageous compared to an on-axis-only run, see fig. 6. After 5 years, DUNE-PRISM limitswill 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 bothon-axis and off-axis as in DUNE-PRISM never hurts the sensitivity. For some decay channels,especially those with large backgrounds like N → νπ , off-axis running is highly beneficial. Forchannels where backgrounds are lower or can be well distinguished from the signal by using spectralinformation, 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 existinglimits for heavy neutral lepton mixing with ν e and ν µ . Significant new territory will be coveredfor 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 numerousextensions of the SM. Their unique capability of moving off-axis will only improve the sensitivityto such scenarios, implying that a rich program of new physics searches can be carried out concur-rently with DUNE’s neutrino oscillation program without requiring an adaptation of the runningstrategy. Acknowledgments We are deeply grateful to Tommaso Boschi for answering numerous questions about his NuShock code [68], to Laura Fields for useful discussions about the DUNE fluxes, and to PedroMachado for helping us in the comparison of our results to the results of ref. [13]. We havebenefited from useful discussions with Gaia Lanfranchi, Albert de Roeck, and the members of theDUNE-HNL working group. This work has been partly funded by the German Research Foun-dation (DFG) in the framework of the PRISMA+ Cluster of Excellence and by the EuropeanResearch Council (ERC) under the European Union’s Horizon 2020 research and innovation pro-gramme (grant agreement No. , “ ν Directions”). The research of L. B. is supported in partby the Swiss National Foundation under Contracts No. and No. IZSAZ2_173357 . 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 signaland background predictions with those of ref. [13]. We have found that our estimate of thebackground due to neutrino–electron scattering agrees well with the results of ref. [13] in both theneutrino- and anti-neutrino mode. In the following, we focus on the estimate of the number ofsignal events.As discussed in sec. II.A a crucial aspect in the simulation of the signal events is the modeling ofthe light meson spectra from the DUNE/LBNF target. In ref. [13], only the production of mesons3 e v e n t s / po s iti on / y ea r bkgsig+bkg primary-only (SoftQCD)sig+bkg beam-dumpsig+bkg refsig+bkg primary-only (HardQCD) m A ' =90 MeV = m φ , α D ε =10 -15 neutrino-mode ( s i g + bkg ) / bkg Δ x OA [m] FIG. 12. Same as fig. 1, but including the result of ref [13] in solid green, and our own primary-only(HardQCD) sample in dotted orange. in the primary proton interaction was considered, and Pythia was used as the main simulationtool. In fig. 12 we compare our results with those of fig. 3 of ref. [13] using the same benchmarkpoint m A (cid:48) = 90 MeV = 3 m φ , (cid:15) α D = 10 − . Contrary to our expectation, we observe that our primary-only (SoftQCD) curve fails to accurately reproduce the results of that paper (solid greenline) at any value of ∆ x OA . As anticipated in section II.A, the source of the discrepancy is thatthe authors of ref. [13] have used a different configuration for their Pythia simulation. Moreprecisely, they activated the HardQCD:All flag, and if we do the same (dotted orange curve infig. 12), we are indeed able to reproduce their results.We remark that using the SoftQCD:All flag is preferable for modeling the proton–protonprimary interactions for the case of relatively low energy fixed target experiments such asDUNE/LBNF. (This was shown in ref. [35], which appeared after ref. [13].) 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 coincidentallysimilar 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 pre-dicted signal and background rates in the analyses of sections II to IV in more detail. As stated inthese sections we use standard frequentist techniques to test the signal + background hypothesisagainst simulated background-only data. We employ a Poissonian log-likelihood function, i.e. − L ( Θ , X ) ≡ n pos (cid:88) j =1 n bins (cid:88) i =1 (cid:20) B ij ( X ) + S ij ( Θ , X ) − B ij (0) + B ij (0) log (cid:18) B ij (0) B ij ( X ) + S ij ( Θ , X ) (cid:19)(cid:21) + (cid:88) c = S,B (cid:40)(cid:18) X norm c, correl σ correl (cid:19) + (cid:18) X tilt c, correl σ correl (cid:19) + n pos (cid:88) j =1 (cid:20)(cid:18) X norm c,j σ pos (cid:19) + (cid:18) X tilt c,j σ pos (cid:19) (cid:21)(cid:41) . (B1)4 -5 -4 -3 -2 -1 mrad < θ OA < mrad / N d N π / d E [ G e V - ] E [GeV] mrad < θ OA < mradE [GeV] mrad < θ OA < mradE [GeV]SoftQCDbeam-dumpHardQCD FIG. 13. Energy spectrum of π s produced in the beam dump in three angular windows: < θ OA < mrad (left panel), mrad < θ OA < mrad (middle panel) and mrad < θ OA < mrad (right panel).The three curves correspond to the three meson samples introduced in the main text: primary-only (Soft-QCD) in red dot-dashed, beam-dump in blue dashed primary-only (HardQCD) in orange dotted. The signal and background rates in the i -th energy bin at the j -th off-axis position includingsystematic biases, S ij ( Θ , X ) and B ij ( X ) , are defined in terms of the rates without biases, S ij ( Θ , and B ij (0) , according to S ij ( Θ , X ) = (cid:2) X norm S, correl (cid:3)(cid:2) X norm S,j (cid:3)(cid:2) X tilt S, correl ( − , . . . , i (cid:3)(cid:2) X tilt S,j ( − , . . . , i (cid:3) S ij ( Θ , ,B ij ( X ) = (cid:2) X norm B, correl (cid:3)(cid:2) X norm B,j (cid:3)(cid:2) X tilt B, correl ( − , . . . , i (cid:3)(cid:2) X tilt B,j ( − , . . . , i (cid:3) B ij (0) . (B2)We collectively denote the vector of physical model parameters Θ and the vector of nuisanceparameters X . The nuisance parameters X norm S, corr and X norm S,j describe systematic normalizationuncertainties 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 correspondingparameters affecting the background ( X norm B, corr , X norm B,j , X tilt B, correl and X tilt B,j ) is analogous. (Notethat tilt errors are not considered in the HNL analysis of section IV.) All systematic errors aretreated as Gaussian and are constrained by the pull terms in the second row of eq. (B1). Assystematic normalization and tilt uncertainties we assume σ pos = 1% (uncorrelated between dif-ferent on-/off-axis positions) and σ correl = 10% (correlated among all positions). The sums inthe first line of eq. (B1) run over n bins = 80 ( n bins = 57 ) energy bins, equally spaced in theinterval [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 ( , / . 45 mrad , 12 m / . 90 mrad , 18 m / . 36 mrad , 24 m / . 81 mrad , 30 m / . 26 mrad , 36 m / . 72 mrad ). For comparison we also present results foran 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 estimatedfrom the log-likelihood ratio [76] Z ( Θ ) ≡ − (cid:18) L ( Θ , ˆˆ X ) L ( ˆ Θ , ˆ X ) (cid:19) , (B3)where ( ˆ Θ , ˆ X ) is the combination of model parameters and nuisance parameters that maximizesthe likelihood (minimizes χ ), and ˆˆ X are the nuisance parameters that maximize the likelihoodfor fixed Θ . Z ( Θ ) follows a χ distribution, with the number of degrees of freedom equal to thedimension of Θ . Consequently, the nσ exclusion region can be estimated by comparing Z ( Θ ) tothe nσ quantile of that χ distribution.5However, experimental collaboration often present their results as constraints on the “signalstrength” µ , which is defined as an overall, energy-independent rescaling factor of the event raterelative to some reference point. In the dark photon model, the signal strength can be taken asthe product of coupling constants µ ≡ (cid:15) α D , while in the leptophobic model it is µ ≡ g Z . In thisway of analyzing the data, Z ( Θ ) is replaced by a one-parameter function Z ( µ ) ≡ − (cid:18) L ( µ, ˆˆ X ) L (ˆ µ, ˆ X ) (cid:19) , (B4)whose values are then compared to the χ 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 signalstrength analysis, however, only µ is varied and all other model parameters are kept fixed. Thestatistical question answered by a signal strength analysis is “What is the constraint on the signalstrength, assuming we already know all other model parameters.” A multi-parameter fit, on theother hand, asks “What are the preferred parameter regions, assuming the chosen model is thecorrect one, but we do not have any prior information on any of its parameters.” To make ourresults comparable to those shown in the literature, we use signal strength analyses throughout. [1] B. Batell, M. Pospelov, and A. Ritz, Exploring Portals to a Hidden Sector Through Fixed Targets , Phys. Rev. D (2009) 095024, [ arXiv:0906.5614 ].[2] P. deNiverville, M. Pospelov, and A. Ritz, Observing a light dark matter beam with neutrinoexperiments , Phys.Rev. D84 (2011) 075020, [ arXiv:1107.4580 ].[3] P. deNiverville, D. McKeen, and A. Ritz, Signatures of sub-GeV dark matter beams at neutrinoexperiments , Phys.Rev. D86 (2012) 035022, [ arXiv:1205.3499 ].[4] MiniBooNE Collaboration, R. Dharmapalan et al., Low Mass WIMP Searches with a NeutrinoExperiment: A Proposal for Further MiniBooNE Running , arXiv:1211.2258 .[5] B. Batell, P. deNiverville, D. McKeen, M. Pospelov, and A. Ritz, Leptophobic Dark Matter atNeutrino Factories , Phys. Rev. D (2014), no. 11 115014, [ arXiv:1405.7049 ].[6] D. E. Soper, M. Spannowsky, C. J. Wallace, and T. M. P. Tait, Scattering of Dark Particles withLight Mediators , Phys. Rev. D (2014), no. 11 115005, [ arXiv:1407.2623 ].[7] B. A. Dobrescu and C. Frugiuele, GeV-scale dark matter: production at the Main Injector , JHEP (2015) 019, [ arXiv:1410.1566 ].[8] P. Coloma, B. A. Dobrescu, C. Frugiuele, and R. Harnik, Dark matter beams at LBNF , JHEP (2016) 047, [ arXiv:1512.03852 ].[9] C. Frugiuele, Probing sub-GeV dark sectors via high energy proton beams at LBNF/DUNE andMiniBooNE , Phys. Rev. D96 (2017), no. 1 015029, [ arXiv:1701.05464 ].[10] P. deNiverville and C. Frugiuele, Hunting sub-GeV dark matter with the NO ν A near detector , Phys.Rev. D (2019), no. 5 051701, [ arXiv:1807.06501 ].[11] G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai, Millicharged particles in neutrino experiments , arXiv:1806.03310 .[12] A. de Gouvêa, P. J. Fox, R. Harnik, K. J. Kelly, and Y. Zhang, Dark Tridents at Off-Axis LiquidArgon Neutrino Detectors , JHEP (2019) 001, [ arXiv:1809.06388 ].[13] V. De Romeri, K. J. Kelly, and P. A. N. Machado, Hunting On- and Off-Axis for Light Dark Matterwith DUNE-PRISM , arXiv:1903.10505 .[14] L. Buonocore, C. Frugiuele, and P. deNiverville, Hunt for sub-GeV dark matter at neutrinofacilities: A survey of past and present experiments , Phys. Rev. D (2020), no. 3 035006,[ arXiv:1912.09346 ].[15] B. Batell, J. Berger, and A. Ismail, Probing the Higgs Portal at the Fermilab Short-BaselineNeutrino Experiments , Phys. Rev. D (2019), no. 11 115039, [ arXiv:1909.11670 ].[16] P. Ballett, T. Boschi, and S. Pascoli, Heavy Neutral Leptons from low-scale seesaws at the DUNENear Detector , JHEP (2020) 111, [ arXiv:1905.00284 ].[17] M. Battaglieri et al., US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report , arXiv:1707.04591 .[18] SENSEI Collaboration, M. Crisler, R. Essig, J. Estrada, G. Fernandez, J. Tiffenberg, M. Sofo haro, T. Volansky, and T.-T. Yu, SENSEI: First Direct-Detection Constraints on sub-GeV Dark Matterfrom a Surface Run , arXiv:1804.00088 .[19] MiniBooNE Collaboration, A. A. Aguilar-Arevalo et al., Dark Matter Search in a Proton BeamDump with MiniBooNE , Phys. Rev. Lett. (2017), no. 22 221803, [ arXiv:1702.02688 ].[20] MiniBooNE DM Collaboration, A. A. Aguilar-Arevalo et al., Dark Matter Search in Nucleon,Pion, and Electron Channels from a Proton Beam Dump with MiniBooNE , Phys. Rev. D98 (2018),no. 11 112004, [ arXiv:1807.06137 ].[21] The DUNE Collaboration, The DUNE Near Detector Conceptual Design Report - Final Draft , .available from https://docs.dunescience.org/cgi-bin/ShowDocument?docid=21267 .[22] E. Izaguirre, G. Krnjaic, P. Schuster, and N. Toro, Analyzing the Discovery Potential for Light DarkMatter , Phys. Rev. Lett. (2015), no. 25 251301, [ arXiv:1505.00011 ].[23] T. Lin, H.-B. Yu, and K. M. Zurek, On Symmetric and Asymmetric Light Dark Matter , Phys. Rev.D (2012) 063503, [ arXiv:1111.0293 ].[24] Planck Collaboration, P. A. R. Ade et al., Planck 2015 results. XIII. Cosmological parameters , Astron. Astrophys. (2016) A13, [ arXiv:1502.01589 ].[25] T. R. Slatyer, Indirect dark matter signatures in the cosmic dark ages. I. Generalizing the bound ons-wave dark matter annihilation from Planck results , Phys. Rev. D (2016), no. 2 023527,[ arXiv:1506.03811 ].[26] T. R. Slatyer, Indirect Dark Matter Signatures in the Cosmic Dark Ages II. Ionization, Heating andPhoton Production from Arbitrary Energy Injections , Phys. Rev. D (2016), no. 2 023521,[ arXiv:1506.03812 ].[27] P. deNiverville and C. Frugiuele, Hunting sub-GeV dark matter with the NO ν A near detector , Phys.Rev. D (2019), no. 5 051701, [ arXiv:1807.06501 ].[28] L. Buonocore, C. Frugiuele, and P. deNiverville, Hunt for sub-GeV dark matter at neutrinofacilities: A survey of past and present experiments , Phys. Rev. D (2020), no. 3 035006,[ arXiv:1912.09346 ].[29] A. Celentano, L. Darmé, L. Marsicano, and E. Nardi, New production channels for light dark matterin hadronic showers , Phys. Rev. D (2020), no. 7 075026, [ arXiv:2006.09419 ].[30] A. Berlin, S. Gori, P. Schuster, and N. Toro, Dark Sectors at the Fermilab SeaQuest Experiment , Phys. Rev. D98 (2018), no. 3 035011, [ arXiv:1804.00661 ].[31] SHiP Collaboration, C. Ahdida et al., Sensitivity of the SHiP experiment to light dark matter , arXiv:2010.11057 .[32] S. Gardner, R. J. Holt, and A. S. Tadepalli, New Prospects in Fixed Target Searches for Dark Forceswith the SeaQuest Experiment at Fermilab , Phys. Rev. D93 (2016), no. 11 115015,[ arXiv:1509.00050 ].[33] Y. Kahn, G. Krnjaic, J. Thaler, and M. Toups, DAE δ ALUS and dark matter detection , Phys. Rev. D91 (2015), no. 5 055006, [ arXiv:1411.1055 ].[34] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O.Rasmussen, and P. Z. Skands, An Introduction to PYTHIA 8.2 , Comput. Phys. Commun. (2015) 159–177, [ arXiv:1410.3012 ].[35] J. M. Berryman, A. de Gouvea, P. J. Fox, B. J. Kayser, K. J. Kelly, and J. L. Raaf, Searches forDecays of New Particles in the DUNE Multi-Purpose Near Detector , JHEP (2020) 174,[ arXiv:1912.07622 ].[36] GEANT4 Collaboration, S. Agostinelli et al., GEANT4–a simulation toolkit , Nucl. Instrum. Meth.A (2003) 250–303.[37] L. Buonocore, C. Frugiuele, F. Maltoni, O. Mattelaer, and F. Tramontano, Event generation forbeam dump experiments , JHEP (2019) 028, [ arXiv:1812.06771 ].[38] C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer, and T. Reiter, UFO - The UniversalFeynRules Output , Comput. Phys. Commun. (2012) 1201–1214, [ arXiv:1108.2040 ].[39] E. Fermi, On the Theory of the impact between atoms and electrically charged particles , Z. Phys. (1924) 315–327.[40] E. J. Williams, Nature of the high-energy particles of penetrating radiation and status of ionizationand radiation formulae , Phys. Rev. (1934) 729–730.[41] C. F. von Weizsäcker, Radiation emitted in collisions of very fast electrons , Z. Phys. (1934)612–625.[42] J. Blümlein and J. Brunner, New Exclusion Limits on Dark Gauge Forces from ProtonBremsstrahlung in Beam-Dump Data , Phys. Lett. B731 (2014) 320–326, [ arXiv:1311.3870 ].[43] P. deNiverville, C.-Y. Chen, M. Pospelov, and A. Ritz, Light dark matter in neutrino beams:Production modeling and scattering signatures at miniboone, t2k, and ship , Phys. Rev. D (2017) Once more on electromagnetic formfactors of nucleons in extended vector meson dominance model , Phys. Rev. C82 (2010) 038201,[ arXiv:0910.5589 ].[45] D. E. Morrissey and A. P. Spray, New Limits on Light Hidden Sectors from Fixed-TargetExperiments , Journal of High Energy Physics (2014) 083, [ arXiv:1402.4817 ].[46] D. Gorbunov, A. Makarov, and I. Timiryasov, Decaying light particles in the SHiP experiment:Signal rate estimates for hidden photons , Phys. Rev. D91 (2015), no. 3 035027, [ arXiv:1411.4007 ].[47] C. Andreopoulos et al., The GENIE Neutrino Monte Carlo Generator , Nucl. Instrum. Meth. A (2010) 87–104, [ arXiv:0905.2517 ].[48] L. F. (for the DUNE Collaboration), “Simulated dune neutrino fluxes.” ROOT Ntuple files availableat https://home.fnal.gov/~ljf26/DUNEFluxes/ . We use the v3r5p4 release which has also beenused for the DUNE TDR.[49] DUNE Collaboration, R. Acciarri et al., Long-Baseline Neutrino Facility (LBNF) and DeepUnderground Neutrino Experiment (DUNE): Conceptual Design Report, Volume 2: The PhysicsProgram for DUNE at LBNF , arXiv:1512.06148 .[50] BaBar Collaboration, J. Lees et al., Search for Invisible Decays of a Dark Photon Produced in e + e − Collisions at BaBar , Phys. Rev. Lett. (2017), no. 13 131804, [ arXiv:1702.03327 ].[51] D. Banerjee et al., Dark matter search in missing energy events with NA64 , Phys. Rev. Lett. (2019), no. 12 121801, [ arXiv:1906.00176 ].[52] J. Beacham et al., Physics Beyond Colliders at CERN: Beyond the Standard Model Working GroupReport , J. Phys. G47 (2020), no. 1 010501, [ arXiv:1901.09966 ].[53] B. Batell, R. Essig, and Z. Surujon, Strong Constraints on Sub-GeV Dark Sectors from SLAC BeamDump E137 , Phys. Rev. Lett. (2014), no. 17 171802, [ arXiv:1406.2698 ].[54] B. A. Dobrescu and C. Frugiuele, Hidden GeV-scale interactions of quarks , Phys.Rev.Lett. (2014) 061801, [ arXiv:1404.3947 ].[55] J. A. Dror, R. Lasenby, and M. Pospelov, Dark forces coupled to nonconserved currents , Phys. Rev. D96 (2017), no. 7 075036, [ arXiv:1707.01503 ].[56] L. Michaels and F. Yu, Probing new U (1) gauge symmetries via exotic Z → Z (cid:48) γ decays , arXiv:2010.00021 .[57] NA61/SHINE Collaboration, N. Abgrall et al., Measurements of Cross Sections and Charged PionSpectra in Proton-Carbon Interactions at 31 GeV/c , Phys. Rev. C (2011) 034604,[ arXiv:1102.0983 ].[58] Particle Data Group Collaboration, M. Tanabashi et al., Review of Particle Physics , Phys. Rev.D (2018), no. 3 030001.[59] R. D. Ball et al., Parton distributions with LHC data , Nucl. Phys. B867 (2013) 244–289,[ arXiv:1207.1303 ].[60] NNPDF Collaboration, R. D. Ball, V. Bertone, S. Carrazza, L. Del Debbio, S. Forte, A. Guffanti,N. P. Hartland, and J. Rojo, Parton distributions with QED corrections , Nucl. Phys. B877 (2013)290–320, [ arXiv:1308.0598 ].[61] MiniBooNE, MINOS Collaboration, P. Adamson et al., First Measurement of ν µ and ν e Eventsin an Off-Axis Horn-Focused Neutrino Beam , Phys. Rev. Lett. (2009) 211801,[ arXiv:0809.2447 ].[62] M. L. Graesser, I. M. Shoemaker, and L. Vecchi, A Dark Force for Baryons , arXiv:1107.2666 .[63] BaBar Collaboration, B. Aubert et al., A Search for Invisible Decays of the Upsilon(1S) , Phys.Rev. Lett. (2009) 251801, [ arXiv:0908.2840 ].[64] A. Atre, T. Han, S. Pascoli, and B. Zhang, The Search for Heavy Majorana Neutrinos , JHEP (2009) 030, [ arXiv:0901.3589 ].[65] M. Drewes and B. Garbrecht, Combining experimental and cosmological constraints on heavyneutrinos , Nucl. Phys. B (2017) 250–315, [ arXiv:1502.00477 ].[66] A. de Gouvêa and A. Kobach, Global Constraints on a Heavy Neutrino , Phys. Rev. D (2016),no. 3 033005, [ arXiv:1511.00683 ].[67] P. D. Bolton, F. F. Deppisch, and P. Bhupal Dev, Neutrinoless double beta decay versus other probesof heavy sterile neutrinos , JHEP (2020) 170, [ arXiv:1912.03058 ].[68] T. Boschi, Nushock: Bunch of tools to study sensitivity of heavy neutrinos , 2019. available onGitHub: https://github.com/tboschi/NuShock , repository cloned on September 15, 2020.[69] C. Lourenco and H. Wohri, Heavy flavour hadro-production from fixed-target to collider energies , Phys. Rept. (2006) 127–180, [ hep-ph/0609101 ].[70] P. Ramana Murthy, C. A. Ayre, H. Gustafson, L. W. Jones, and M. J. Longo, Neutron Total Cross-Sections on Nuclei at Fermilab Energies , Nucl. Phys. B (1975) 269–308.[71] ZEUS Collaboration, H. Abramowicz et al., Measurement of charm fragmentation fractions inphotoproduction at HERA , JHEP (2013) 058, [ arXiv:1306.4862 ].[72] S. Aoki et al., Study of tau-neutrino production at the CERN SPS , arXiv:1708.08700 .[73] Dawson-Haggerty et al., “trimesh.” available from https://trimsh.org/ .[74] P. Ballett, M. Hostert, S. Pascoli, Y. F. Perez-Gonzalez, Z. Tabrizi, and R. Zukanovich Funchal, Neutrino Trident Scattering at Near Detectors , JHEP (2019) 119, [ arXiv:1807.10973 ].[75] J. Beacham et al., Physics Beyond Colliders at CERN: Beyond the Standard Model Working GroupReport , J. Phys. G (2020), no. 1 010501, [ arXiv:1901.09966 ].[76] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Asymptotic formulae for likelihood-based tests ofnew physics , Eur. Phys. J. C71 (2011) 1554, [ arXiv:1007.1727arXiv:1007.1727