Optimizing Energetic Light Dark Matter Searches in Dark Matter and Neutrino Experiments

Neutrino and dark matter experiments with large-volume ($\gtrsim 1$ ton) detectors can provide excellent sensitivity to signals induced by energetic light dark matter coming from the present universe. Taking boosted dark matter as a concrete example of energetic light dark matter, we scrutinize two representative search channels, electron scattering and proton scattering including deep inelastic scattering processes, in the context of elastic and inelastic boosted dark matter, in a completely detector-independent manner. In this work, a dark gauge boson is adopted as the particle to mediate the interactions between the Standard Model particles and boosted dark matter. We find that the signal sensitivity of the two channels highly depends on the (mass-)parameter region to probe, so search strategies and channels should be designed sensibly especially at the earlier stage of experiments. In particular, the contribution from the boosted-dark-matter-initiated deep inelastic scattering can be subleading (important) compared to the quasi-elastic proton scattering, if the mass of the mediator is below (above) $\mathcal{O}$(GeV). We demonstrate how to practically perform searches and relevant analyses, employing example detectors such as DarkSide-20k, DUNE, Hyper-Kamiokande, and DeepCore, with their respective detector specifications taken into consideration. For other potential detectors we provide a summary table, collecting relevant information, from which similar studies can be fulfilled readily.


Introduction
The dark matter puzzle is a clear motivation for physics beyond the Standard Model (SM). While the evidence for the existence of dark matter is all rooted in its gravitational interaction with ordinary matter, searches via its non-gravitational couplings are actively underway. Of those trials, the strategy of dark matter direct detection has been playing a role of the major driver in searching for relevant signatures. Most of direct detection experiments are designed to observe a recoil of target material (henceforth called primary signature) which is induced by the elastic scattering of non-relativistic dark matter (see e.g. Ref. [1] for a modern review on dark matter direct detection experiments).
A variation in this search scheme is to look for inelastic scattering signals -which was originally motivated by the DAMA annual modulation signature [2] -, imagining the process that a dark matter particle scatters off to an excited state along with a target recoil whose energy spectrum differs from that in the elastic scattering mentioned above [3].
By construction, the excited state may be de-excited back to the dark matter state as previously expected in [4], potentially leaving visible signals at the detector (henceforth called secondary signature): for example, X-ray photon in neutrinoless double beta decay experiments [5]. One may attempt to observe both primary and secondary signatures, but it is usually challenging for both of them to overcome the relevant detector threshold simultaneously due to inadequate dark matter kinetic energy. Indeed, the rich structure of inelastic dark matter models have endowed themselves with the potential to explain a diverse range of astrophysical phenomena [4,[6][7][8] and have inspired novel LHC search strategies [9][10][11][12][13].
A myriad of experimental efforts to observe dark matter signals have been devoted under the search schemes discussed above. However, none of them have recorded solid darkmatter-induced signatures yet, so they merely sets stringent bounds on parameter space of associated dark matter models. The null observation motivates alternative approaches. One possible direction to pursue is to look for similar experimental signatures invoked by dark matter relics of different mass scales [14,15]. Since most of the existing direct search experiments aim at weak-scale dark matter, hence the associated detectors are designed accordingly, new detector material and/or technology are often demanded in order to perform relevant experiments [14,15]. Alternatively, one may search for any scattering signatures of relativistically incoming dark matter, usually having in mind the mass scale of standard thermal dark matter.
A straightforward production mechanism is to obtain relativistic dark matter at particle accelerators, where a certain fraction of initial-state beam energy is transferred to the dark matter. The elusive nature of typical dark matter often requires highly intensified particle beam essentially to increase signal statistics, e.g., fixed target experiments [14][15][16]. The larger fraction of literature studied elastic scattering of relativistically produced dark matter, e.g., Refs. [17][18][19][20][21], but the energetic nature of such dark matter essentially allows decent cross section for its "up"-scattering to an excited (or equivalently heavier unstable) state, under the framework of inelastic dark matter. Reference [22] pointed out the potential of detecting both the primary recoil induced by such relativistic dark matter and visible decay product(s) of the excited state and showed that it allows relevant signal searches to suffer from significantly less background contamination, hence inducing recent development in phenomenological investigations [23,24].
Another class of mechanisms invokes production of energetic (light) dark matter in the present universe. We emphasize that our study here is straightforwardly applicable to several physics scenarios, models, or frameworks involving relativistically produced dark matter, but it is instructive to develop our argument in the context of a concrete example. Such an example is the scenario of boosted dark matter (BDM) [25]. In BDM models, it is assumed that a dark sector containing two dark matter species with a hierarchical mass spectrum [26]. Separate symmetries are usually employed to stabilize the two dark matter species, e.g., Z 2 × Z 2 or U(1) × U(1) . The overall dark matter relic is set by the so-called "assisted freeze-out" mechanism [26]. Suppose that the heavier and the lighter species are denoted by χ 0 and χ 1 , respectively. Typical models hypothesize that χ 0 does not directly interact with SM particles but pair-annihilates into a χ 1 pair which directly couples to SM particles, i.e., the χ 0 relic abundance is determined by "assistance" of χ 1 . As a result, χ 0 usually remains as the dominant relic component, whereas χ 1 becomes subdominant. Under this setup, it is hard to detect relic χ 0 at standard dark matter direct detection experiments due to its tiny coupling strength to SM particles, while it is again hard to detect relic χ 1 due to its small statistics in the current universe. However, the model setup allows for χ 1 production via pair-annihilation of χ 0 in the galactic halo. The produced light dark matter χ 1 acquires a significant Lorentz boost factor due to the mass gap between χ 0 and χ 1 , which opens a novel physics opportunity, the search for signatures of relativistic dark matter scattering. The signal detection prospect is deeply related to the total χ 1 flux F 1 [25]: where m 0 denotes the mass of χ 0 which consists of most of the dark matter relic in wellmotivated regions of parameter space. Here the reference value for σv 0→1 , the velocityaveraged annihilation cross section of χ 0 to χ 1 , corresponds to a correct dark matter thermal relic abundance of χ 0 , under the assumption that an s-channel process dominates the annihilation process. The above relation implies that for χ 0 of weak-scale mass (i.e., ∼ 100 GeV), the incoming flux of lighter dark matter χ 1 (near the earth) is as small as O(10 −8 cm −2 s −1 ). Thus, large-volume neutrino detectors such as Super-Kamiokande (SK), Hyper-Kamiokande (HK), Deep Underground Neutrino Experiment (DUNE), and IceCube Neutrino Observatory are preferred in the search for experimental signatures induced by χ 1 elastic scattering-off detector material [25,[27][28][29][30][31][32][33] or those by χ 1 inelastic scattering (which can be viewed as a combination of the original boosted dark matter scenario and the inelastic dark matter model) [34]. We call the former and the latter scenarios eBDM and iBDM, respectively, as shorthand throughout the rest of this paper. The SK Collaboration has conducted the search for high-energy electron recoil ( > ∼ 0.1 GeV) induced by χ 1 elastic scattering and reported the first results [35]. On the other hand, it was shown that sub-GeV/GeV-range m 0 can increase the χ 1 flux substantially, while keeping the resultant relic abundance consistent with the current measurement, so that ton-scale dark matter direct detection experiments such as Xenon1T and LUX-ZEPLIN can be sensitive enough to eBDM/iBDM signals [36][37][38]. Recently, the COSINE-100 Collaboration has searched for the electron-positron pair in coincidence with the primary electron signal by χ 1 inelastic scattering as a signature of an iBDM interaction and reported the first results of direct search for iBDM [39].
We remark that in terms of recoiling target particles, two channels are considered: an electron target and a proton target. The former is an elementary particle so that predicting a χ 1 scattering cross section with electrons in the detector material is rather straightforward. On the other hand, the latter is a composite object and is typically bound in nuclei, and as a result, the corresponding scattering cross section involves form factors. Moreover, since incident light dark matter χ 1 is relativistic, a deep inelastic scattering (DIS) process may arise. When it comes to iBDM, the proton target is often advantageous for χ 1 to "up"-scatter off to a heavier (unstable) dark-sector state [34,37]. Considering these factors altogether, therefore, it is of great importance to choose a better channel for a given parameter region to explore. This strategical approach is highly motivated, in particular, at the earlier stage of experiments. For example, several sub-kiloton-scale neutrino detectors in Short Baseline Neutrino Program (SBN) [40,41] and prototypical DUNE (ProtoDUNE) [42,43] are in operation or ready to take data within a few months, and physics opportunities in terms of iBDM [44] and eBDM [45] have been proposed recently.
In light of this situation, we perform a dedicated study to provide useful guidance for boosted dark matter searches in this paper which can be taken as a reference search for energetic light dark matter coming from the universe. In more detail, we first show that the DIS contribution out of an entire proton scattering cross section is negligible, as long as the mass of the particle mediating the interaction between energetic light dark matter χ 1 and SM particles is not much larger than the energy scale inducing a DIS process. Therefore, in many of the well-motivated scenarios, it is sufficient to consider only contributions by genuine proton scattering. We then compare the electron channel and the proton channel through their respective scattering cross sections for both eBDM and iBDM signals, including realistic factors such as detector energy threshold, cuts, and angular resolution at several benchmark detectors.
To deliver the main ideas efficiently, our paper is organized as follows. In Section 2, we briefly discuss the scenario of boosted dark matter and an example model to describe the interactions between lighter dark matter species χ 1 and SM particles, followed by listing up several benchmark detectors and their key characteristics. We then look into scattering cross sections of an incoming χ 1 with electron and proton targets in Section 3.1, putting a particular emphasis on the proton DIS scattering. Detailed comparison between electron and proton scattering channels for both eBDM and iBDM follows in Sections 3.2 and 3.3 at the theory level and the (semi-)detector level, respectively, while we scan over the mediator mass and the χ 1 mass for a given m 0 . Example phenomenology will be discussed and demonstrated in Section 4, and our concluding remarks will appear in Section 5. For the sake of reference, we provide a couple of appendices. In Appendix A, we provide our derivation for various scattering cross section formulas and some detailed description of our data analysis. A summary of key specifications of detectors other than the benchmark ones is presented in Appendix B.

Benchmark Models and Detectors
We begin with setting up benchmark models and detectors with which our detailed analysis will be demonstrated. The method can be straightforwardly extended to other models and detectors, but the case studies that we are performing in this paper will provide a baseline for other applications.

Dark matter models and experimental signatures
The simplified model under consideration is divided into two parts, the one delineating the production mechanism of boosted dark matter at the universe today and the one describing the interaction between SM particles and the boosted dark matter. Expected experimental signatures follow once the latter part is defined.
We note that there are several ways to create relativistic (or at least fast-moving) dark matter particles in the universe: for example, two-component dark matter scenario [25,26,46,47], models with a Z 3 symmetry which may induce semi-annihilation processes [48], models involving anti-baryon-numbered dark matter-induced nucleon decays inside the sun [49], scenarios with decaying super-heavy particles [29,30,50], or energetic cosmic-ray induced (semi-)relativistic dark matter scenarios [51][52][53]. One can also think of various places from which boosted dark matter dominantly comes. Examples include the galactic center (GC) [25], the sun [27,28], and dwarf galaxies [31]. Among those possibilities, we simply choose the two-component dark matter scenario in which the dominant flux of boosted dark matter comes from the GC, as our benchmark model. In this scenario, the boost factor of BDM is a free parameter, determined by the mass gap between the two dark matter species. Thus, we can change the boost factor freely and study/show the resulting phenomenological effect without any model restriction.
As briefly explained in the Introduction, one of the two dark matter species (usually the heavier χ 0 ) indirectly couples to SM-sector particles via the other dark matter species (usually the lighter χ 1 ) which directly communicates with SM particles. Assuming that dark matter relic abundance is determined thermally, we see that the χ 0 relic is set with the aid of χ 1 [26]. In more detail, χ 0 is in (indirect) contact with the thermal bath via χ 0 χ 0 → χ 1 χ 1 followed by sufficiently large χ 1 χ 1 → SM SM, in annihilation models described by an effective operator, e.g., where Λ parameterizes some high-scale physics. The thermally averaged annihilation cross section of χ 0 χ 0 → χ 1 χ 1 at the present universe is fixed to be ∼ 10 −26 cm 3 s −1 yielding the observed value of dark matter relic abundance, with the assumption that an s-wave annihilation dominates χ 0 χ 0 → χ 1 χ 1 . The calculation procedure to find the flux of χ 1 is quite similar to that for the photon flux via dark matter pair-annihilation, so we write it as follows: where ρ describes the χ 0 density distribution in terms of the line of sight (l.o.s.) s and solid angle Ω with θ being the angle between the direction of l.o.s. and the axis connecting the GC and the earth. We essentially reproduce the result in Eq. (1.1) in the second line of the above formula, assuming that the χ 0 is distinguishable from its anti-particle, saȳ χ 0 , and that the dark matter halo is distributed according to the Navarro-Frenk-White profile [54,55]. In the other case in which χ 0 andχ 0 are indistinguishable, one can simply drop the prefactor 1/2 in the formula. It is noteworthy that the flux formula (2.2) is independent of the type of interactions between the dark sector and SM sector as we have not assumed any particular type. As for the interactions between the dark sector and the SM sector, we adopt a dark gauge boson scenario as our reference model. For a general analysis, we take the framework of inelastic boosted dark matter [34,37,44] which includes an additional (unstable) dark sector state χ 2 heavier than χ 1 . We allow both diagonal interaction of χ 1 with target (i.e., eBDM) and off-diagonal interaction (i.e., iBDM). Minimal ingredients forming the relevant sector of our reference model are χ 1 , χ 2 , and a hidden massive gauge boson X µ . The interaction Lagrangian includes the following operators: where the first term describes the kinetic mixing between U(1) EM and U(1) X [56][57][58][59][60][61][62][63][64], that is, field strength tensors for the ordinary photon and a hidden gauge boson F µν and X µν are mixed by parameter . The diagonal and off-diagonal gauge interactions are parameterized by the couplings g 11 and g 12 . The example realization of such interactions (in particular, flavor-changing currents) in a model construction was discussed in Refs. [3,37]. 1 Of course, one can alternatively consider Higgs portal type [70][71][72][73][74] or dipole type [37] interactions in replacement of X µ and other types of dark sector particles such as scalars and vectors. However, the analysis method in this paper still goes through even in those alternative scenarios. Given the interactions in (2.3), three classes of experimental signatures can arise as shown in Figure 1.
• The topmost diagram depicts nothing but an ordinary elastic scattering process of χ 1 involving an electron or proton recoil.
• The second diagram sketches a scattering signature where the secondary process happens "promptly" so that a fermion pair ff from the decay of the dark gauge boson comes out of the primary vertex, within detector vertex resolution.
• On the other hand, in the last diagram, either X or χ 2 is relatively long-lived so that a pair of time-correlated primary and secondary vertices appear displaced.
Obviously, the first one is relevant to eBDM, whereas the other two are iBDM-induced.
Here we denote the decay products of X by generic fermion f , but in the detection level some clarification is needed. In particular, if f is either τ lepton or quark, its appearance can be traced back by its decay products or hadronized objects. We enumerate several well-motivated truth-level final states and briefly discuss their experimental features.
• e + e − : Electron (positron) is one of the easily reconstructible objects in various types of detectors including Cherenkov and Liquid Argon Time Projection Chamber (LArTPC) ones. It actively emits electromagnetic radiation while traveling through a detector medium. Such radiation manifests as a set of Cherenkov light in Cherenkov detectors and a series of showering in LArTPC detectors.
• µ + µ − : Muon is the cleanest object in typical detectors. Due to its massiveness relative to electron, electromagnetic radiation is less vigorous. So, muons have Cherenkov radiation rings with sharp outer edges unlike electrons having rings with blurred outer edges. In LArTPC detectors, it does not leave much showering, but ionizes atoms in the vicinity of its trajectory, creating a clean track. Furthermore, if it stops inside a detector, it decays and emits a Michel electron resulting in a "kink" at the end of the track. 2 • π + π − : This channel opens up once the mass difference between χ 2 and χ 1 is greater than twice the mass of charged pion. A pion behaves similar to a muon as both have the same electric charge and similar mass values. Since a pion undergoes nuclear interactions, it typically travels a shorter distance than a muon for a given energy. Nevertheless, distinguishing a pion from a muon (or vice versa) is highly nontrivial, so that some non-conventional techniques may be desired. For example, Ref. [75] takes convolutional neutral networks to differentiate charged pions from the others, in particular, muons, in a LArTPC detector, and reports 70 − 75% pion tagging efficiency (modulo ∼ 20% muon contamination).
As the mass gap between χ 2 and χ 1 increases, more variety of modes (e.g., π + π 0 π − , π + π − π + π − , etc.) become available. In other words, they often accompany multiples of the above bulleted particles, increasing the complexity of particle identification (PID). Before closing this subsection, we shall make a few comments for our study with iBDM-initiated processes. First, while the analysis method is completely applicable to any SM charged leptons and quark-induced final states, we take electron and positron for simplicity. Kinematically, they just require a small mass difference between χ 1 and χ 2 ( > ∼ 1 MeV); however, the minimum mass gap to open a dimuon channel is 211 MeV. A large boost factor of χ 1 , γ 1 is usually demanded to access a χ 2 with such a mass gap [34]. Therefore, lighter χ 2 (if exist) are preferred for a given γ 1 , or a too heavy χ 2 is inaccessible if γ 1 is not large enough. In this sense, the e + e − pair is highly motivated over the other signatures. Second, all three final state particles may be collimated, if an iBDM process is initiated by a χ 1 with a large boost factor. Angular resolution of the detector becomes crucial to separate the three particles. Depending on detectors, available are some features to identify or tag merged objects; for example, dE/dx in LArTPC detectors [76]. In our analysis, we consider the issue of angular separation when reporting the results for iBDM signals. Third, given the experimental signatures shown in Figure 1, potential backgrounds should be carefully identified and assessed especially for the sensitivity calculation. While iBDM suffers far less from background contamination, eBDM actually does because of its simple signature. Our main focus is comparisons among different signal channels. So, we for the moment pretend to be safe from background issues while referring to e.g., [35,44,76] for more systematic discussions. Finally, even eBDM models would give rise to signatures similar to the iBDM ones if a secondary process is accompanied via dark gauge gauge boson radiation-off of initial/final-state χ 1 (i.e., dark-strahlung) [77]. iBDM search strategies are essentially relevant to the dark-strahlung channel, but the interpretations of experimental results need to be conducted under a proper model hypothesis.
• DarkSide-20k [78]: As a unified dark matter direct detection program of the four LArbased projects (ArDM, DarkSide-50, DEAP-3600, and MiniCLEAN), the DarkSide-20k experiment has been approved due to the successful experience in operating the DarkSide-50 detector. In the framework of the DarkSide-20k experiment, a dualphase LArTPC with an active (fiducial) mass of 23 t (20 t) will be deployed at Laboratori Nazionali del Gran Sasso [a depth of ∼ 3800 meter water equivalent (m.w.e.)] in Italy, where the rate of cosmic rays is reduced to ∼ 1.1 m −2 hr −1 , in 2021. The LArTPC, i.e., the active (fiducial) volume, will be an octagonal shape with a height of 2.39 m (2.27 m) and a distance between parallel walls of 2.9 m (2.78 m). 3 • DUNE [79][80][81][82][83]: The far detector of DUNE will consist of four LArTPC modules and be located about 1500 m (≈ 4300 m.w.e.) underground at the Sanford Underground Research Facility in South Dakota, USA, where the expected rate of cosmic rays is ∼ 0.6 m −2 hr −1 . The first (second) module is planned to be ready for operations in 2026 (2027). Two LArTPC technologies, single-phase (SP) and dual-phase (DP), are planned. Each module will consist of a cryostat with internal dimensions 15.1 m (width) × 14.0 m (height) × 62.0 m (length) which will contain a total (fiducial) LAr mass of about 17.5 kt (at least 10 kt). These LArTPC detectors will have excellent angular resolution (θ res ), good PID capability, and relatively low energy threshold (E th ), e.g., θ res ∼ 1 • and E th ∼ 30 MeV for an electron.
• HK [84][85][86]: HK is a next generation underground water Cherenkov experiment based on the very successful operation of SK and will serve as a far detector of T2HK/T2HKK, a long baseline neutrino experiment for the upgraded J-PARC beam. It will consist of two detectors, each of which will have a cylindrical water tank with 60 m (51.8 m) in height and 74 m (67.8 m) in diameter and holds a total (fiducial) water mass of 258 kt (187 kt). The first detector will be hosted at Tochibora mine (a depth of 650 m ≈ 1750 m.w.e.) near the current SK site in Japan with the expected cosmic-ray rate of ∼ 27 m −2 hr −1 . The second one is currently considered to be constructed under the Mt. Bisul or Mt. Bohyun with about 1000 m overburden (≈ 2700 m.w.e.) in Korea with the reduced cosmic-ray rate of ∼ 5.7 m −2 hr −1 . The operation of the first detector in Japan will be ready in 2027. For electrons, HK will be able to have quite low energy threshold of ∼ 5 MeV, but much higher energy threshold is required to achieve good angular resolution, e.g., E th = 100 MeV for θ res ∼ 3 • .
• DeepCore [87,88]: DeepCore is a subarray of IceCube, an ice Cherenkov experiment, that has an approximately five times higher detector module density than that of the original IceCube array. It has been fully installed between 2100 and 2450 m below the surface of the icecap at the South Pole, where the expected cosmic-ray rate is ∼ 21 m −2 hr −1 , and taking physics data since May 2010. Due to a denser module array, DeepCore can lower the energy threshold by over an order of magnitude (E th ≈ 10 GeV) than that of IceCube. The effective target mass of DeepCore varies from ∼ 5 Mt to ∼ 30 Mt as the signal energy does from ∼ 10 GeV to ∼ 100 GeV since it has no clear boundary of target material. In this study, we take a cylindrical lump of ice with 350 m in height and 70 m in radius holding a total mass of 5 Mt as a conservative effective target mass with E th ≈ 10 GeV. PID is only good for the muon, and θ res ∼ 1 • for a muon-track event while θ res 10 • for a shower event.
A summary of key characteristics of various relevant detectors including benchmark ones are tabulated in Appendix B for convenience of reference.

Signal Cross Sections
The cross section formulas of the primary scattering of χ 1 with a target T -which are categorized as the quasi-elastic e-scattering, p-scattering, and DIS by the types of the primary signature -for the benchmark model are summarized in Appendix A. In this section, we compare these primary cross sections in a wide range of parameter space. Then, we take the detector effects into account and expect the preferred primary signatures in the reference experiments for any given parameter set within the range of our consideration.

p-scattering vs. DIS
We first contrast the (quasi-elastic) p-scattering with DIS in this scenario. In our terminology (as in much of BDM literature), the p-scattering means that the incoming light dark matter χ 1 scatters off proton elastically, χ 1 p → χ 1 p (eBDM), or inelastically, χ 1 p → χ 2 p (iBDM), in which the proton does not break apart, hence not accompany additional hadronized objects. In reality, it is probabilistic, competing with other scattering channels including DIS, and can be given by a function of the modulus of the spatial momentum of the recoiling proton, p p ≡ | p p |. In this paper, we assume that it is a step function around p p = 2 GeV for simplicity, i.e., the events with p p < 2 GeV is categorized as p-scattering.
Here the value p p = 2 GeV is the boundary value for which the probability of producing a pion or a charged secondary in the water is 50% based on the Monte Carlo study by the SK Collaboration [89]. 4 As we will argue shortly, the differential cross section in p p peaks around p p 1 GeV in the parameter space of interest (see also for example, FIG. 2 of Ref. [34]), and thus the precise p p value differentiating the DIS and p-scattering regimes does not significantly affect our analysis results. In addition, we require energy threshold E th for any recoil proton, whereas no corresponding cuts are imposed on visible particles in the DIS process as it may involve unnecessary complication. Indeed, as we shall show, the p-scattering is more important than the DIS in most of parameter space of interest. Therefore, this seemingly "unfair" treatment on the p-scattering nevertheless does not affect our final conclusion.
In the first two panels of the upper row of Figure 2, we show the contours of σ cut χ 1 p /σ DIS in the m 1 − m X plane for E 1 (= m 0 ) = 10 GeV (upper-left) and 50 GeV (upper-middle) fixing m 2 = m 1 (i.e., eBDM). Here σ cut χ 1 p denotes the primary p-scattering cross section with the energy threshold cut imposed. We specifically require the recoil kinetic energy of the proton to be above 21 MeV [31,90] which would be adopted in a "DUNE-like" detector in the next section. In calculating the DIS scattering cross section, σ DIS , we take the MSTW2008NNLO parton distribution functions [91] and require the energy transfer Q > 1.5 GeV. We confirm that our σ DIS reproduces the shape of the differential cross section of neutrino-induced DIS in Ref. [92] after replacing the dark gauge boson by a W boson. For the neutrino scattering, it is well known that the DIS dominates over the quasi-elastic scattering (and Σ resonance) once the energy of incoming neutrino is large enough, E ν > O(10) GeV [92]. On the other hand, we find that the p-scattering cross section in BDM scenarios dominates over the DIS cross section even when the energy of incident χ 1 is as large as 50 GeV (and up to 100 GeV which is not shown here), as long as the mediating particle (i.e., the dark gauge boson in our study) is lighter than O(1) GeV. Regarding the fact that we imposed a "penalty" to p-scattering cross section by an NDIS/yr/kt/Y, E1=50 GeV, m2=1.5m1 Figure 2. Upper panels: Contours of σ cut χ1p /σ DIS expected at a "DUNE-like" detector. The first two plots show the results for eBDM (i.e., m 2 = m 1 ) with E 1 being 10 GeV (leftmost) and 50 GeV (middle). The rightmost panel displays the σ cut χ1p /σ DIS = 1 lines corresponding to various δm values, with E 1 being 50 GeV. Lower panels: The expected number of DIS events. The first two plots are for eBDM (i.e., m 2 = m 1 ) with E 1 being 10 GeV (leftmost) and 50 GeV (middle), while the last one is for iBDM with m 2 = 1.5m 1 and E 1 = 50 GeV. E th cut, our statement is rather robust. In the upper-right panel of Figure 2, we display the results corresponding to non-zero δm, i.e., iBDM processes. Contours are defined by σ cut χ 1 p /σ DIS = 1 with E 1 = 50 GeV. The overall behavior is similar to the cases with δm = 0. However, the regime where the p-scattering is in favor widens as δm/m 1 increases, because an up-scattering process is easier with proton itself than with a parton inside the proton.
We can semi-analytically understand these behaviors of p scattering and DIS in this regime in the following manner. First of all, the differential p-scattering cross section in proton recoil momentum is peaking toward small p p ( m p ) due to the t-channel exchange of X (see Appendix A for details). From Eqs. (A.1) through (A.4), we obtain in the limit of p p m p . Here we omit the overall factor from the phase space integral and terms in the numerator of the relevant matrix element. This relation implies that the p-scattering cross section rises in decreasing m X ( m p ≈ 1 GeV) where the quartic dependence is kept as long as p p m X . By contrast, the energy transfer Q in the DIS should be larger than ∼ 1.5 GeV, and in turn, much larger than m X under consideration. The DIS differential cross section is given by for m X 1 GeV (see also Appendix A.2). This implies that the DIS cross section does not vary much for m X 1 GeV. Our numerical study further suggests that σ DIS be less than σ cut χ 1 p for m X ≈ 0.1 GeV and E 1 50 GeV, and therefore, the ratio of σ cut χ 1 p to σ DIS should be greater than 1 for m X 1 GeV. Moving onto the lower panels of Figure 2, we now show the contours of the theorylevel number of DIS-induced events N DIS per year·kt·Y with Y defined by (atomic number)/(atomic weight). We calculate N DIS as follows: where F 1 is the χ 1 flux coming from the GC, given in Eq. (2.2), N i is the number of protons or neutrons inside a fiducial volume of a detector, and t exp is the amount of time exposure. The leftmost and middle panels show results for eBDM case with E 1 = 10 GeV and 50 GeV, respectively. The dark-sector coupling g 11 is assumed unity, and the kinetic mixing parameter is set to be 10 −4 for illustration. By contrast, in the rightmost panel, we set m 2 = 1.5 m 1 together with g 12 = 1 and = 10 −4 for E 1 = 50 GeV. For simplicity, the contours are shown only in the parameter region that σ cut χ 1 p /σ DIS 1, i.e., where σ DIS can be at least comparable or larger. Nevertheless, we observe that N DIS does not change much when m X 1 GeV and decreases in increasing m X and m 1 in the entire parameter space. The former two results are consistent with our argument that the differential cross section dσ DIS /dxdy ∝ 1/(Q 2 + m 2 X ) 2 . The last result (dependence on m 1 ) is due to the fact that the incoming χ 1 actually scatters off a parton instead of a proton for the DIS case.
In case of HK which uses water as the target material, i.e., Y = 10/18, with the fiducial volume of 380 kt, the contour N DIS · yr −1 kt −1 Y −1 = 10 −3.3 corresponds to 0.1 yr −1 . Since this is the expected number without considering any detector effects such as energy threshold, angular and position resolutions, the actual N DIS would be negligible even after running the HK for more than 10 years. Consequently, it is generically hard to observe DIS-induced signal events unless a much larger detector or flux of χ 1 is introduced. 5 Reminding that the iBDM cross section is further suppressed for DIS, as stated previously, the aforementioned conclusion applies both for eBDM and iBDM signals.
In summary, we find that the p-scattering is larger than the DIS as far as the mediator X is lighter than ∼ 1 GeV even in the case of E 1 ≈ 50 GeV. Furthermore, we anticipate that it is unlikely to observe DIS-initiated signal events in the near future belonging to the parameter regime where the DIS cross section dominates over the p-scattering cross section, modulo the χ 1 flux given by Eq. (1.1). These series of observations are truly contrasted to neutrino-initiated events. The main difference stems from the fact that the weak gauge bosons W and Z involved in neutrino scattering processes are much heavier than typical dark gauge boson considered in this work. In this context, one may wonder whether or not examining the DIS channel is motivated essentially to probe the model space where m X is sizable (say, a few tens of GeV). In the vector portal scenario, an 2 suppression is unavoidable in any scattering channels, on top of the suppression by ∼ 1/m 4 X . Therefore, the signal sensitivity relevant to the DIS channel, in general, quickly drops as reduces. Again this is clearly supported by our exercise shown in the lower panels of Figure 2. We henceforth focus on the (quasi-elastic) p-scattering and e-scattering processes throughout the rest of this paper unless specified otherwise.

e-scattering vs. p-scattering: basic considerations
Given our observation in the previous section that the contribution from the boosted χ 1induced DIS is negligible as far as a light mediator (much lighter than say, W boson mass) is concerned, we can now compare the e-scattering (i.e., involving an electron recoil) with p-scattering (i.e., involving a proton recoil). The beginning discussion is devoted to pure theoretical results which would have been obtained with "perfect" detectors; that is, we include neither characteristics of a detector such as resolutions and energy threshold nor nuclear effects. Although this setup is unrealistic, the relevant exercise allows us to develop intuitions and paves a road to truth-level understanding. Once done with the theoretical investigation, we discuss, in the next section, how the preferred channels are affected by the inclusion of detector effects such as energy threshold, angular resolution, and the acceptance of the secondary signatures. Elastic BDM processes require the first effect as they accompany a target recoil only, whereas inelastic BDM ones are affected by all three. To demonstrate the difference according to detector types, we consider two representative ones, a DUNE-like LArTPC detector and a HK-like Cherenkov light detector.
We show our findings with respect to E 1 , the incoming energy of χ 1 ; again, results will be relevant for the same E 1 value, irrespective of the production mechanism of boosted χ 1 . Figure 3 demonstrates comparisons between e-scattering cross section and p-scattering cross section for four example E 1 (or equivalently the mass of χ 0 in the annihilating twocomponent dark matter scenario) values: E 1 = 0.5, 1, 5, and 50 GeV in the upper-left, upper-right, lower-left, and lower-right panels, respectively. Results are shown in the plane of the χ 1 mass, m 1 , vs. the dark gauge boson mass, m X , with seven different mass gaps from δm = 0 (red) to δm = 50m 1 (purple). Note that the case with a vanishing mass gap corresponds to the elastic BDM scenario. We essentially scan over the two-dimensional m 1 − m X parameter space, and at each scan point we compute the total scattering cross section within allowed phase space. 6 For the proton channel, we consider it up to recoil proton spatial momentum being 2 GeV beyond which the proton may mainly break apart and give rise to a DIS, as stated previously. While performing the scan, we observe that the proton scattering cross section is always (at least, slightly) greater than the one corresponding For the proton channel, the phase space within p p = 2 GeV is calculated, if phase space is allowed beyond it. In most of parameter space, the proton scattering is showing a greater cross section than the electron one is, so the contours are drawn along σ χ1e = 0.9σ χ1p . The proton (electron) channel becomes advantageous (comparable) as m X and m 1 increase (decrease). The black dot-dashed lines are defined by m X = 2m 1 . See the text for more detailed discussions.
to electron in the scanning ranges of interest, unless m X ∼ a few keV. So, we choose to divide the region to the p-preferred and the e-comparable along the boundary defined by σ χ 1 e = 0.9 σ χ 1 p . Above (below) the boundary, the p-scattering (e-scattering) channel comes with more (comparable) number of signal events.
Together with the boundary curves, we add black dot-dashed diagonal lines to mark We show which processes (Feynman diagrams) are allowed (check mark) or not allowed (cross mark) due to the dark sector mass spectrum in each scenario. Gray horizontal lines mean that the phenomenology is relatively independent of the respective process. The column on the right provides a comment on possible ways of having a visible signature in each scenario.
the border defined by m X = 2m 1 . Above the line an on-shell X decays invisibly into dark matter pair χ 1χ1 , whereas below the line X is allowed to decay to visible SM particles assuming that no lighter dark-sector particle exists. We henceforth call the former and the latter scenarios I and II. In both scenarios, the mass spectrum of the dark sector will dictate how and under which circumstances the visible decays of χ 2 can take place. We list all possibilities that give rise to three visible particles below and provide a pictorial summary in Figure 4.
• Scenario I-ii: The mass spectrum satisfies 2m 1 < m X and δm < m X , but 3m 1 < m 2 . If g 11 is suppressed or vanishes (i.e., a model maximizes the off-diagonal gauge interactionχ 2 γ µ χ 1 X µ as in the examples in Ref. [37]), the channel of χ 2 → 3χ 1 is not allowed, hence a visible decay is possible.
• Scenario II-i: The mass spectrum satisfies m X ≤ 2m 1 and m X < δm (or ≤ m 1 + m 2 if g 11 = 0). χ 2 emits an on-shell X as a decay product, and the X decays visibly.
• Scenario II-ii: The mass spectrum satisfies δm < m X ≤ 2m 1 (or ≤ m 1 + m 2 if g 11 = 0). χ 2 decays visibly via a three-body process just like Scenarios I-i and I-ii.
Looking into plots in Figure 3, first of all, we see the trend that the region, where the e-scattering channel stays competitive, expands as δm becomes smaller (i.e., less inelastic) and incident χ 1 comes with more energy. The revealed dependence on the former is not surprising because smashing a (lighter) χ 1 on a (heavier) proton target is much more advantageous in transiting to a much heavier χ 2 state, considering the maximally allowed m 2 for a given pair of E 1 and m T (T = e, p). The relation is given by which can be approximated to The relation in (3.5) corresponds to the electron target case (i.e., m T = m e ), whereas that in (3.6) corresponds to the proton target case (i.e., m T = m p ). These formulas for the two limiting cases clearly support the fact that for a given E 1 the proton target is better to access much larger values of m 2 than the electron target. Now defining δm = r m 1 , taking the electron target, and solving the above inequality (3.4), we obtain the below inequality for m 1 in term of E 1 and r: This implies that vertical drops take place below the saturation point defined in the righthand side (henceforth called kinematic "barrier"). In reality, as E 1 increases, the actual drop happens close to this kinematic "barrier". When it comes to the case of lager E 1 , the value in (3.7) increases so that the e-comparable area extends toward higher m 1 . There is another barrier in the m X direction mainly due to the faster drop of σ χ 1 e in increasing m X , compared to σ χ 1 p . Indeed, our numerical study finds that the p-scattering cross section is roughly constant or independent of E 1 within the parameter space of interest, whereas the e-scattering one is enhanced in rising E 1 . Therefore, the "e-comparable" area is widened along the direction of m X in growing E 1 . We now close this section, summarizing the main observations as follows (see Table 1).
• If a BDM search hypothesizes a heavy dark gauge boson (but still not much above the GeV scale), the proton scattering channel expedites discovery.
• If a model conceiving iBDM signals allow for large mass gaps between χ 1 and χ 2 , the proton channel is more advantageous.
• On the other hand, the electron channel becomes comparable/complementary in probing the parameter regions with smaller m 1 and m X .
• As the boosted χ 1 comes with more energy, more parameter space where the escattering is comparable opens up.
Again we emphasize that even if these observations are predicated upon simplified calculations, they are not too far away from the reality. As we will find out in the next section, many of generic trends are retained even in the presence of various realistic effects. Furthermore, depending on detector designs, some effects may be negligible; for example, energy threshold could be extremely small, highly collimated objects could be resolved, and so on. We will come back to this issue and make a brief discussion in Section 4.

e-scattering vs. p-scattering with realistic effects
In order to compare the e-scattering and the p-scattering channels with detector effects taken into account, a detector type is first specified because different ones are characterized by different dimensions, energy threshold, resolution, and so on. Our study in this section begins with a LArTPC detector. More specifically, we employ a "DUNE-like" experiment, meaning that it has four modules each of which comes with a cubic-shaped 17.5 kt total volume formed by {width × height × length} = {15.1 m × 14.0 m × 62.0 m}. In our analysis, we assume a fiducial volume of 10 kt for each module in a "DUNE-like" experiment. More complete DUNE detector specifications and properties are collected in Table 2 in Appendix B.
We essentially scan over a set of mass points in the m 1 − m X plane for a given pair of E 1 and m 2 /m 1 . The values of m 1 are varied from 10 −4 GeV to 10 −0.5 GeV with 70 points equally spaced in the logarithmic scale, while m X ones are varied from 10 −4 GeV to 10 0 GeV with 80 points equally spaced in the logarithmic scale. Namely, we perform simulation for 5,600 points. In each point, we generate 5 million events using the TGenPhaseSpace class implemented in the ROOT package and reweight them with respective matrix element values. In this study, we consider the case that the secondary vertex accompanies only an e − e + pair for simplicity. However, a sequence of procedures described later on are straightforwardly applicable to other channels. For the final state particles, we require the following set of criteria: ii) ∆θ e−i > 1 • , ∆θ p−i > 5 • with i denoting the other visible final state particles, and iii) both primary and secondary vertices should appear in the detector fiducial volume. 7 We warn the reader that the minimum spatial momentum for proton (i.e., p th p = 200 MeV) corresponds to E th = 21 MeV, predicated upon the assumption that the level of E th in ArgoNeuT [90] can be achievable in DUNE. The last criterion is applied if any intermediary particle (either χ 2 or X) is long-lived. However, we do not calculate the laboratory-frame decay length event-by-event because it would take an enormous amount of time. To save Figure 5. Comparisons between the e-scattering and the p-scattering channels at a DUNE-like detector, including various realistic effects such as energy threshold, angular separation, acceptance with respect to the displaced vertex. Contours are the numbers of signal events calculated by Eq. (3.8) with statistics of 40 kt·yr. The upper row is for E 1 = 0.5 GeV while the lower row is for E 1 = 5 GeV. From the left to the right panels, δm are set to be 0, m 1 , and 10m 1 , correspondingly. Red and blue dotted contours show the expected number of signal events in the proton and electron channels with 40 kt·yr statistics. Red-shaded and blue-shaded regions denote p-preferred and epreferred ones, respectively, and the boundaries between them are given at σ χ1e = σ χ1p . The white regions are where decay process χ 2 → χ 1 e + e − is not kinematically allowed, while the gray regions are where some or all of the final-state particles are not isolated or detected. See the text for more detailed explanations. the time and computing resource, we instead adopt a "shortcut" but rather conservative strategy. A detailed description can be found in Appendix A.3.
The results with E 1 = 0.5 GeV are shown in the upper row of Figure 5. Let us first take a look at the leftmost one corresponding to the eBDM case. Due to the fact that the proton differential cross section dσ/dp p is quickly rises toward smaller p p and the cuts harder on the proton, the fiducial cross section of the e-scattering now can be the same as that of the p-scattering unlike the comparisons in the previous section. So, we divide the entire space of interest into the e-preferred (blue shaded) and the p-preferred (red shaded) regions by σ χ 1 e = σ χ 1 p . Just for developing the intuition on the expected number of signal events, we add contours by red and blue dotted curves corresponding to p-scattering and e-scattering events, respectively. The number of events, N sig , can be calculated by where F 1 is the signal flux shown in Eq. (2.2), t exp is the amount of time exposure, and N p(e) is the number of total target protons (electrons) inside a detector fiducial volume. We assume statistics of 40 kt·yr. Here A denotes the final signal efficiency including A (see Appendix A.3). As in the previous section, the diagonal black dot-dashed line describes the boundary beyond (below) which an on-shell dark gauge boson X predominantly decays invisibly (visibly). Therefore, any limit or interpretation with respect to the mass points above and below the line should be associated with corresponding exclusion limits in the dark gauge boson mass vs. kinetic mixing parameter. We will come back this point in Section 4. Basically, the e-preferred region expands, compared to the corresponding one displayed in the upper-left panel of Figure 3. As argued in Appendix A, this is because the distribution of the differential cross section dσ/dp p of the p-scattering is steeper than that corresponding to the e-scattering (see also Ref. [34] for generic shapes of proton and electron recoil energy spectra). This implies that the application of the E th cut could reject the region around which dσ/dp p is peaking, in particular when m X gets much smaller than m p [see Eq. (A.13)] and m 1 is too small to boost up the recoiling proton to overcome the energy threshold. Possible search strategies inferred from this exercise are not much different from what is summarized in the previous section; the electron channel is better in probing smaller m 1 and m X , and the other way around. However, care should be taken here. Although many points in the above parameter space allow for quite a few signal events, eBDM suffers from large backgrounds such as atmospheric neutral current neutrino scattering events since it has only target recoil in the final state. Therefore, search strategies targeted at (seemingly) well-motivated parameter points should be designed with the potential of background contamination taken explicitly into account.
When it comes to the iBDM scenario, relevant phenomenology becomes even richer. Now all selection criteria i) through iii) are imposed. The third one regarding the displaced vertex requires explicit values for kinetic mixing parameter and dark-sector gauge coupling g 12 . Unfortunately, these factors are not completely canceled out in the ratio of σ χ 1 e to σ χ 1 p for a given {E 1 , m 1 , m 2 } because the decay length also depends on the boost factor of the slowly decaying particle which is in turn a function of target mass as well. So, we fix and g 12 to be 10 −4 and 1 throughout this section for illustration, but the expected results with other and g 12 choices differ not much (i.e., the same intuition should go through).
We start by looking at the white regions. Since we focus on the case of an electronpositron pair coming out of the secondary vertex, the immediate parent particle for them should be massive enough to create them. This is translated as follows. If χ 2 decays to χ 1 and an e − e + pair via an off-shell X (i.e., m 2 < m 1 + m X ), the mass gap between χ 2 and χ 1 , δm should be greater than 2m e . The relevant region of δm < 2m e is shown by the wider one sweeping through the space vertically. 8 By contrast, if χ 2 decays to an on-shell X along with χ 1 (i.e., m 2 > m 1 + m X ), the dark gauge boson mass should be at least twice the electron mass to disintegrate to the e − e + . 9 This case of m X < 2m e is mapped to the horizontally-lying white region at the bottom of the plot. Another type of white region arises in the case of δm = 10m 1 (see the rightmost panel). As discussed before, there is the kinematically allowed maximum m 2 value for a given set of E 1 , m 1 , and target mass m T . Although the proton scattering generously allows for a wide range of δm, it also encounters the associated kinematic "barrier" which is actually realized in the right-hand side of the plot.
The gray-colored regions are somewhat different in the sense that iBDM events can be produced but some or all of final state particles are not isolated or detected. For the vertically-stretching gray region, the associated mass spectra marginally allows for the creation of the e − e + pair in the χ 2 decay, but they are too soft to get over the threshold for electron. However, for the horizontally-extended gray region, the (too) light X is so significantly boosted that its decay products e + and e − are too close to each other, hence not isolated by criteria ii). The remaining narrow gray band shown in the upper rightmost panel of Figure 5 emerges because of the threshold for protons. Remember that nearby space is already close to the kinematic "barrier" for the p-scattering process. As a result, most of incoming E 1 is spent for forming a massive χ 2 , but only a tiny fraction of energy is transferred to a target, inducing a soft recoiling proton.
Speaking of diagonal lines, there arises a different type denoted by the black dashed. One can easily see that with respect to the line χ 2 decays to χ 1 e + e − via from on-shell X to off-shell X as moving bottom to top, or vice versa. Therefore, the signal acceptance may be affected substantially due to this kinematical transition, along the line.
Just like the eBDM case, the e-preferred region becomes wider than the corresponding theoretical prediction again because cuts on energy and angle diminish the signal acceptance in the p-scattering channel. The lower vertical drop arises because the up-scatter of χ 1 becomes kinematically challenging beyond it. Indeed, the associated N sig contours show that the e-scattering cross section falls rapidly near the border. By contrast, the vertical boundary toward bigger m X was absent in the corresponding theory prediction, but now emerges. The reason for this lies in the angular separation between the e + and e − from the χ 2 decay. Since the proton is usually heavier than the other particles involved in the process, the proton recoil is not likely to carry out a large amount of energy, but χ 2 is. Thus, χ 2 is significantly boosted so that its decay products are inclined to be too collimated to satisfy the angle cut. Of course, a similar behavior is anticipated in the e-scattering channel. However, the recoiling electron and the e − e + pair are more likely to share the incoming energy E 1 (together with the outgoing χ 1 ). Three electron tracks are often moving in the same direction, but some fraction of events can pass the angle cut. We tion [37], but it requires an introduction of additional particles to induce the relevant operator. We do not examine this possibility for the sake of minimality. 9 We do not consider the lighter on-shell X dominantly decaying to three photons, due to the Landau-Yang theorem [93], by the 8-dimensional operators with four field strength tensors. In this case, the mean decay time of X is too long to provide a secondary signature in a terrestrial detector and various cosmological and astrophysical bounds can apply [94]. Comparisons between the e-scattering and the p-scattering channels at a HK-like detector, including various realistic effects such as energy threshold, angular separation, acceptance with respect to the displaced vertex. Contours are the numbers of signal events calculated by Eq. (3.8) with statistics of 380 kt·yr. See the caption in Figure 5 for more details.
see that this sort of e-preferred region becomes even wider as δm increases (see the upper rightmost panel of Figure 5). Basically, a wider range of low-mass m 1 is accessible so that the resultant low-mass χ 2 will give a merged e − e + pair in the proton channel.
We remark that care should be taken in interpreting the results in the rightmost panels where m 2 = 11m 1 > 3m 1 . In this case, the result above m X = 2m 1 line would lose its meaning, if g 11 in Eq. (2.3) were sizable. χ 2 would decay invisibly to three χ 1 particles. Therefore, the comparison in that regime is relevant to models with vanishing or suppressed g 11 (see Scenarios I-ii and II-ii described in Section 3.2).
We next show the results corresponding to E 1 = 5 GeV in the lower panels of Figure 5. Similar behaviors and trends are essentially retained as in the case of E 1 = 0.5 GeV. A larger value of E 1 simply opens more accessible phase space in the χ 1 mass, so one can understand the results as some sort of "stretching" of those in the E 1 = 0.5 GeV case. In summary, the e channel is complementary to or even favored over the p channel toward smaller m 1 or smaller m X . On the other hand, the p scattering will be a discovery channel when a search is targeting at the regime where χ 1 and X are heavier. 10 Finally, we emphasize that background-free searches for iBDM are possible as relevant signal events come up with several unique features. Hence, hunting for iBDM signals will provide complementary information in investigating dark-sector physics.
Lastly, we discuss how the phenomena that we have observed at a DUNE-like LArTPC detector are affected by different detector specifications. As advertised before, we take a "HK-like" Cherenkov detector with a fiducial mass of 380 kt, as an example. Again we refer to Table 2 in Appendix B for detailed specifications and properties. Similarly to the DUNE-like detector, we demand the final state particles to satisfy the following conditions: i) p e > 100 MeV, 1.07 GeV < p p < 2 GeV, 10 The expected number of events in the p-scattering channel looks small because a small value of and one-year data collection are assumed. One can reach the discovery-level events by increasing the hypothesized value of and/or data collection.
ii) ∆θ e−i > 3 • (∆θ e−i > 1.2 • ) for p e < 1.33 GeV (p e > 1.33 GeV) and ∆θ p−i > 3 • for all p p with i running over the other visible final state particles, and iii) both primary and secondary vertices should appear in the detector fiducial volume.
Note that the minimum required energy to observe an electron in SK/HK is expected to be ∼ 5 MeV, but the associated angular resolution is too poor. So, we instead adopt 100 MeV to keep a decent level of angular resolution as per the search for elastic BDM interacting with an electron performed by the SK Collaboration [35]. One should also notice that a momentum threshold of p p 1.07 GeV, which corresponds to E th 485 MeV, is required for the relatively heavy proton to emit Cherenkov light. Finally, the acceptance A associated with criterion iii) is evaluated in the same way as explained earlier.
The comparisons between the e-scattering and p-scattering channels are shown in Figure 6. We report the results only with E 1 = 5 GeV, because E 1 = 0.5 GeV is so small that entire events in the proton channel would not be selected due to the threshold for protons. Hence, the associated comparison is meaningless. Nevertheless, the results with E 1 = 5 GeV will suffice to discuss phenomenological differences. Looking at the eBDM case first (the leftmost panel), we see that the e-preferred region is significantly extended in comparison with the corresponding theory prediction in Section 3.2. This is because a proton should have enough kinetic energy to create Cherenkov radiation while the proton scattering cross section is typically peaking toward small energy/momentum transfer especially for m X < ∼ O(1) GeV. Therefore, the electron channel may be considered more important in the search for BDM signals. When it comes to iBDM cases (see the middle and the rightmost panels of Figure 6), similar understanding goes through. A crucial difference is the fact that the gray regions become much wider for a given E 1 . This is purely due to the larger thresholds and angular resolution. Therefore, HK-like Cherenkov detectors are not ideal for probing parameter space with small m X and/or small m 1 through the iBDM channel. However, at the same time, this implies that devising a search strategy getting around the threshold and resolution issues would enable to explore the above-mentioned regime.

Example Data Analysis
Furnished with various observations and guidelines discussed in the previous section, in this section we demonstrate how actual analyses would be conducted, focusing on iBDM signals. Since our benchmark model described in Eq. (2.3) involves interactions with a dark gauge boson X, it is natural to study expected experimental sensitivities in dark gauge boson parameter space with respect to the four benchmark experiments listed in Section 2.2. We essentially contrast parameter coverages for several reference model points. By doing so, we develop the intuition, what experiment and channel would be better motivated or more advantageous for a given reference point. Conversely, each experiment would design optimized search strategies and choose best-motivated parameter space, performing similar exercises demonstrated here.
A parameter scan is done in the following way. We first fix values of E 1 , m 1 , and m 2 , and then divide the expected m X coverage into 50 segments equal-sized in the logarithmic scale. For each parameter point defined by (E 1 , m 1 , m 2 , m X ), one million events reweighted by the matrix element are generated and energy and angle cuts are applied in order to calculate a cut-related acceptance. We next scan space from 10 −6 to 3 × 10 −3 after chopping the range into 350 intervals equal-spaced again in the logarithmic scale. For each , the maximum laboratory-frame mean decay length of the long-lived particle (either χ 2 or X depending on the mass spectrum)¯ max lab is calculated with g 12 set to be unity for simplicity. Plugging¯ max lab to the relevant fitted function, we extract an A value. The final acceptance A can be obtained by multiplying the cut-related acceptance by A . One can calculate the number of expected signal events, utilizing Eq. (3.8). Here N p(e) and t exp can be easily determined, once one chooses a detector and the amount of time exposure. Cross sections (before imposing any selection criteria) are calculable according to the formulation detailed in Appendix A, and the flux value for a given E 1 (= m 0 ) comes from Eq. (2.2) in the annihilating two-component dark matter scenario. Labeling the 90% confidence level (C.L.) exclusion limit by N 90 which depends on the number of expected background events, we calculate it with a modified frequentist construction [95,96]. An experiment is said to be sensitive to a given signal if N sig ≥ N 90 , i.e., parameter points whose N sig is greater than N 90 will be excluded with 90% C.L. under the associated background assumption. We hereafter take a zero-background assumption for simplicity, as iBDM-induced events come with distinctive features that known SM processes (e.g., atmospheric neutrino-initiated ones) are hard to mimic (see also discussions in Section 2). 11

DarkSide-20k vs. DUNE
Our first exercise is to compare DarkSide-20k and DUNE experiments. Before showing relevant results, we need to explain our event selection scheme for DarkSide-20k. It is basically designed to be best sensitive to the signatures induced by the scattering process between slowly-moving WIMP or WIMP-like dark matter and a nucleon. The scale of typical recoil kinetic energy is expected to be order of 1 − 100 keV, so that detectors are most sensitive to this energy range. In particular, the ionization signal (called S2) is sort of "amplified" in the gas-phase Argon to measure such a small energy deposit. However, in many well-motivated cases, E 1 (or equivalently m 0 ) is of sub-GeV range, and as a consequence, the typical energy scale of final state visible particles in our signal processes will be several tens of MeV or greater. A potential problem with a large energy deposit is that photomultiplier tubes -which are the agents to measure the amount of deposited energymay be saturated due to too much amplification beyond their measurement capacity, and thus the associated data analysis would face severe challenges. To avoid this potential issue, we here take a rather tricky strategy to look at only scintillation signal (called S1). 12 The iBDM process under consideration often involves a displaced secondary vertex as discussed 11 Nevertheless, one may argue with a few plausible possibilities such as neutrino-induced resonance and DIS events involving a handful of mesons which subsequently decay to charged particles. A more systematic discussion will be elaborated in Ref. [76]. 12 We expect that iBDM-induced S1 are rather strong unlike WIMP-induced ones.
in Section 2. Therefore, it would be possible to identify two separated vertices by some dedicated S1-pattern analysis. 13 For illustration, we assume that a ≥ 30 cm separation is discernible at the DarkSide-20k detector. 14 So, the corresponding acceptance A at the DarkSide-20k detector is calculated by checking if the secondary vertex is not only within the associated detector fiducial volume but separated from the primary vertex by more than 30 cm. We find that the resultant A is well accommodated by a different, empirical fit template 15 : where c i (i = 1, 2, 3, 4) are fit parameters. No other cuts associated with energy and angle are imposed, so that A can be directly translated to final signal acceptance A.
The results are displayed in Figure 7, where we compare experimental sensitivities which would be achieved in the DarkSide-20k and DUNE experiments in the plane of m X vs. . The brown-shaded covers the current excluded regions whose boundary values are extracted from Refs. [101] and [102] for Scenario I and Scenario II, respectively. We assume 1-year exposure for concreteness and apply the same selection criteria listed in the previous section, for DUNE. The upper-left panel is for Scenario I, contrasting DarkSide-20k-20t·yr (red lines) and DUNE-40kt·yr (blue lines) in the proton (solid lines) and electron (dashed lines) channels for a reference point defined by (E 1 , m 1 , m 2 ) = (250, 5, 15) MeV. We see that DarkSide-20k can cover a decent range of parameter space in spite of 2,000 times smaller fiducial mass than DUNE. Being aware that the scattering cross section is quadratic in and A decreases as 2 , we estimate that a 2,000 times smaller target mass would lead about 6.7 (≈ 2000 1/4 ) times degradation in probing values. This sort of rough estimation goes through the electron channel, as DUNE can explore down to ∼ 10 times smaller for a given m X . E th for electrons at DUNE is sufficiently low relative to the chosen E 1 and effectively no minimum required length for displaced vertices (due to its mm-scale position resolution) unlike DarkSide-20k together with a bigger detector length scale. However, when it comes to the proton channel, the reach at DarkSide-20k is degraded only by a factor of ∼ 2 at m X = 0.01 GeV. The reason for this is that no energy and angular cuts are imposed for DarkSide-20k whereas an application of sizable energy threshold and angular separation cuts away a certain fraction of events at DUNE. On the other hand, recoil proton gets harder in increasing m X so that threshold affects less and the difference in the coverage gets gradually larger. Another phenomenon to find here is that the proton channel is preferred over the electron one in both DarkSide-20k and DUNE experiments. In the previous section we observe that the p-scattering is more advantageous as E 1 decreases. Especially, in order for DarkSide-20k to attain some signal sensitivity, smaller E 1 values 13 Here we do not claim or justify at all that this way of analysis is viable, which is certainly up to experimentalists' effort and not within the scope of this paper. However, we bring readers' attention to e.g., DEAP-3600 [97][98][99] in which an interaction vertex is identified only with a set of S1 signals.
14 For example, DEAP-3600 [97][98][99] claims a resolution of ∼ 10 cm for a single vertex. As there are two vertices in a signal event, we believe that our assumption here is reasonable. 15 It is inspired by the fit template introduced in Ref. [100]. 2)]. Therefore, the search for BDM signals in the proton channel is better motivated at DarkSide-20k.
In the lower-left panel, we examine the expected experimental sensitivity for E 1 = 100 MeV. Because of energy threshold, DUNE has no sensitivity to this reference point. Since a smaller E 1 implies a larger χ 1 flux, the resultant parameter coverage becomes wider. We further investigate the dependence of experimental sensitivities on the vertex resolution (denoted by d res ), employing two more different values, d res = 90 cm and d res = 250 cm. We see that d res = 90 cm still allows to explore a similar range of parameter space. By contrast, d res = 250 cm rapidly reduces the regions to probe because d res becomes comparable to the detector length scale of DarkSide-20k.
Moving onto the experimental sensitivities for Scenario II, we demonstrate our results for (E 1 , m 1 , m 2 ) = (250, 50, 60) MeV, in the upper-right panel of Figure 7. Since the chosen parameter values forbid χ 1 from up-scattering to a χ 2 of 60 MeV in the electron channel, we report only results in the proton channel. Unlike the previous case, DarkSide-20k can explore a wider range of parameter space than DUNE can, mainly due to the threshold for proton. As we argued in Section 3, recoil proton is likely to be softer in decreasing m X so that many events at DUNE do not pass the E th criterion for proton. Again this clearly shows how the energy threshold affects the experimental sensitivity.
An interesting phenomenon worth mentioning is that there exists a finite range of m X near m X = δm for which no signal sensitivity is essentially available in both experiments. But the underlying reasons for the individual ones are quite different. Note that χ 2 decays to an on-shell X below the m X = δm line so that the potentially long-lived particle is X (χ 2 ) below (above) the line. For a given value, χ 2 (undergoing a three-body decay) is much more long-lived than the on-shell X. For DarkSide-20k, the on-shell X usually decays rather "promptly" with respect to the given d res , and as a result, two vertices are not resolvable. However, in decreasing m X , the on-shell X becomes significantly boosted, so that the laboratory-frame life time can be significantly dilated. Hence, DarkSide-20k restores the signal sensitivity for lighter m X . When it comes to DUNE, the electron and positron coming from the on-shell X are not energetic enough to overcome the energy threshold. For the chosen set of parameter values, the scattered χ 2 is not much boosted. Just below m X = δm, the produced X is (roughly) as boosted as χ 2 so that its energy is no more than m X × E 1 /m 2 ≈ 40 MeV with m X = 10 MeV and the maximal boost of χ 2 assumed (i.e., γ 2 < E 1 /m 2 ). Therefore, each of the positron and electron carries away ∼ 20 MeV which is below the threshold. However, as m X gets smaller, the electron (positron) emission direction relative to its own boost direction and/or the boost direction of X comes into play. Therefore, a certain fraction of events can pass the selection criteria and the signal sensitivity starts to grow again.
We next show the expected sensitivity for E 1 = 100 MeV in the lower-right panel of Figure 7 just as in the case of Scenario I and see the potential of probing a wider range of parameter space. The dependence of experimental sensitivities on the vertex resolution is also investigated for completeness. We again find that experimental sensitivities mildly depend on d res unless it is comparable to the detector length scale of DarkSide-20k.
One should be mindful that in this analysis we did not include a potential contribution from the coherent scattering of χ 1 off target nuclei. The reason is that the energy of incoming χ 1 is 100 MeV or more so that we expect that the majority of events happen via an incoherent scattering process. Nevertheless, it is fair to say that the chance that the momentum transfer (to the target) falls below O(10 MeV) may not be negligible even for E 1 = 250 MeV. 16 In this context, our reference E 1 values may be in a regime far from concluding either a full coherent or full incoherent scattering. However, we believe that our approach using the incoherent scattering cross section σ χ 1 p is rather conservative since the coherent nuclear scattering would be enhanced by the atomic number of a target nucleus for dark gauge boson scenarios.

DUNE vs. HK
The next exercise is to contrast DUNE and HK experiments. The total mass of HK detectors is about 10 times larger than that of DUNE far detectors. On the other hand, the DUNE detectors are expected to fulfill preciser measurements thanks to the LArTPC technology. Therefore, it is expected that they both would show a similar level of performance in the search for iBDM signals for a given amount of time exposure. Our main results are displayed in Figure 8.
First of all, the upper-left panel shows sensitivities with respect to Scenario I, which would be achievable in the electron (dashed lines) and proton (solid line) channels in the DUNE (blue lines) and HK (green line) experiments, for a reference point defined by (E 1 , m 1 , m 2 ) = (500, 5, 15) MeV. Event selections for DUNE and HK are done according to selection criteria described in Section 3.3, and again 1-year data collection is assumed.
Note that here the proton channel is not available in HK because E 1 is too small for any single recoil proton to overcome the relevant threshold. Comparing the electron channel, we see that both experiments show similar coverage. The chosen m 1 is close to m e so that the four final state particles (i.e., e − , e − , e + , and χ 1 ,) tend to equally share the incoming energy. Therefore, in a large fraction of events, all three electrons have an energy greater than 0.1 GeV. This implies that other choices of m 1 would result in smaller coverage of HK because of harder selection cuts. Now comparing the e-scattering and p-scattering in DUNE, one should notice that the latter wins over the former as m X increases whereas the trend is reverse toward smaller m X . This agrees with the observation made in Sections 3.2 and 3.3 that the e-scattering channel becomes competitive as m 1 and m X decrease.
We next raise E 1 , m 1 and m 2 , reporting the associated result in the lower-left panel of Figure 8. Increasing E 1 essentially means a reduction of the χ 1 flux, so one would naively expect that HK shows wider coverage due to its 10 times bigger mass. We find, however, that selection cuts indeed matter over the detector fiducial mass. Relatively larger energy threshold often forces events to fall in the collimated regime where the angular separation plays a crucial role, in turn. On the other hand, keeping a sufficient angular distance often forces (some) final state particles to be soft relative to threshold. Therefore, a right balance between the two factors, in general, would be desired like the previous case in order to accomplish better experimental sensitivities. Now the proton channel becomes available because E 1 is large enough to allow for p p > 1.07 GeV. However, the same, above-described argument applies to the proton scattering more stringently, and as a result, it enables us to probe only excluded regions.
We perform similar exercises for Scenario II and present the results in the right panels of Figure 8. As before, the upper and lower ones are for E 1 = 0.5 GeV and E 1 = 2 GeV, respectively. Basically, similar interpretations are relevant here. Unlike the study in DarkSide-20k vs. DUNE, no discontinuity around the line of m X = δm arises because both E 1 values are sufficiently large. Nevertheless, we observe some sort of transition about the line. Interestingly, sensitivities in the proton channel get worse in decreasing m X . For the p scattering, χ 2 typically takes away more energy than the proton because m p m 1 , m 2 . This implies that the on-shell X, which is now long-lived, becomes more boosted in decreasing m X , so that not only resultant A drops rapidly but the e − e + pair from the X decay becomes more merged. In particular, HK is affected by the latter effect more than DUNE. On the contrary, χ 2 in the electron channel is not as much boosted as that in the proton channel so that the above-mentioned effects are not substantial.

HK vs. DeepCore
As the final exercise, the DeepCore experiment is compared to the HK experiment. Unlike DUNE and HK, DeepCore does not see charged particle-induced tracks (except muons) but measures cascade-type signatures. Therefore, if all three final state particles come out of a single point (within the detector vertex resolution), it is (almost) impossible to distinguish it from a true recoil-only event. To get around this issue, we take a rather simplified strategy; any event showing double (displaced) cascades will be tagged as signal.
(See Ref. [103] for a double-cascade signal in the context of heavy neutrino-accompanying processes.) A τ neutrino would give rise to a double-cascade event, but it would have to carry (at least) sub-PeV energy which is much away from the energy scale of interest here. From a comparison between IceCube and DeepCore detector arrays, we estimate that ∼ 5 meters distance is resolvable in DeepCore. The following selection criteria are applied: i) p recoil e > 10 GeV, p secondary e + e − > 10 GeV, and ii) the secondary vertex should appear in the detector fiducial volume and be at least 5 meters away from the primary vertex.
The second requirement in i) means that the energy deposit by the e − e + pair should be collectively greater than 10 GeV threshold. The main spirit behind criterion ii) is similar to that for DarkSide-20k, so we adopt the same strategy to calculate the corresponding A .
In this analysis, we slightly modify the selection criteria for HK as well; for a fairer comparison we do not demand an angular separation of 3 • for the electron-position pair emitted from the secondary vertex. As we will show shortly, the E 1 value under consideration is as large as order several tens of GeV. In such an energy regime, potential backgrounds are expected scarce so that some selection cuts may be relaxed appropriately. However, a more detailed discussion is not a focus of this work. We simply assume that the aforementioned selections and modifications will allow for a zero-background analysis. 17 Unlike the previous two exercises, we report the results only in the electron channel. In particular, any (isolated) recoil proton cannot exceed the energy threshold of the DeepCore experiment. Furthermore, considering the required energy scales relevant to DeepCore, the proton-involved (incoherent) scattering processes are not dominant. The χ 1 -induced DIS would come into play. According to our discussion in Section 3.1, however, we may have to collect the data for > ∼ 10 years in much of the parameter space of interest in order to reach the expectation of having a single event, before multiplying the acceptance. Nevertheless, it is worth to look into DIS-induced events with large values of E 1 as the effective mass of DeepCore is at least ∼ 5 Mt. We leave the relevant dedicated analysis to the future [104].
We finally demonstrate experimental sensitivities for Scenarios I and II in the left and right panels of Figure 9, respectively. DeepCore-5Mt·yr and HK-380kt·yr are color-coded by purple and green. In the left panel, we first contrast the sensitivities with (E 1 , m 1 , m 2 ) = (50, 0.01, 0.03) GeV. The incoming energy is enough to overcome the threshold in DeepCore. On the other hand, the associated χ 1 flux is not large, so DeepCore is better motivated due to its ∼ 13 times larger effective mass. However, we see that as smaller E 1 (here E 1 = 25 GeV) is assumed, DeepCore quickly loses signal sensitivity because of the energy threshold, while HK achieves improved sensitivity due to a higher χ 1 flux.
Similar phenomenon can be observed for Scenario II (see the right panel). For larger E 1 , DeepCore would be better motivated, but for smaller E 1 hypotheses, the sensitivity of DeepCore falls rapidly, motivating detectors which are sensitive to smaller energy deposit.
Here we see that DeepCore has no signal sensitivity below the m X = δm line. The reason is that in this regime an on-shell X is created and decays quickly within 5 meters, i.e., no double-cascade signature arises. The chosen mass spectrum indeed does not let χ 2 decay late even when m X > δm. So, the acceptance A at HK does not change drastically across the m X = δm line and both of the HK curves do not show an appreciable kink structure around the line.
Before moving on to the concluding section, we make a couple of noteworthy comments. Since DeepCore comes with an effective mass of 5 Mt, a few years of data taking could make it sensitive to parameter space with E 1 > ∼ 1 TeV, which is far beyond the reach of other experiments. More interestingly, such a large E 1 allows χ 1 to up-scatter to a much heavier state so that visible particles other than an electron-positron pair can be emitted from the secondary interaction vertex (see Section 2). For example, the e-scattering channel with E 1 > 100 GeV may accompany µ + µ − in the final state. This is certainly an interesting enough possibility, as it may induce a signature of a cascade plus a track (from different vertices). We reserve the subject for future work [104].

Conclusions
The search for boosted dark matter has received rising interest as an alternative approach to probe dark sector physics including cosmological dark matter. Several dark matter model frameworks have been proposed in order to give rise to boosted dark matter in the universe today: for example, two-component dark matter scenario [25,26,46,47], Z 3stabilized dark matter models carrying semi-annihilation processes [48], models involving dark matter-induced nucleon decays [49], models with decaying super-heavy particles [29,30,50], or cosmic-ray induced energetic dark matter scenarios [51][52][53]. Many ongoing and future dark matter direct detection and neutrino experiments can observe signals induced by boosted dark matter, and a host of phenomenological studies have proposed search strategies, channels, and sources of boosted dark matter [25, 27-34, 36-38, 44-47, 49-53, 76, 77], with regard to those experiments. The relevant field of research is rather new and systematic approaches toward devising clever search schemes are absent.
A more efficient design to dig out elusive boosted dark matter signals is of great importance especially at the earlier stage of experiments, improving sensitivity to BDM models.
In light of this situation, we performed a dedicated study to address the issue above, taking a benchmark model, namely, a two-component dark matter scenario with boosted (lighter) dark matter scattering off either electron or proton. We allowed the boosted dark matter χ 1 to scatter off not only elastically (i.e., eBDM) but inelastically (i.e., iBDM) to a heavier unstable dark-sector state χ 2 which further decays back into the lighter dark matter and additional secondary (visible) particles. For each of the two possible scenarios, the scattering cross sections of boosted dark matter via electron and proton targets were carefully investigated. Deep inelastic scattering cross sections were included in association with the proton scattering. 18 We have found that the DIS contribution is negligible in the total proton scattering cross section, as far as the mass of the particle mediating the interaction between boosted dark matter and SM particles is not significantly larger than the typical energy scale to induce a DIS. A parameter choice for which the DIS dominates over the proton scattering cross section typically results in a tiny total cross section so that the proton channel becomes irrelevant unless the detector of interest comes with a huge volume to compensate for the small signal rate. Therefore, treating the proton itself as a target suffices when one considers the cross section of boosted dark matter scattering-off a proton. We emphasize that the generic conclusions here readily hold for other scenarios involving boosted/relativistic dark matter (e.g., scenarios with cosmological dark matter decaying to lighter dark matter, fixed target experiments), as our study was done with respect to the energy of incoming boosted dark matter.
We then compared electron scattering and proton scattering channels in both of elastic and inelastic boosted dark matter signals. We first performed simplified comparisons assuming a perfect detector (i.e., perfect energy/angular resolutions, no energy thresholds, etc.) simply to develop our intuition. We found that the proton scattering channel is better motivated as hypothesized m 1 and m X increase. The proton channel is also advantageous if the underlying model assumption is an iBDM process with a large mass difference between χ 1 and χ 2 . We then took some realistic effects into account such as energy threshold, angular resolution, and tagging efficiency for the displaced vertex. To show the dependence on detector types, we performed separate studies on a DUNE-like LArTPC detector and a HK-like Cherenkov detector. For the former type, most of the observations made in the simplified study are relevant. On the other hand, for the latter type, the electron channel becomes competitive in a broader range of kinematically allowed parameter space. We, however, observed that (relatively) hard cuts associated to the HK-like detector significantly reduce not only the range of accessible parameter space but the expected number of signal events in the search for iBDM. Therefore, it seems that ameliorating experimental performance relevant to the event selection is an inevitable task toward improved signal sensitivities.
We finally conducted example analyses in our benchmark detectors: DarkSide-20k, DUNE, Hyper-K, and DeepCore. Experimental data, which would be collected at those 18 Recently, a DIS module for eBDM is implemented in GENIE [105].
detectors, was interpreted in terms of coverage in the plane of dark gauge boson mass vs. kinetic mixing parameter, as our benchmark model to describe interactions between χ 1 and SM particles involves a dark gauge boson. We found that relevant experimental sensitivities highly depend on various detector specifications such as geometry, energy threshold, resolutions, and so on. Thus, an optimized (BDM or BDM-like signal) search scheme can be established up on multilateral considerations about those factors. Detectors with a small volume and small threshold (e.g., DarkSide-20k) are better for searching for signatures induced by less energetic boosted χ 1 , whereas those with a large volume and large threshold (e.g., DeepCore) can be sensitive to signal events coming from high energetic boosted χ 1 although the associated χ 1 flux is often small. A preciser measurement/reconstruction of final state visible particles can help smaller-volume detectors achieve experimental sensitivities comparable to larger-volume detectors.
where δm ≡ m 2 − m 1 . In the laboratory frame, the differential scattering cross section is where E T = m T + E 1 − E 2 is the energy of the recoiling target and λ(x, y, z) = (x − y − z) 2 − 4yz. The kinematically allowed maximum and minimum of E T , defined as E + T and E − T , are where p χ 1 = E 2 1 − m 2 1 . Let us explain how to calculate the terms in M 1 × (form factors), following Ref. [106]. In the matrix element, the vertexū(p p )γ µ u(p 0 ) should be replaced byū(p p )Γ µ u(p 0 ), where p p is the 4-momentum of the recoiling proton while p 0 is the 4-momentum of the proton at rest before the scattering. We then express Γ µ in terms of F 1 and F 2 as follows: where q ν is the transferred energy and momentum, i.e., p ν 1 −p ν 2 = p ν p −p ν 0 , and the anomalous magnetic moment κ = 1.79 (-1.91) for the proton (neutron). For the proton, the form factors are related with the Sachs electric and magnetic form factors G E and G M such that for Q 2 = −q 2 = 2(E p − m p )m p . Note that q 2 is negative for E p > m p due to the metric (1, −1, −1, −1). The Rosenbluth formula in conjunction with experimental measurements shows the following dipole approximation up to around Q 2 ∼ 10 GeV 2 with µ p = (1 + κ p )e/(2m p ) = 2.79 [107]. 19 The form factors are also related as (A.10) 19 Note, however, that the ratio of the electric and magnetic form factor µpGE/GM decreases from 1 for larger Q 2 > 1 GeV 2 from the high precision double polarization experiments [108,109]. In our analysis, we assume that keeping the dipole approximation for Q 2 4 GeV 2 is reasonable in the region of elastic scattering of our proton target. Note that the term p e / m 2 e + p 2 e is almost one as long as p e m e . Hence, the differential cross section peaks around p e m e . Moreover, we numerically find that the distribution is rather smooth for m e p e E 1 , in contrast to that of dσ χ 1 p /dp p . Interestingly, for E 1 m 1 , δm and p e m e , m 1 , the value of the differential cross section is almost flat up to p e = E 1 + as long as it is smaller than the maximum E max e in Ref. [34]. The flat region becomes wider in increasing E 1 .

A.2 Deep inelastic scattering
We can write the amplitude squared for deep inelastic χ 1 + N → χ 2 + X, as where D µν and W µν are the tensors related to the χ 1,2 and hadronic vertices, respectively. We call p 1,2 the four-momentum of χ 1,2 and q µ ≡ p µ 1 −p µ 2 is the momentum transfer. Taking into account non-vanishing masses for the dark sector fermions, we obtain D µν = g µν q 2 − (m 1 − m 2 ) 2 − 2(p µ 1 q ν + q µ p ν 1 ) + 4p µ 1 p ν 1 , (A.20) while we parameterize (k µ is the four-momentum of the target nucleon and k 2 = M 2 ) The structure functions W 1,2 are written in terms of the parton distribution functions (PDFs): where x is the momentum fraction carried by the parton and f (x) is the associated PDF.
A sum over all partons is implicit in Eq. (A.22). Proceeding with a further calculation we obtain

A.3 "Short-cut" scheme
We first prepare an acceptance table for the 10 kt DUNE-like detector module as a function of mean decay lengths of arbitrary long-lived particles. To this end, an arbitrary coordinate (regarded as a primary vertex) is generated in the detector fiducial volume, a random direction is generated, 20 and then an actual flight distance is also generated according to the decay law with respect to a given mean decay length. From the coordinate, the direction, and the actual flight distance, the secondary vertex point is identified. If it is within the fiducial volume, it is accepted, otherwise rejected. For each of 10,000 decay length values, 100 million pairs of the coordinate and the direction are produced to obtain the acceptance (denoted by A ). Not all decay length values are probed, but we have found that the decay length vs. inverse acceptance is rather well described by a combination of a second order polynomial and a straight line jointed at some point which can be determined by fit: where c i (i = 1, 2, 3, 4) denote fit parameters. Therefore, the acceptance for any noninvestigated decay lengths can be calculated with this empirical fit function. Now for each mass point, we calculate the maximum laboratory-frame mean decay length of the associated long-lived particle (denoted by¯ max lab ), and find the corresponding acceptance by feeding¯ max lab to the fit function. Finally, the eventual signal efficiency is given by the product of efficiencies attached to the above-bulleted items.

B Information of Detectors
Here we collect specifications and useful information of detectors adopted at various largevolume ( > ∼ 1 ton) dark matter and neutrino experiments, and tabulate in Table 2. One can take advantage of this table in combination with the main messages of our data analyses, when designing experimental strategies for searching for BDM-induced or similar signatures.  Table 2. A summary of detector specifications in various large-volume dark matter and neutrino experiments. The numbers in the Active Volume entry for DUNE, SK, and HK are the values of their total volume. Note that p p = 1.07 GeV corresponds to E th 485 MeV for SK/HK. [PID: particle identification, LAr/LXe: liquid argon/xenon, DP/SP: dual/single phase, TPC: time projection chamber, S1: prompt primary scintillation, LS: liquid scintillator, m.w.e.: meter water equivalent].