Prospects for Detecting Boosted Dark Matter in DUNE through Hadronic Interactions

Boosted dark matter (BDM) is a well-motivated class of dark matter (DM) candidates in which a small component of DM is relativistic at the present time. We lay the foundation for BDM searches via hadronic interactions in large liquid-argon time-projection chambers (LArTPCs), such as DUNE. We investigate BDM-nucleus scattering in detail by developing new event generation techniques with a parameterized detector simulation. We study the discovery potential in a DUNE-like experiment using the low threshold and directionality of hadron detection in LArTPCs and compare with other experiments.


Introduction.
Despite the overwhelming gravitational evidence for the existence of dark matter (DM), its microscopic nature remains a profound puzzle. A leading DM paradigm is that of Weakly Interacting Massive Particles (WIMP), consisting of a single species of deeply nonrelativistic particles. Over the past few decades, however, DM detection experiments have excluded large swaths of the parameter space for WIMPs [1][2][3][4][5][6][7][8], motivating serious consideration of non-minimal models and alternative candidates.
In a class of models beyond the minimal WIMP scenario, a small relativistic component of DM, boosted dark matter (BDM) [9][10][11], is produced at the present time and can be detected via its interactions with the Standard Model (SM) particles. The detection of BDM could be a smoking gun for DM discovery in cases where the dominant component of DM is hard to detect, yet it requires new experimental strategies beyond the current DM searches.
BDM may originate from scenarios of dark sectors, with multiple components of DM or with non-minimal stabilization mechanisms, such as semi-annihilating DM [11,12], selfannihilating DM [13,14], decaying DM [15,16], DM induced nucleon decay [9,17], or cosmic ray acceleration [18][19][20][21]. A minimal two-component scenario described in the original works [9][10][11] includes a cold component ψ as the dominant component of DM with very small scattering cross sections with the SM particles, and a relativistic, less massive secondary component χ produced by the annihilation of ψ, that effectively interacts with the SM particles. Thermal freezeout via processes such as ψψ → χχ annihilation in the early Universe may determine DM relic abundance as a new realization of WIMP miracle [10,11,22]. Meanwhile, present-day annihilation in DM-concentrated regions, such as the Galactic Center (GC) or the Sun, generates BDM χ that can be detected via their interaction with electrons or hadrons.
The phenomenology of BDM features a relatively small flux and typically (semi-)relativistic outgoing SM particles upon BDM-SM particle scattering. As conventional DM direct detection experiments focus on the detection of lowenergy recoils in a detector mass up to a few tons, they generally do not have the best sensitivity for BDM searches (with the exception of very low mass DM, see [16,[23][24][25]). On the other hand, massive neutrino detectors, sensitive to energetic SM particles, stand out as ideal facilities.
Experimentally, BDM can be observed in interactions with electrons or with nuclei. In order to cover all possibilities, both kinds of interactions should be studied. While previous work has extensively studied the parameter space probed by BDM models when χ scatters off electrons [10,[26][27][28][29][30][31][32][33], BDM detection via interaction with nuclei, which is complicated by nuclear effects, is much less studied. Nuclear scattering is, however, the dominant process in numerous well-motivated models [34][35][36][37][38][39][40][41][42], making this interaction a potential discovery channel. It is possible that leptonic interactions are not even present at all.
We expect the main background for BDM interactions in a detector to be from atmospheric neutrinos interacting via the neutral current, leaving activity in an energy range similar to the signal. Unlike this background, all the BDM signal comes from a single source, and the signal contribution can therefore be enhanced by selecting events aligning with the source's location. This is done by selecting events based on the total momentum of all the detectable particles produced in the interaction.
The exploration and development of the novel technology of liquid-argon time-projection chambers (LArTPCs) as neutrino detectors has ramped up in the last decade, and will culminate in the upcoming next generation neutrino facility, Deep Underground Neutrino Experiment (DUNE) [19, 23, 27-29, 32, 33, 43-50]. Recent studies have shown that search for BDM via interactions with electrons would benefit from DUNE's excellent particle identification [27-29, 47, 51-63]. LArTPCs are, however, expected to most significantly im-prove the sensitivity to BDM in the hadronic channels. The accessible kinematic range for hadrons in water Cherenkov detectors is limited by the Cherenkov threshold (a momentum of 1.07 GeV for protons) and, in the case of inelastic scattering, by the quality of the reconstruction of overlapping rings [11,64]. Massive detectors based on liquid scintillators do not provide directionality, and segmented liquid scintillators, which can offer the details of an event such as directionality, have either a small volume or a relatively coarse granularity, due to cost. LArTPCs have millimeter resolution, leading to a low detection threshold of hadrons and an ability to reconstruct recoil direction. They are scalable and have excellent capabilities in calorimetry and thereby particle identification. These features combine to alleviate the limitations of current experiments. Multi-kiloton LArTPC experiments like DUNE, which features a fiducial mass of 40 kilotons, therefore hold great potential for BDM searches.
Liquid argon, a dense fluid with moderately large nuclei, is a consummate target candidate of BDM detection through hadronic scattering, granting higher interaction rates. Conversely, the nuclear effects of argon, and in particular the propagation and interaction through the nucleus (known as final-state interactions, or FSI) of the hadrons produced in the BDM-nucleon interactions, alter the kinematics of the particles revealed in the detector [44,65], jeopardizing the reconstruction, and in particular that of the direction, of the signal candidates.
In this Letter, we study the observed hadronic signatures of BDM in LArTPC detectors evaluating the nuclear effects using a novel Monte Carlo (MC)-based analysis, and obtain the search sensitivity taking into consideration the atmospheric neutrino background.
We consider the following representative BDM model as a benchmark for our study [11]. The model consists of two components of DM. The dominant component, ψ, scatters off of hydrogen in the Sun, gets captured and builds up, then annihilates into the relativistic lighter component, χ, i.e. the BDM: (1) The modeling of this annihilation is not particularly relevant to the phenomenology at hand, but we assume that this is the dominant annihilation process for ψ. As we discuss shortly, so long as the annihilation cross section is sufficiently large, it will not enter into the determination of the BDM flux. Although annihilation is also possible in the Galactic Center, for a broad range of parameters the flux from the Sun will dominate over that from the Galactic Center. The large solar flux makes it possible to have observable signals with scattering cross sections of weak scale size or even smaller. Equilibrium between DM capture and DM annihilation is generically reached in the solar core [9,11], eliminating the parametric dependence on the DM annihilation cross section. As a minimal assumption, we do not introduce leptonic interactions, although the model can accommodate both leptonic and hadronic couplings as independent interactions. More details of this model as well as of semi-annihilation scenarios can be found in [11].
The BDM χ produced in the above process then emerges from the Sun at high velocity and scatters off of nuclei in the detector, via a process of the form where N is a nucleon and X is any number of hadrons. Hadronic DM interactions share similarities with neutral current neutrino scattering, which makes it natural to perform simulations in the framework of neutrino MC software suites, such as GENIE [65][66][67]. For this analysis we introduce in GE-NIE a BDM module [67] to perform all of the cross section calculations and event generation. GENIE simulates several nuclear effects, such as nucleon motion, Pauli blocking, and final state interactions of hadrons escaping the nuclear remnant after scattering. It further includes parton distributions, fragmentation, and hadronization in deep inelastic scattering, with some corrections to deal with the relatively low energy regime of interest. The hA final state model [66,68,69] is employed to model nuclear effects, though it can be changed to compare with other models. In the energy regime relevant to LArTPC neutrino detectors, E ν > ∼ 100 MeV, coherent nuclear scattering is highly suppressed and scattering is dominantly off nucleons, that become unbound from the nucleus, or electrons, that have negligible binding energies compared to the momentum transfer.
We include elastic scattering, yielding a recoiling nucleon, and deep inelastic scattering (DIS), yielding multi-hadron final states, in our modeling, while conservatively neglecting resonant inelastic scattering processes during which an excited baryon is produced and decays [95]. The diagrams of all the three processes are shown in Fig. 3 in the Supplementary Material. Elastic scattering off nucleons can be described by an axial form factor. As for neutrino scattering in GENIE, the axial form factor is assumed to have a dipole form. The normalization of this form factor is given by the spin form factors, which are currently best determined by lattice QCD calculations [70]. The hadronic component of the DIS scattering cross section is described by a hadronic tensor, which depends on parton distribution functions (PDF). GENIE uses a PDF that includes corrections for the relatively low energy regime of interest. The fragmentation and hadronization of the final state depends on the invariant mass of the final state hadronic system. At low invariant masses, an empirical model is used [71], in which we assume that DM scattering is similar to neutrino scattering. At high invariant masses, a model based on PYTHIA [72] is used. Further details can be found in the Supplementary Material.
Analysis strategy.
For concreteness, we focus on a benchmark in which both components of DM are scalars. Both are required to interact with quarks in order to enable solar capture for the heavy component, and terrestrial detection for the light component. The interactions with the SM particles are mediated by a spin-1 vector boson, Z , with a gauge coupling g Z . We assume that the quark current is axial. Both DM species, as well as the SM quarks, have a charge under this interaction, which is a free parameter of the model. Without loss of generality we take the BDM charge Q χ = 1. As a simple benchmark, we take the quark charges Q f = 1 for all quark flavors, and consider m Z = 1 GeV. For m Z > ∼ 1 GeV, the effect of the Z on the BDM scattering kinematics is small. We leave the heavy DM charge Q ψ and gauge coupling g Z , as well as the masses of the DM species m ψ , m χ as free parameters. Note that the lighter DM emerges from the Sun with a Lorentz boost Even with a mild hierarchy of masses, the velocity of χ from the Sun can be much larger than that of the virialized DM in the Solar System, such that χ can escape and reach Earth as BDM.
Within this model, we determine the flux of BDM χ through a detector on Earth. The flux depends on three sequential processes: capture, annihilation, and rescattering. The DM capture rate in the model considered has been determined in [11] and we use a similar calculation. Given a heavy DM mass, m ψ > ∼ 4 GeV, and a large enough annihilation cross section, (σv) ann > ∼ 3 × 10 26 cm 3 /sec, DM loss through annihilation and DM gain through capture reach the equilibrium within the lifetime of the Sun over the entire parameter space to which a multi-kiloton LArTPC is sensitive. In this case, BDM flux is simply determined by the DM capture rate. In addition, direct detection experiments and Super-Kamiokande exclude the parameter region in which BDM rescatters as it exits the Sun [11]. Rescattering is thus negligible for the parameter range of interest for this study, leading to a nearly monoenergetic flux of χ. Combining the above processes, we find that the magnitude of the flux where C is the ψ capture rate, proportional to g 4 Z , and D is the distance from the Sun to the Earth, i.e. 1 AU.
We scan over the parameter space of four BDM masses m χ in the range of (5 -40) GeV and three boost factors γ = 1.1, 1.5, 10.0, while probing the coupling constant g Z . For a mass m ψ below 5 GeV, evaporation of captured dark matter would lead to drastically reduced flux on the Earth, while above 40 GeV, the DM mass no longer has a significant effect on the detection efficiency. Our three choices of γ, in order, represent the benchmark cases where the BDM-hadron interaction is all elastic scattering, a mixture of elastic and inelastic scattering, and mostly inelastic scattering.
We use GENIE to simulate the BDM signals and the atmospheric neutrino background, where the BDM signal simulation, discussed above, is developed for this analysis. The direction of the Sun with respect to the detector, evaluated based on the SolTrack package [73] and on the geographical coordinates of DUNE [74] as an example, is encoded in the samples. Based on the Bartol atmospheric neutrino flux [75] at Soudan, the atmospheric neutrino samples include neutralcurrent (NC) neutrino events and events where ν τ interacts with the detector target via charged current (CC) and the outgoing τ leptons decay into hadrons. The rest of the CC neutrino interactions is not included in the background sample, as we assume with the information offered by LArTPCs, it can be efficiently rejected by discarding the events which contain muons or electrons as final state particles [96]. Charged particles produced in the χ-Ar interactions, as well as in the propagation of neutral particles in liquid argon, are the observables in LArTPCs [44,76]. The four-momenta of the stable final-state particles, including protons, electrons, photons, and charged pions, are convolved with the detector resolution reported in the DUNE Conceptual Design Report (CDR) [77], while the ones with kinetic energy below the detection threshold in DUNE CDR and all the neutrons are excluded. The distribution of cos θ, where θ is the angle between the Sun's direction and the total momentum of the final-state stable particles obtained from the procedure described above, is shown in Fig. 1. In spite of the smearing from nuclear effects and detector resolution, Fig. 1 quantitatively demonstrates that the angular correlation in the BDM signal events persists, and can be exploited to distinguish from the uniformly distributed background. Selection criteria on cos θ are optimized to different signal samples, and efficiency ( Ar ) and expected number of background events (b) in each selection are used to obtain the sensitivity of the BDM search. Details of this analysis are described in the Supplementary Material. Discussion.
The projected sensitivity for 10 years of livetime with a DUNE-like detector with 40 kilotons of LAr The expected reach is compared to current constraints from BDM (χ) sensitivity with Super-Kamiokande and spin-dependent direct detection searches for non-relativistic DM (ψ) at PICO-60L [5]. The leading constraints from direct detection of neutron scattering by PandaX [7] are subdominant for the models we consider.
is shown in Fig. 2, in which the gauge coupling constant in the benchmark model, g 4 Z , is excluded at two standard deviations over the range of parameters considered. Under the assumption that the χ relic abundance is negligible and undetectable by direct detection experiments, we compare the sensitivity of this benchmark model with the current constraint at Super-Kamiokande by reinterpreting their atmospheric neutrino measurement [64], as detailed in the Supplementary Material, and spin-dependent direct detection searches for ψ [5]. Since Super-Kamiokande does not have sensitivity to BDM at γ = 1.1 due to its high threshold for protons, it is absent in the first panel of Fig. 2. It is worth noting that in a supporting study we also find that the fermionic BDM shows kinematic characteristics similar to the scalar BDM, and similar sensitivity can be achieved, while its parameter space is more constrained by direct detection experiments [44].
We demonstrate that underground, massive LArTPC detectors can have unique, complementary capability of searching for BDM, taking into consideration the realistic nuclear effects, detector resolution, and background for the first time. Notably, we show that the sensitivity of this search technique is not compromised by the nuclear effects. The method and the dedicated GENIE package we develop for this study offer the means to characterize the BDM interactions in hadronic channels event by event, and is straightforward to adapt to different models of nuclear effects and atmospheric neutrino fluxes. Our simulations offer the most accurate description of the signal and backgrounds for hadronic interactions to date.
BDM models are well-motivated and are gaining attention. The framework of BDM [78][79][80][81][82] has been explored for interpreting the new observation from the Xenon1T experiment [83]. In that same context, this analysis pipeline will help narrow down the parameter space possible if a complementary observation occurs in LArTPC detectors; it provides kinematic information to distinguish the signal across different classes of BDM models, cross match with potential leptonic interactions also detectable with DUNE, and from neutrino detection of models of ψψ → νν. Additional studies that combine potential DUNE results with those from direct detection or other experiments can help to further determine the properties of a BDM model.
We presented the first dedicated study of BDM search via hadronic interactions in underground, massive LArTPCs, paving the avenue for future, sophisticated analyses. The work can also be extended to investigate the requirements of detector specifications and reconstruction criteria for BDM and similar astroparticle searches.

Acknowldgments.
We thank Costas Andreopoulos, Robert Hatcher, and Marco Roda for support related to GE-NIE. We are grateful to Jonathan Asaadi, Mark Convery and Hirohisa Tanaka, for all the discussion about the features of different detectors, and Aaron Higuera, for the conversation regarding the rate of atmospheric neutrinos. We also thank Jesse Thaler and Kaustubh Agashe for discussion. JB is supported in part by PITT  In this Supplementary Material we present certain details of our analysis that may be of interest to the reader, but are not essential to understanding our work. In Sec. I, we describe the physics entering the GENIE BDM event generation module used in our analysis. Sec. II discusses the effects of nuclear interactions on the observable particles. We describe the background event generation procedures specific to our analysis in Sec. III. Our detector simulation procedure is outlined in Sec. IV. Finally, in Sec. V we describe our analysis strategy in detail.

I. BOOSTED DARK MATTER EVENT GENERATION (PARTICLE LEVEL)
The liquid argon target is comprised, fundamentally, of electrons, quarks, and gluons. Existing studies on BDM scattering so far have focused on BDM-electron or BDM-nucleon (or BDM-H) scattering. A detailed study on BDM-nucleus scattering is lacking and involves complex nuclear effects. In this section, we translate quark-level interactions into cross-sections and event generation for interactions of BDM with the nuclei in a target. These interactions are implemented in the GENIE Monte Carlo event generator and are now a part of GENIE version 3. The full details of these interactions and the software package are presented in Ref. [67].
There are several regimes for this interaction depending on the kinematically allowed momentum transfers for the interaction. As is standard with neutrino-nucleus scattering, we parameterize the momentum transfer with Q 2 = −q 2 = −(k − k) 2 , where k and k are the four-momenta of the incoming and outgoing BDM particle respectively. For an elastic scattering on a nucleon at rest, one can relate this Q 2 to the outgoing nucleon kinetic energy by Q 2 = 2M N · E k,N , where M N is the nucleon mass and E k,N is the outgoing nucleon kinetic energy.
At low momentum transfers (Q 2 (100 MeV) 2 ), only coherent scattering off the nucleus is possible. Since isotopes of argon with an odd number of neutrons are very rare and we are considering models where spin-dependent interactions dominate, this process is highly suppressed in argon and we neglect it entirely.
For (100 MeV) 2 < ∼ Q 2 < ∼ (1 GeV) 2 , the only significant process is dark matter elastic scattering off of nucleons χ + N → χ + N , where N refers to a nucleon, that is a proton p or a neutron n. We will refer to this process simply as elastic scattering, as is conventional in studying neutrino scattering. In this regime, nuclear effects such as Fermi motion and Pauli blocking are relevant. Furthermore, at higher momentum transfers, the nucleon form factor becomes an important effect. All of these effects are described in detail below, in Sec. I A.
For (800 MeV) 2 < ∼ Q 2 < ∼ (1.8 GeV) 2 , inelastic scattering begins to become an important process. At these threshold momentum transfers, inelastic scattering is dominated by resonant production of excited baryons N * and ∆, χ + N → χ + N * and χ + N → χ + ∆. This process, called resonant scattering, is rather complicated to describe and suffers from large modeling uncertainties. In the present analysis, we omit these processes, rendering the limit projections we derive somewhat more conservative. Description and modeling of these interactions will be performed in future work.
For Q 2 > ∼ (2 GeV) 2 , deep inelastic scattering off partons in the nucleons, χ + q → χ + q, becomes an increasingly good description. Diagrams illustrating a typical interaction for each of these processes are shown in Fig. 3. We now provide a detailed description of the elastic and deep inelastic scattering cross-sections and other relevant physics for each process.
Event generation for scattering in GENIE proceeds via a series of modules that implement the relevant nuclear and particle physics. Most of these implement nuclear physics effects, such as Fermi motion, Pauli Blocking and final state nuclear interactions. These remain unchanged from their neutrino scattering implementation for BDM scattering. We therefore focus below on the determination of the differential scattering cross-section as well as the BDM kinematics, which are the points at which BDM scattering differs from neutrino scattering. We work here exclusively in the nucleon rest frame, which is not the same as the lab frame because of nucleon Fermi motion, but is reached by a trivial boost of the BDM-nucleon system.

A. Elastic Scattering
The differential cross-section for elastic neutrino scattering in GENIE follows the calculation of Ahrens et. al. [84], though the formalism has been developed elsewhere in the literature and is standard. As discussed in the Letter, we focus on the case where BDM interacts via a spin 1 boson that has axial couplings to the quarks. The amplitude for elastic scattering then depends on hadronic matrix elements of the form: where F A and F P are the axial and pseudo-scalar form factors respectively, q = k − k is the momentum transfer four-vector, and m N is the nucleon mass. For scalar BDM, it is straightforward to show that the term involving F P vanishes in the amplitude. The differential cross-section in the single kinematic variable Q 2 can be written as following the construction of Ref. [84]. The parameters are given by with Here, Q χ is the Z charge of the scalar BDM and s, u are Mandelstam variables, and E χ is the energy of the incident BDM. The form factor F A is assumed to have a dipole form, where M A is a parameter that needs to be fit to data. The default value for this parameter in GENIE, which we keep, is 0.99 GeV.
The normalization of this form factor is given, in general, by a combination of the spin form factors of the nucleon. Assuming isospin symmetry, the form factors for the proton and neutron are with the quark axial charges Q f . The spin form factors need to be either extracted from data or calculated on the lattice. We take them to be [70] ∆u = 0.84, ∆d = −0.43, ∆s = −0.09.
Note that the range of momentum transfers is given by The phase space for deep inelastic scattering (DIS) is described by two, rather than one, variables, in addition to the complex hadronic phase space determined by the hadronization procedure. One intuitive way of breaking down the phase space here is in terms of the momentum transfer Q 2 and the total invariant mass of the final state hadronic system W .
While the variables Q 2 and W are physically intuitive, it is simpler to describe the cross-section in terms of variables x and y, where x is the usual Bjorken variable and y is the fractional energy loss of the incoming DM particle, where E χ and E χ are the energy of incoming and outgoing DM, respectively. These variables can be written in Lorentz invariant form, related to Q 2 and W 2 as Note that these variables range in a subset of 0 < x, y < 1 that can be solved for numerically. To proceed and calculate the cross-section, we follow closely the notation of Ref. [85]. We define the hadronic tensor as the initial spin averaged, final state summed squared hadronic matrix element at fixed Q 2 and W 2 , summed and integrated overall all possible final states. By Lorentz invariance, the hadronic tensor has the form where p is the four-momentum of the initial nucleon. The F i are structure functions that are related to the quark PDFs below. For scalar DM scattering, we find The structure functions F i here are given in terms of the quark PDFs by the following relations by where f f are the parton distribution functions for quark flavor f , combined with the Callan-Gross relation and the Albright-Jarlskog relations The parton distributions used in GENIE are a patched version of the GRV98lo PDFs [86].
Once the x and y of a DIS event are selected, the hadronic final state phase space is then populated using one of two hadronization models. At low energies, an empirical Koba-Nielson-Olesen (KNO) model is used in the neutrino. Absent empirical observation of BDM, we must make an assumption of the empirical behavior to implement this model for DM. We assume that BDM scattering behaves like neutrino scattering in the KNO model. At high energies, PYTHIA is used to hadronize the final state hadronic system. This procedure remains unchanged for BDM scattering.
For details of DIS interactions for fermionic BDM, see [67].

C. Resonant scattering
Resonant scattering via excited baryon states is implemented for neutrino scattering in GENIE, but the implementation of their models for BDM is challenging to validate. This process is not studied in the present analysis, though it can only increase the sensitivity to BDM.

II. IMPACT FROM NUCLEAR EFFECTS
In the recent years, interest in the interactions of hadrons produced within the nucleus on their way out of the nuclear remnant ("final state interactions") has surged within the community owing to their significant impact on precision measurements of neutrino oscillations and search for nucleon decays, especially with detectors based on large nuclei like oxygen and argon.
In this analysis, the default hA final state interaction (FSI) model in GENIE [66,68,69] is used to model the nuclear FSI, but it is straightforward to switch to different models in the future data analysis with the tool we developed for this study. This model uses empirically determined total cross sections for various processes that hadrons propagating through the nuclear remnant can undergo, such as pion absorption, elastic and inelastic scattering, and charge exchange. The cross sections are extrapolated to high energies where data is unavailable. The focus of this model was on iron for the MINOS experiment. Alternate models currently include the hN model, which implements a more complex intranuclear cascade designed for situations with multiple scattering. On the other hand, it currently does not include important medium corrections [87]. More recent iterations of both models have been developed as well.
In addition to FSI, there is some modification of the kinematics due to nucleon motion within the nucleus and Pauli blocking. GENIE models these phenomena with a Fermi gas. Figure 4 illustrates the impact of the nuclear effects on the distribution of cos θ, where θ is the angle between the total momentum of the final state visible particles (i.e. excluding neutrons and the outgoing BDM) and the incident BDM. The kinematic feature at cos θ ≈ 0.25, which originates from the elastic component of the scattering, gets smeared out dominantly by the effects of nucleon Fermi motion, which misaligns the detector and nucleon rest frames.

III. BACKGROUND
As stated in the Letter, we consider the neutral current interactions of neutrinos produced in the atmosphere as the main background. The absolute rate of interactions from atmospheric neutrinos is calculated integrating the neutrino cross sections with argon with the atmospheric neutrino fluxes using Bartol model [75] at the MINOS location in Soudan, as in [44] p. . The simulation of the background processes is performed using the LArSoft toolkit [88] interfaced with GENIE, including generation of ν µ ,ν µ , ν e andν e distributed according to the reference Bartol flux. Additional corrections to account for a larger energy spectrum and for neutrino flavor oscillations are described below.
As with the signal, each generated neutrino interaction is assigned a direction toward the Sun randomly extracted from the unbiased distribution of the Sun position with respect to the DUNE far detector. This direction is solely used to estimate the angle cos θ used as an observable in this study.
A. Solar magnetic activity Solar magnetic activity affects cosmic ray deflection and as a consequence the rate and spectrum of atmospheric neutrinos. The activity oscillates between minimal and maximal with a period of about 11 years. Atmospheric neutrino fluxes include the effects of this activity, and are provided separately in the minimum and maximum activity scenarios. In our study we mix for each process two samples simulated with the two scenarios. Since the period of this activity is close to the 10 years of duration of data taking we consider in this study, we assume one full cycle and therefore we mix the two samples with equal weight.

B. Extension of background estimation to higher neutrino energy
Our atmospheric neutrino background estimation includes only neutrino energies between 100 MeV and 10 GeV, being based on Bartol flux at Soudan [89]. With our choice of parameters, BDM interactions can cover a larger energy range and we need to extend the background estimation to cover that range; because the neutrino flux rapidly decreases with the neutrino energy E ν , roughly as E −2 ν , we elect to extend the coverage only up to 100 GeV. To do so, we use the Honda flux (at Homestake), which extends up to 10 TeV, scaling it so that the flux integrated in the energy range 1 to 10 GeV matches DUNE background estimation. The approximation implied by this procedure is that the two atmospheric neutrino models, Bartol and Honda, scale with energy in the same way. We estimate this approximation to carry an error of about 20%. The choice not to employ Bartol fluxes for this extension is purely technical, due to a temporary issue in the LArSoft software. Likewise, the choice to use samples with narrow energy ranges is technical. Due to the steep decrease of flux with energy, the generation of a single sample with large energy range uses computational resources very inefficiently.
High energy neutrino events inherently present kinematics different from lower energy ones. In our simple analysis events are selected according to a single quantity: the angle (cos θ) between the direction of the reconstructed particles and direction of the Sun. These two directions are uncorrelated for the atmospheric neutrino background and cos θ is mostly independent from the neutrino energy. We confirm that the Honda high energy neutrino NC interaction sample shows the same distribution in cos θ as the Bartol atmospheric neutrino sample that constitutes our reference background (see Fig. 1) and that both cos θ distributions are consistently uniform. Because of this, we simply retain the reference background sample in this analysis, scaling its size up by a factor to account for the high energy contribution. The small size of the resulting correction, 3.8%, suggests that interactions with even higher energy neutrinos above 100 GeV will contribute negligibly to this background.

C. Tau neutrino background
Under the assumption of being able to identify and discard background events where charged-current interactions produce electrons or muons, our background is constituted mainly of atmospheric neutrino interactions via neutral current. An exception is a charged-current interaction where a τ lepton is produced that decays into a neutrino and hadrons. This happens with a branching fraction close to 60%. While neutrinos of τ flavor are rarely produced in the interaction of cosmic rays with the atmosphere, it is still possible for a muon neutrino to transform ("oscillate") into ν τ . The probability of this transition for a muon neutrino of energy E ν is described by the formula where the parameters θ 13 , θ 23 and ∆m 2 31 have been measured [90]. An accurate computation of the rate is complicated by the dependency on L, the distance from the point in the atmosphere where the neutrino is produced to the point in the detector where it interacts. This distance can be as short as a few kilometers for neutrinos produced right above the detector, to more than ten thousand kilometers for the ones produced at the opposite side of the Earth; this is compared to the factor ∆m 2 31 /4 ≈ 3 MeV/km. We simplify the problem by the very conservative approximation of oscillation probability being maximal independently from E ν , by setting the last of the three factors of the expression above to 1, yielding P (ν µ → ν τ ) ≈ 95%, with the understanding that this represents a significant overestimation of this component of the background.
Oscillation does not have observable consequences on neutral-current interaction backgrounds. We generate the tau neutrino sample using the same Honda flux as for muon neutrinos. To ensure that the size of the ν τ samples is consistent with the other background samples, we impose the same rate of interaction via neutral current for ν µ and ν τ of the same energy, by properly scaling the tau neutrino interaction rate. The rate of interaction of ν τ via charged current, the one relevant for this part of the background, is scaled with the same factor, but it remains much smaller than for ν µ at low energy, being suppressed by the larger mass of τ lepton. The charged-current ν τ interactions are mostly suppressed for E ν < 10 GeV, whereas the charged-current ν µ interactions with the same energy range, E ν < 10 GeV, constitutes 90% of those in the full energy range we consider, 100 MeV < E ν < 100 GeV. For this reason, the ν τ charged-current background is much smaller than the atmospheric neutrino neutral-current background, 2.8% in our estimation. As for the contribution to the background from the high energy extension, the reference background sample, i.e. the atmospheric neutrino sample based on Bartol flux and with no oscillation, is still used for the analysis, and the effect of oscillation into ν τ is included as a 2.8% correction factor on the total background rate.

IV. DETECTOR SIMULATION
For each event, we use GEANT4 [91][92][93] to simulate the propagation of the final-state SM particles in liquid argon until any short-lived particles have decayed. The four-momentum of the stable particles, protons, neutrons, charged pions, muons, electrons, and photons, are convolved based on the parameters characterizing the detector response. The convolution of the four-momenta accounts for the energy and angular resolution of the detector. Only the particles with their convolved energy greater than the detector threshold are taken into account in the subsequent steps of the analysis. The baseline scenario deployed in this analysis consists of the detector response and threshold reported in the DUNE Conceptual Design Report (CDR) [77], as listed in Table I, and, of no neutron detection. The resulting sensitivity on BDM search is presented in the Letter.

Particle type Detection Threshold (KE) Energy Resolution
Angular Resolution To evaluate the impact from the detector response and threshold, as well as the capability of reconstructing neutrons, an alternative set of energy resolution is deployed. This set of tabulated energy resolution was obtained and studied by the authors of Ref. [94]. In addition, we study the cases where 90% of neutrons can be detected and reconstructed. Further, we lower the detection threshold to 20 MeV in kinetic energy (KE) for protons, neutrons, and to 30 MeV in KE for charged pions, labeled as the "optimistic" scenario for detection thresholds. All the scenarios are outlined in Table II.
We obtain similar sensitivity on BDM search from all the scenarios being tested. This is owing to the fact that we deploy a simple analysis approach, as depicted in Section V, and do not utilize plenty of information to which better energy resolution, lower detection threshold, or capability of neutron reconstruction is relevant.   , and expected number of signal (s) and background events (b) after cuts for our benchmark models (mass mχ and boost γ) assuming an exposure of 40 kiloton and 10 years.

V. ANALYSIS
The BDM signal events are expected to have final-state particles roughly aligned with the incoming BDM particle, which we take to be coming from the Sun. We use this feature to select events with enhanced the signal-to-background ratio.

A. Baseline Analysis
We develop selection criteria based on θ, which, as defined in the Letter, is the angle between the total momentum of the final-state stable SM particles and the incident BDM (aligned with the Sun). The detector response and thresholds are taken into account, as described in Sec. IV. The single variate selection is optimized to the minimal signal strength, s , for which we could obtain a sensitivity to BDM signal at 5 standard deviations, The factor Ar represents effectively the product of the acceptance and efficiency of the signal selection, while the expected number of selected BDM events, s, can be written as s = Ar s . We evaluate Ar and the number of the background events b respectively from the BDM and atmospheric neutrino MC samples, and the selection criterion on cos θ is individually optimized to each benchmark BDM signal sample, as tabulated in Table III. Note that Eq. 21 is used for optimizing the selection criteria, but not for extracting the sensitivity to the BDM signal.

B. Alternative Selection Oriented to Moderate Boost Signals
Owing to the kinematics of the elastic scattering, the cos θ distributions in the moderate boost signal samples (e.g. γ = 1.1) are more widely spread, and, as a consequence, the single variate selection based on cos θ is less efficient, resulting in a smaller Ar and a greater b in Table III. To improve the signal-to-background ratio in these cases, the kinematic correlation between cos θ and P is studied, where P denotes the value of the total three-momentum of the final-state SM particles. In the limit that the nucleon is at rest when it is struck by the incident BDM, these variables are perfectly correlated for a given model at fixed invariant mass for the final-state hadronic system. This correlation does not hold for the background and should allow for further separation of signal and background. A few cut-based analyses using the two variables, cos θ and P , are explored. With the simple statistic estimate used in this analysis, the sensitivity is comparable to the baseline analysis; however, we expect more sophisticated analyses, with better understanding on nuclear effects and detector response, to significantly improve the sensitivity of the BDM search.

C. Impacts from detector response and threshold
We study the impacts from different scenarios of detector response and threshold by performing the baseline analysis with the convolved four-momentum from all the scenarios listed in Table II. In addition, we compare the results combining different detector response and analysis strategies (baseline analysis versus alternative analysis). Similar to the conclusion obtained from Sec. V B, to significantly improve the sensitivity of BDM search requires better understanding on the BDM signal and atmospheric neutrino background, including nuclear effects and the flux of atmospheric neutrinos, as well as more sophisticated analyses.

D. Statistical Method
We obtain the projected sensitivity for 10 years of livetime with a DUNE-like detector with 40 kiloton of LAr. Since we expect a large number of signal events for the parameter space at the boundary of the discovery reach, the expected significance is evaluated with a large statistics estimate [90], where s and b are the numbers of expected signal and background events respectively. We find that our LArTPC reference detector is sensitive to g 4 Z = (1.54 − 22.0) × 10 −7 at two standard deviations over the range of parameters considered, as shown in Fig. 2.

E. Super-Kamiokande Data Comparison
A reinterpretation of the NC elastic ν + p → ν + p measurement from the atmospheric neutrino events collected in Super-Kamiokande is performed for comparison with our LArTPC analysis. The BDM events scattered on hydrogen and oxygen atoms in Super-Kamiokande are simulated by the same BDM module in GENIE, and, accounting for the Cherenkov threshold and the Cherenkov cone selection and efficiency in [64], the events containing a single proton with momentum between 1.07 and 2.62 GeV are selected and scaled. The sensitivity of the BDM signals in Super-Kamiokande is thereby evaluated based on the simulated BDM events and the atmospheric neutrino data in Tables I and II in [64].