Inelastic Dark Matter at the Fermilab Short Baseline Neutrino Program

We study the sensitivity of the Fermilab Short-Baseline Neutrino (SBN) experiments, MicroBooNE, ICARUS, and SBND, to MeV- to GeV-scale inelastic dark matter interacting through a dark photon mediator. These models provide interesting scenarios of light thermal dark matter, which, while challenging to probe with direct and indirect detection experiments, are amenable to accelerator-based searches. We consider production of the dark sector states with both the Fermilab Booster 8 GeV and NuMI 120 GeV proton beams and study the signatures of scattering and decay of the heavy excited dark state in the SBN detectors. These distinct signatures probe complementary regions of parameter space. All three experiments will be able to cover new ground, with an excellent near-term opportunity to search for cosmologically motivated targets explaining the observed dark matter abundance.

A host of disparate gravitational phenomena provides seemingly incontrovertible evidence for dark matter (DM), yet it is striking how little we know about its fundamental properties. This stark situation has spurred extensive explorations for novel theories and new experimental probes of DM over the past decade. While many interesting ideas have emerged, the hypothesis of a light dark sector weakly coupled to the Standard Model through a portal interaction is particularly compelling. Viable models of thermally produced light DM below the Lee-Weinberg bound [1] can be realised in this framework [2][3][4]. Experiments at the intensity frontier have a vital role to play in probing the low masses and feeble couplings characteristic of such theories. In particular, proton beam fixed-target experiments, including accelerator neutrino beam experiments, offer an interesting testing ground for new light, weakly coupled particles and DM . Due to the large collision luminosity and forward lab-frame kinematics, a substantial forward flux of relativistic dark sector states may be produced in the primary proton-target collisions. These dark particles may subsequently be observed through their decays or scatterings as they traverse a detector located downstream of the target. These experiments provide a critical component of a broader experimental program of searches at medium-energy e + e − colliders, fixed-target missing energy/momentum experiments, low mass direct detection experiments, rare particle decay experiments, and the LHC, which will provide powerful sensitivity to light dark sectors and DM in the MeV-GeV mass range [30][31][32][33][34].
A particularly promising near-term venue for proton fixed target dark sector searches is provided by the Fermilab Short-Baseline Neutrino (SBN) program. The SBN program comprises three liquid argon time projection chamber (LArTPC) detectors -SBND [35], Mi-croBooNE [36], and ICARUS [37] -installed along the Booster Neutrino Beam (BNB) [38]. Interestingly, dark sector particles can be produced both from the medium energy (8 GeV) on-axis Booster proton beam and from the higher energy (120 GeV) Main Injector NuMI proton beam. Concerning the latter, MicroBooNE and ICARUS are located approximately 8 • and 6 • off-axis with respect to the NuMI beam, implying that the associated flux of dark particles passing through these detectors can be substantial. Several studies have highlighted the exciting sensitivity of the SBN experiments to light exotic new particles, including long-lived heavy neutral leptons [39], DM tridents [40], Higgs portal scalars [41], and elastically scattering DM [42], and the first experimental searches have been carried out by MicroBooNE [43,44]. With SBND coming on-line soon and ICARUS now operating, significant improvements in reach are expected over the arXiv:2106.04584v1 [hep-ph] 8 Jun 2021 next several years.
In this work we investigate the near-term prospects for the SBN experiments to probe models of MeV-GeV scale inelastic dark matter (iDM) coupled to the SM through a kinetically mixed dark photon. First proposed in the context of the DAMA anomaly [45,46], iDM models as advocated in [16,47] have received broader attention in recent years as a compelling scenario for sub-GeV thermal DM [19,20,23,[48][49][50][51][52][53][54][55][56][57][58][59][60][61]. The off-diagonal nature of the DM coupling to a heavier excited state provides a simple mechanism to evade otherwise stringent bounds from direct detection experiments and cosmic microwave background observations. On the other hand, simple iDM models also feature new potential signatures associated with DM produced in accelerator experiments, including inelastic up-and down-scattering and semi-visible decays of the excited state to DM and a pair of charged leptons. While the sensitivity of accelerator-based neutrino experiments to iDM has been studied in other contexts [16,19], here we consider the specific experimental issues related to their production in the Booster and NuMI beam lines and their detection in the SBN LArTPC detectors. In particular, this is to our knowledge the first time that the sensitivities of the three upcoming SBN experiments to iDM scenarios have been estimated. Furthermore, we have for the first time singled out and quantified the neutrino-induced background relevant for long-lived decay for experiments with angular reconstruction capability, as well as estimated the background rejection efficiency of angular cuts. In particular, searches for onaxis dark sector production from the Booster at SBND and off-axis production from NuMI at ICARUS can explore significant new regions of parameter space, including those predicting the observed DM relic density.
The paper is structured as follows. In Section II we present the iDM model to be studied in this work, including a discussion of the long-lived excited dark state and the DM relic abundance. Our modelling of the iDM signal and neutrino-induced backgrounds at the SBN experiments, as well as our analysis strategy, is discussed in Section III. Our projections for the SBN sensitivity to the iDM model are presented in in Section IV, and we provide our conclusions and outlook in Section V.

II. LIGHT INELASTIC DARK MATTER
In this section we discuss a simplified model of sub-GeV iDM interacting via a dark photon mediator that can be probed by Fermilab SBN experiments. After introducing the model, we discuss the thermal relic abundance predictions and the existing experimental constraints on this scenario.

A. Model
We consider a model of fermionic iDM based on a spontaneously broken U (1) D gauge symmetry with a massive dark photon mediator, V µ [47]. 1 The dark photon interacts with the SM via kinetic mixing with the ordinary photon. The Lagrangian describing the dark photon is where M V is the dark photon mass, ε is the kinetic mixing strength, and F µν is the SM photon field strength. The kinetic mixing generates a coupling of the dark photon to the electrically charged particles with strength suppressed by ε. In the physical basis, the dark photon interactions with the SM are thus described by the coupling ε eV µ J µ EM , where J µ EM is the SM electromagnetic current.
The dark matter sector consists of two Majorana fermion states χ 1 and χ 2 , with Lagrangian Here M χi denotes the physical fermion masses, and we assume M χ1 < M χ2 such that χ 1 is the stable DM candidate. Two important parameters that will appear in our phenomenological considerations are the fractional mass splitting ∆ χ , defined as and the dark photon-to-DM mass ratio R χ : It is clear that M χ1 , ∆ χ and R χ parameterise the spectrum of the dark sector. The dark sector particles are assumed to dominantly couple to the dark photon mediator through an off-diagonal interaction, where g D ≡ √ 4πα D is the U (1) D gauge coupling. The simplified model presented here can be realised as the low energy limit of a renormalisable completion involving a dark Higgs field with a U (1) D charge of 2 and a Dirac fermion field with U (1) D charge of 1. Gauge and Yukawa couplings of the dark Higgs field generate both a dark photon mass term and small Majorana masses for the Weyl components of the Dirac fermion after spontaneous U (1) D breaking. Provided the Dirac mass is much larger than these Majorana masses, the dominant interaction of the DM is given by equation (5). Diagonal couplings of the form V µ χ i γ µ γ 5 χ i are suppressed by the ratio of the difference in Majorana masses to the Dirac mass. Due to an enhanced charged conjugation symmetry, it is technically natural for the Majorana mass difference to be small or even vanishing [51]. We note that the dark Higgs boson cannot in general be completely decoupled from the spectrum at large α D ∼ 0.5 while maintaining perturbativity [18]. However, as long as its mass remains of the order of the dark photon mass it does not significantly affect the phenomenology [23], and thus it will thus not be considered in this work.
The iDM scenario enjoys several attractive features. First, for M χ1 + M χ2 M V , the observed DM relic abundance can be obtained through thermal freeze-out of DM co-annihilation to SM particles for certain choices of model parameters, providing predictive cosmologically motivated targets for experiments. This will be discussed in detail below in Section II B. Second, the heavier Majorana fermion state χ 2 decays rapidly for ∆ χ M χ1 > 2m e , implying that DM annihilation is strongly suppressed for temperatures below the mass splitting. Thus, stringent constraints on late-time DM annihilation from precision studies of the cosmic microwave background temperature anisotropies do not apply to this scenario [64,65]. Lastly, the non-relativistic scattering of DM off of SM particles in direct detection experiments, χ 1 SM → χ 2 SM, faces a strong kinematic suppression even for very small mass splittings due to the inelastic (endothermic) nature of the transition 2 . Therefore, accelerator experiments can provide a unique probe of this simple and well-motivated sub-GeV DM scenario.

B. Relic density for iDM
The underlying mechanisms relevant for the cosmological production of light iDM have been thoroughly studied in the literature. Our goal here is to obtain a simple recasting procedure of the relic density targets obtained from the full evolution of the Boltzmann equations from Refs. [16,50,54] for generic values of ∆ χ and R χ . The relic density arises from χ 1 χ 2 co-annihilation to SM particles, which can be resonantly enhanced when This requires a modification of the standard instantaneous freeze-out approximation [73]. In particular, we will build our recasting procedure in the resonant region following the approach developed in [74].
The effective thermally-averaged cross-section takes the form: where we have included the number of degrees of freedom of each DM component, and the second term reads: with R = M 2 V (Mχ 1 +Mχ 2 ) 2 − 1 and (σv) lab the total DM annihilation cross-section times velocity to any SM states χ 1χ2 → SM. The first case corresponds to the resonant limit with Γ V M V where most of the annihilation occurs directly on the resonance, while the second is the standard result for a s-wave cross-section. In order to include all the potential hadronic channels, we use the R-ratio approach [47], defining which we have extracted from DarkCast [75,76]. The final relic density scales as: where g eff is the SM effective number of degrees of freedom evaluated at the freeze-out temperature x f . We can estimate x f following the standard equation [73]: with m pl = 1.22 × 10 19 GeV the Planck mass and g eff the SM effective number of degrees of freedom. Note that as shown in Eq. (6), σv eff contains another exponential factor of x f . Based on the above ingredients, we can obtain a simple estimate of the relic density. Since numerous recent works have presented full solutions of the Boltzmann equations for different choices of parameters [16,50,54], we will instead use the ingredients above to interpolate between the full relic density lines of [54]. 3 In more details, the steps are the following.
• Solve Eqs. (9) and (10) for both the reference relic density line and for the new values ∆ χ and R χ , leading to (Ωh 2 ) ref , (Ωh 2 ) new , and to the ratio: where we used the latest DM relic density estimate from the Planck collaboration [77], and since our estimate is based on the instantaneous freeze-out the ratio is in general not precisely 1, particularly for large ∆ χ and strongly resonant cross-sections.
• Determine the final relic density target in ε f from solving (Ωh 2 ) new = 0.12 × r ref .
In practice, the above procedure is thus simply a "smart" interpolation between the existing full calculation results using the instantaneous freeze-out results. We have crosschecked our interpolations with the independent results from Ref. [50].

C. Existing constraints
A large variety of intensity frontier experiments have the potential to search for iDM-like signatures. Analyses from past experiments have been re-interpreted in recent years, and new results from modern experiments have also been presented. We list below the most important ones along with the relevant limits we include in our final result.
a. CHARM The CERN-Hamburg-Amsterdam-Rome-Moscow experiment, in addition to its main on-axis detector, used an off-axis calorimeter module to search for new long-lived particles, including heavy neutral leptons [78,79] and axion-like particles [80]. These limits have been re-interpreted in the context of various dark sector models in the literature. In particular, the limits from CHARM on the iDM scenario have recently been considered in [53]. Note that we extend the previous re-interpretation by further including the direct parton-level production, improving the CHARM limit for the large mass regime. We have simulated parton-level production of dark sector states using the same strategy as for the SBN program described in Sec. III B.
Although the original search from CHARM relied on observing two charged particles, the analysis did not resolve each track. Instead they required that the observed pulse height in a scintillator plate, placed directly in front of the main calorimeter, should be larger than 1.5 minimum ionizing particles (corresponding to 9 MeV). Coupled with the requirement that no more than 2 tracks should be seen in the calorimeter, these cuts lead to no detected events.
b. NuCal Using the 70 GeV proton beam from the U70 accelerator, NuCal accumulated 2 · 10 18 protons on target (PoT) [81,82]. Similarly to CHARM, it featured a large decay volume and could therefore operate in a near-zero background mode. We use the limits based on [53].
c. LSND Using the Los Alamos Meson Physics Facility 0.8 GeV proton beam, the LSND experiment accumulated ∼ 10 23 PoT [83]. A measurement of neutrinoelectron elastic scattering analysed in [83,84] has been recasted to place limits on the iDM model, and we use the limits derived in [16] and [53] depending on the model parameter sets used. It is worth mentioning that this result was derived from the reinterpretation of neutrino scattering results, assuming that the e + e − pair would be reconstructed as a single electron for small enough angular separation with a similar efficiency.
d. NA64 The NA64 experiment, a fixed-target experiment employing a CERN SPS secondary 100 GeV electron beam, has accumulated 2.84 × 10 11 electron-ontarget (EoT), with plans for collecting 5 × 10 12 EoT in the coming years. Based on an active target, it searches for large missing energy due to a dark photon escaping the target [85]. It mostly constrains the low mass region.
e. BaBar Based on the e + e − PEP-II B-factory, the BaBar experiment has accumulated a total luminosity of 53 fb −1 , which was analysed for single photon events with large missing energy [86]. This analysis places strong constraints on dark photons that decay invisibly to DM particles. However, for iDM with large mass splittings, the heavy state χ 2 instead decays semi-visibly within the detector, which significantly weakens the bounds for large dark photon masses [54,87].
f. Experiments sensitive to DM scattering The MiniBooNE-DM collaboration has recently used the MiniBooNE detector in an "off-target" run configuration to search for light elastic sub-GeV DM [88,89] as described in more detailed below. For large dark gauge coupling α D ∼ 0.5, the corresponding limit is competitive with NA64 and we will include it when relevant. We also note that re-interpretations of analyses from other past experiments sensitive to DM-electron scattering, such as BEBC [42,90], E137 [91][92][93][94], NOvA [22,95], can have competitive sensitivity. See Ref. [42] for a recent study and comparison.

III. SEARCHING FOR IDM AT THE SBN PROGRAM
The SBN experiments at Fermilab include the Short-Baseline Near Detector (SBND), MicroBooNE, and ICARUS. These detectors all employ liquid argon time projection chambers (LArTPCs) to detect tracks of charged particles. The detectors are placed along the Booster Neutrino Beam (BNB), which is generated by striking a 1.7 nuclear interaction length beryllium target with 8 GeV kinetic energy protons. The BNB also has the capability to run in "off-target" mode, during which the proton beam is steered past the beryllium target and directed onto the iron absorber at the end of the decay volume, 50 m downstream of the target. In off-target mode the beam-related neutrino flux is significantly reduced compared to the normal running mode, allowing for enhanced sensitivity to new physics signals. This off-target run configuration was used by the MiniBooNE-DM collaboration to search for light elastic sub-GeV DM [88,89]. While there are currently no plans to perform an off- . OA indicates "off-axis" and λ is the interaction length of the corresponding material. The distance (to the centre of the experiment) D and typical detector length L are further indicated. Note that the geometry used is the active volume, around 25 − 30% larger than the fiducial volume.
target run with the SBN experiments, we will investigate the potential sensitivity of such a run to the iDM model below In addition to this beam, the ICARUS detector is less than 6 degrees off-axis from the Neutrinos at the Main Injector (NuMI) beam, generated by 120 GeV protons impacting a 2-interaction-length graphite target. We thus also consider the sensitivity of the ICARUS detector to off-axis production of iDM in the NuMI beam (ICARUS-OA) We have summarised the various experimental setups considered in this work in Table I.

A. iDM signatures in neutrino experiments
Dark sector particles interacting through the vector portal are typically produced in proton beam dump experiments through the following three different mechanisms [97]: • Pseudoscalar meson decay. Low mass dark photons may be efficiently produced through the decays π 0 , η, η → γV .
• Proton bremsstrahlung. Dark photons can be directly produced via bremsstrahlung, pp → pV X, which dominates for M V ∼ 1 GeV [97,98] • Drell-Yan. For higher beam energies and masses, dark sector particles can be directly produced via the parton-level process qq → V → χ 1 χ 2 .
The iDM model leads to several potential signatures in proton fixed-target experiments including up-scattering, down-scattering, and semi-visible χ 2 decay [16]. Which signature dominates is primarily dictated by the lifetime of the excited state.
• Long-lived. If the lifetime of χ 2 is long enough for it to reach the detector, it may then either downscatter to a DM state χ 1 or decay to a DM particle χ 1 and visible SM particles through an off-shell dark photon. In this case the scattering signatures (both χ 1 e → χ 2 e and χ 2 e → χ 1 e) would give rise to a subdominant contribution with respect to the decay signatures.
• Quasi-stable. In the region where the χ 2 lifetime is too long, the dominant signatures are χ 1 e → χ 2 e and χ 2 e → χ 1 e, which both mimic the standard DM elastic scattering signature.
• Short-lived. If χ 2 is too short lived to reach the detector the only signature is the up-scattering χ 1 e → χ 2 e. Depending on the χ 2 lifetime, χ 2 can decay back into the neutrino detector, leading to a displaced vertex signature. Otherwise the upscattering mimics the standard DM elastic scattering.
Hence there are three distinct regions corresponding to three characteristic phenomenological signatures: (i) upscattering followed by decay χ 2 in the detector, (ii) decay of beam-produced χ 2 in the detector, and (iii) scattering. In order to understand which signatures dominate in different regions of the parameter space we need to study the χ 2 decay width. In the regime M V M 1 ∆ χ M 1 > 2m e , the excited state decays via χ 2 → χ 1 e + e − . The partial decay width for this process is given by [47] The strong dependence on the small mass splitting leads to macroscopic decay lengths even for relatively large couplings. For the experiments of interest here, the lab frame decay length is .
In order to fully discuss the iDM phenomenology, a comment regarding the decay modes of χ 2 is necessary to assess both the pure decay and the upscattering followed by decay phenomenology. For mass splittings less than the dimuon threshold, ∆ χ M χ1 < 2M µ , the excited state decays via χ 2 → χ 1 e + e − . For larger χ 1 − χ 2 mass splittings, new leptonic and hadronic channels open up. We have included these additional channels based on the R-ratio approach, obtaining the differential hadronic decay rate χ 2 → χ 1 + (hadrons) via: where √ s is the momentum of the off-shell dark photon (as done in e.g. [55]). We present the leptonic branching ratios for several values of the splitting parameter ∆ χ in Fig. 1. In particular, we find that the branching ratio suppression for the e + e − final state is always fairly mild, reduced by at most a factor of 3 for the largest masses accessible in beam dump experiments. We note that in the opposite regime of very small splittings below the di-electron threshold, ∆ χ M χ1 < 2m e , the excited state decays via χ 2 → χ 1 + 3γ and is stable on cosmological time scales [63,66].
Overall, Eq. (14) provides insight into the characteristic phenomenology for given parameter choices, as we now discuss in further detail.
a. Long-lived decay χ 2 → χ 1 e + e − For small splittings, small kinetic mixing and intermediate range masses the dominant signature consists of a beam-produced excited state χ 2 , which is long lived enough to travel to the detector and decay semi-visibly. A particularly important aspect of the considered three-body decay χ 2 → χ 1 e + e − is that the angle of the final e + e − is correlated with the energy of the pair and the mass splitting ∆ χ . More precisely, since the total exchanged squared momentum in the off-shell dark photon s satisfies s < (∆ χ M χ1 ) 2 we can decompose s in the lab frame in terms of the energy and angle of the e + e − pair, leading to [16] where, in the laboratory frame, we used an angular cut θ c on the e + e − angle and and energy cut E c on both e + and e − . For small angular separation and a typical energy threshold E c = 30 MeV, we easily get an approximation for the lower accessible mass threshold: b. Down-scattering and up-scattering signatures We now turn to the region of very small ∆ χ and small M χ1 , for which both up-scattering and down-scattering signatures resembles those of regular elastic scattering since the χ 2 is quasi-stable on detector length scales. In this parameter region our focus will be on scattering with electrons, which yields a clean forward energetic single electron signature. The cross section for both up-and down-scattering is the same as for standard elastic scattering up to small corrections of order ∆ χ M χ1 /E and ∆ χ M χ1 /m V [16]. The approximate differential cross section with respect to the outgoing electron energy E f in the laboratory frame is [92] dσ where E is the incoming χ 1 or χ 2 energy and f (E f ) = 2m e E − m e E f + m 2 χ + 2m 2 e . c. Up-scatter followed by decay Finally, for the large ∆ χ and large ε region, the main signature is χ 1 upscattering followed by a prompt or displaced χ 2 decay. We focus on χ 1 -nucleon up-scattering in the large ∆ χ region for a couple of reasons. First, up-scattering off electrons is kinematically suppressed for large mass splittings, particularly for DM produced with the lower energy BNB. The χ 1 energy threshold E min to initiate an up-scattering on a target particle of mass m is We see that for large splittings the typical energy for upscattering off electrons is much greater than for scattering off nucleons. Second, while elastic DM-nucleon scattering searches typically must contend with large neutrinoinduced backgrounds [88,89], for the case of iDM each up-scattering process can lead to an associated visible decay vertex, providing an additional handle to reject the backgrounds and leverage the large DM-nucleon scattering rate. We have included both the above electronscattering cross-section (for the small splitting regime) as well as up-scattering on nucleons via the public code BdNMC [97], adapting the implemented cross-section to the Dirac fermion case using the results of Ref. [10]. As the visible charged particles from the χ 2 decay are expected to constitute the main signal, we have not implemented cuts on the outgoing nucleon kinematics. Our projections for the expected number of events should be considered as conservative estimates, particularly for ICARUS in an off-axis configuration where the typical DM energy ranges between 5 to 30 GeV. Note that we do not include the potentially relevant deep-inelastic scattering processes, but these would also be worth exploring in a full study.

B. Signal simulation chain
It is clear from the discussion above that an accurate description of our signal events requires a complete simulation of the iDM production at the target and its subsequent propagation to the SBN detectors, followed by the decay and/or scattering processes. The main steps of our simulation chain are presented below.
For meson-induced production in the BNB, we simulate π 0 and η meson distributions using the Sanford-Wang empirical distribution [99] using the fits made by the MiniBooNE collaboration [100] as implemented in BdNMC [97]. When considering the BNB-OT and the NuMI-OA configurations, we rely instead on the meson distributions stored in the database of Ref. [101,102]. The latter are based on a GEANT4 [103] simulation of the meson production in the current NOvA target for NuMI-OA and in the iron dump used during the Mini-BooNE run [100] for BNB Off-target. For larger masses, we use the proton bremsstrahlung production as implemented in BdNMC [97], which includes resonant ρ/ω vector meson mixing using the time-like form factor derived from [104]. When relying on parton-level production, we instead use the MadDump [24] plugin for Mad-Graph5 aMC@NLO [105] to generate the distribution of both χ 1 and χ 2 from an s-channel dark photon.
An important subtlety is that the proton bremsstrahlung process as implemented in BdNMC relies on the assumption that where E beam is the incoming proton energy, m p the proton mass, E V the outgoing dark photon energy in the lab frame, and p V T its transverse momentum [97]. For the on-axis BNB beam, we use the cut For the NuMI-OA, the situation is peculiar in that one needs p V T ∼ 0.1E V in order to point toward the off-axis ICARUS detector. In order to satisfy Eq. (20), we use All production mechanisms are then centralised to a modified version of BdNMC [18,23], which further includes the propagation of the dark states to the detector and the decay of the heavy dark state χ 2 → χ 1 e + e − via an off-shell dark photon. For the large splitting regime where additional decay channels open, we include them in the estimate for the χ 2 lifetime and in the corresponding branching ratio suppression as we focus on e + e − final states in this work. At the decay stage, the final state is directly analysed and the various event selection cuts described in the next sections are applied to the Monte Carlo truth. When considering up-scattering and down-scattering signatures, we use the built-in routines from BdNMC [10,97], adapting the cross-sections to the pseudo-Dirac DM case considered in this work as discussed above.

C. Backgrounds and selection cuts
There are three sources of background for the signal we consider. Cosmic rays deposits that are not associated to the cosmic ray tagger can mimic signal or interfere with reconstruction. We do not consider this background in our analysis and assume that cosmic rays signals can be efficiently removed. We also neglect interactions of beam neutrinos in the dirt surrounding the detector hall. The third background, which we do include in our analysis is interactions of beam neutrinos in the fiducial volume.
The dominant signature expected in the SBN experiments being a e + e − pair arising from a χ 2 decay, we focus first on the neutrino-induced background relevant for this final state. We then discuss the case of down-scattering and up-scattering before concluding with some remarks on the use of timing information.
a. χ 2 decay signature and angular separation The e + e − pair produced in decays of the heavy dark matter state χ 2 → χ 1 e + e − undergo bremsstrahlung in the detector, leading to a shower of e + and e − that are detected by the TPC. Provided that the separation between the original e + e − is sufficiently large, the individual four-momenta of these instigating particles can be reconstructed. This allows for detailed information about potential signal events to be measured and for separation from the backgrounds as we discuss below.
There are several ways in which neutrino interactions can mimic the signal of an e + e − appearing with no other activity. These include misidentified hadrons or muons, γγ of eγ with misidentified photons, and single γ conversion, along with no other reconstructed activity. We assume perfect identification of hadrons and muons in our analysis. Furthermore, we assume that γγ and eγ can be distinguished by the increased charge of one or both of the resulting showers and the displacement of the photon shower vertices from the neutrino interaction point by one radiation length. We focus on the photon conversion background.
To study this background, we apply the following event generation scheme. The flux of beam neutrinos is taken from Ref. [106]. The interaction of the beam neutrinos with argon is simulated using GENIE [107,108].
The unstable hadrons other than charged pions and neutrons from GENIE final state particles are decayed using Pythia8 [109]. Photons are injected into a very large volume of liquid argon in GEANT4 [103]. If the first interaction they undergo is conversion, which is the dominant process at the relevant energies, the resulting e − e + pair is added to the event. If they undergo Compton scattering instead, the resulting hard electron is added to the event. For any other type of interaction, which is an extremely rare occurrence, no new particles are added. For example, a tiny number of events have photons that undergo the photoelectric effect and those photons are ignored. The initial photon is then removed from the event.
We proceed to study the kinematics of this final set of particles. We do not apply energy or angular smearing to the particles. We veto events containing the following particles if they appear above the corresponding kinetic energy threshold: protons above 50 MeV, charged pions and muons above 30 MeV. This is based on Ref. [110], but we use a lower charged pion threshold as calorimetric reconstruction of charged pions will not be required here. Neutrons are not considered in our analysis. Events are required to have exactly two electrons above a kinetic energy of 30 MeV. We then examine the total number of background events from this procedure, as well as those that have an electron pair with an opening angle larger than 5 • and 10 • . The results of these selection criteria are shown in Table II. We further apply these kinematic criteria to the signal to determine a signal reconstruction efficiency. It is worth noting that electrons at an opening angle of even 5 • may not be easily reconstructed as a pair of particles. Since we apply these criteria to both the signal and background events, it is assumed that even if distinct particles are not reconstructed, an analysis can be developed to find events that look like a merged pair of showers. Since we apply our criteria to both signal and background, we should have a good estimate of the total expected number of such events.
In our analysis, we rely only on kinematic information about the particles produced in the signal heavy state decays and the background neutrino scattering in argon. Additional information may be used in an experimental analysis that could help to further reduce the background. Neutrino scattering can lead to production of a recoiling neutron. Neutrons are rather challenging to reconstruct as they tend to bounce around the detector, leaving small localised deposits whenever they hit a nucleus. The ability to reconstruct neutrons is still under study in LArTPCs and we ignore neutrons entirely in our analysis. The ability to veto neutrons, as well as short proton tracks (i.e. those from protons below 50 MeV), would help to reduced the background significantly. Furthermore, neutrino scattering produces a nuclear remnant that will decay into a stable nucleus. The de-excitation of the nuclear remnant leads to low energy photons and is challenging to detect, but is certainly a distinguishing feature of the background compared to the signal.
b. Down-scattering and up-scattering background We next discuss the case of up-and down-scattering signatures. For the small splitting regime with standard DM-electron scattering signature, we require the scattered electron to be extremely forward in order to mitigate the beam neutrino backgrounds (see e.g. [88]), and we also impose an energy cut similar to that for the decay signature: where θ is the angle of the outgoing electron with respect to the interaction point direction and E e its energy in the laboratory frame. This cut leaves the neutrino-electron scattering as the dominant source of background. In our projections we will consider the electron scattering signature for both ICARUS-OA and the SBN off-target run. Given the small neutrino-electron scattering rate, the offaxis location of ICARUS with respect to the NuMI beamline, and the suppressed neutrino flux in off-target run mode, one can expect a near background free search in these configurations after Eq. (22) is imposed [12,88]. Up-scattering off electrons with subsequent decay in the detector is a distinctive signal with few relevant backgrounds. If the decay is prompt, then the signal effectively is three electromagnetic showers with one being at narrow angle with respect to the beam. While this signal could be challenging to reconstruct if the showers overlap, there are few backgrounds from either neutrinos or cosmic rays that could fake it. The most plausible is charged current electron neutrino interactions with inelastic π 0 production and a short-distance photon conversion. Every piece of this is rare (electron neutrinos, inelastic pion production, and short-distance photon conversion), so the background is expected to be very small. It is possible for parameters where the decay is a bit displaced with respect to the up-scatter that this ν e CC + π 0 background could be more significant, but it is likely still too rare to pose a significant challenge. For an offtarget run, this potential background would be even further reduced. It is therefore plausible that such a search would be effectively background free. As discussed above, up-scattering off electrons is typically kinematically suppressed in the large splitting regime where χ 2 undergoes rapid decay (see Eq. (19)), and for this reason we will not investigate it further in this work. However, we note that it could still be relevant for lower DM masses and the higher energy DM produced in the NuMI beam, particularly given the distinctive characteristics of the signal.
Similar signals with up-scatter off nucleons face more potential backgrounds from neutrino scatter off nucleons that are less well studied. Neutrino-nucleon scattering with additional production of π 0 in particular are a potential background for χ 2 decay lengths in the tens of cm. For short χ 2 lifetimes, charged-current production of pe − could also pose a challenge if the electron shower is misreconstructed as two separate showers. As there are no backgrounds that exactly mimic the signal, there are many handles to enhance the signal-to-background ratio. Furthermore, an off-target run would further reduce possible neutrino backgrounds. Although a full study of reconstruction and background reduction in this scenario is beyond the scope of this work, background should not be a major impediment to this signal. c. Timing In some parts of parameter space, it could be advantageous to exploit timing information about events. Neutrinos travel at speeds very close to the speed of light, while heavier states may travel appreciably slower, leading to a lag relative to the backgrounds. LArTPC detectors have relatively coarse timing capabilities, with microsecond time resolution, leading to delayed signals only for very slowly moving states. The LArTPC detectors are, however, equipped with photomultiplier tubes (PMT) used to detect scintillation light associated with charged particles passing through the liquid argon that can yield nanosecond timing resolution. This could dramatically increase the power of timing to distinguish signal from background, but is as yet untested and so we do not exploit this possibility.
Timing information has been partially used by the Mi-crobooNE collaboration to study detection of heavy neutral leptons [43]. In their search, cuts were applied to look for events that fell outside of the entire batch window of 1.5 µs. This sort of timing analysis is somewhat simpler to develop, but is not suitable for more highly boosted particles and discards a significant number of signal events. On the other hand, it allows for drastically reduced backgrounds.

IV. RESULTS
In this section, we present our sensitivity projections of the SBN experiments to the iDM model. These projections are derived based on our simulations for iDM production, χ 2 decays, and χ 1,2 scattering, as well as our treatments of detector thresholds and backgrounds (in the case of the χ 2 decay signature), as discussed in Section III. The SBN experiments will be able to probe two qualitatively distinct parameter regions that predict the observed DM relic density and are currently unconstrained: 1) DM masses in the 100 MeV -1 GeV scale with relatively large mass splittings ∆ χ 0.2, and 2) DM masses in the 10 -100 MeV scale with relatively small mass splittings ∆ χ 0.1. The former is viable due to a sensitivity gap between beam dump and colliders searches that emerges for moderate χ 2 decay lengths, cτ χ2 ∼ O(mm − cm), while the latter corresponds to the opposite case of long-lived or quasi-stable χ 2 . We will discuss the prospects of SBN experiments in both regimes.

A. Large mass splitting region
When the mass splitting between χ 2 and χ 1 becomes large enough to lead to a relatively short-lived χ 2 , the sensitivities from beam dump experiments, such as CHARM and NuCal, is limited by the distance between their detector and the beam target. Moreover, for moderate values of ε, e + e − colliders such as BaBar have limited sensitivity due to smaller production rates and luminosity, or, in some cases, the lack of dedicated searches for displaced decays. This creates a sensitivity gap in the iDM parameter space corresponding to sub-ns χ 2 lifetimes (similar to the well-known "Mont's gap" in visibly decaying dark photon searches). Based on this discussion, it is clear that the SBN experiments, with baselines of order hundreds of meters, will not have sensitivity to primary beam-produced χ 2 decays in this large mass splitting parameter region. Additionally, we expect that for such large values of ∆ χ , χ 1 up-scattering off electrons will typically be kinematically suppressed, particularly for iDM produced with the Booster beam. On the other hand, the SBN experiments can instead efficiently search for χ 1 up-scattering off nucleons, χ 1 N → χ 2 N , followed by the fast semi-visible decay χ 2 → χ 1 e + e − .
We illustrate the potential reach of the SBN experiments in the large mass splitting regime in Fig. 2. Our projections assume the full anticipated datasets (see Table I): 6.6 · 10 20 PoT from the BNB for MicroBooNE and SBND and 7.7 · 10 21 PoT from NuMI beam for ICARUS-OA. For each experiment we display lines corresponding to 3 χ 1 -nucleon up-scattering events, assuming that the excited state χ 2 subsequently decays in the detector, leading to a prompt or displaced e + e − pair. In absence of dedicated background estimates for this particular signal channel, the results in Fig. 2 should be considered as qualitative possible targets for the SBN experiments in this range of the parameter space. It is worth emphasising that for the relevant coupling ranges ε ∼ 10 −3 , the typical decay length times boost factor of χ 2 is submetric, implying that most of the up-scattering events will be followed by a displaced decay in the detector, offering a clean way of further reducing the background. In view of this, we restrict our projections Fig. 2 to parameter regions corresponding to relatively short decay lengths, cτ χ2 smaller than around 10 m. We observe that all three SBN experiments can probe the thermal relic density target, which is indicated as a solid black line. In particular, for R χ = 2.5, we find that a significant part of the unconstrained parameter space predicting the observed relic abundance can be probed by the SBN experiments, corresponding to DM masses in the several hundred MeV range. Note that in the largest mass range, the state χ 2 can also decay to final states containing a pair of muons. While beyond our current scope, the di-muon final state represents another interesting experimental signature that could be leveraged by the SBN experiments.
Along with our SBN projections, Fig. 2 displays the existing constraints from DELPHI [111,112], CHARM [53,[78][79][80], NuCal [53,81,82], and BaBar [54,86,87], as discussed in Section II C. We have additionally derived a reinterpretation of the CHARM result including in particular the parton-level production not considered in [53,79], significantly extending the reach to the larger mass range. Along with these existing constraints, there are several ongoing/proposed experiments which can probe the iDM model. We will defer our discussion of these experiments to the end of this section.
Finally, we note that the large mass splitting region is particularly interesting since the short χ 2 decay length may open a region of parameter space that is compatible with both an explanation of the (g − 2) µ anomaly [113] and the relic density target, see e.g., Ref. [87]. While recasts of existing BaBar searches already probe the (g−2) µ interpretations shown in Fig. 2 [54,114], we find the SBN experiments can provide an independent dedicated tests of these important parameter regions. Furthermore, for very large splittings ∆ χ ∼ 1, it was recently pointed out in [56] that a fraction of the (g − 2) µ -favoured parameter space remains unconstrained by BaBar mono-photon searches. We illustrate in Fig. 3 the potential sensitivity of SBN experiments by showing the 3-event lines for DM-nucleon up-scattering followed by a χ 2 decay, with the parameter choice R χ = 4, α D = 0.1 and ∆ χ = 1. In the (g − 2) µ -favoured range of kinetic mixing, we expect around 10 4 − 10 5 such events depending on the experiment, significantly above any potential backgrounds 4

B. Small splitting region
We now turn to the second DM-relevant regime in which the mass splitting between the two states χ 2 and χ 1 is small, leading to long-lived or even quasi-stable χ 2 . We will first discuss the prospects for observing the semivisible decays of beam-produced χ 2 particles in the SBN detectors. As was noted in previous sections, the potential signals associated with such decays will have to compete with a mild neutrino-induced background. In order to obtain accurate sensitivity projections, we include the background estimates described in the previous section, corresponding to around 20k single photon We have estimated the CHARM sensitivity (grey purple region) at 95% C.L. (3 events) results based on the procedure described in the main text. The upper grey region is excluded by DELPHI [111,112], and BaBar [54,86]. The lower grey region is the NuCal exclusion from [53]. The green region is the 2σ range for the (g − 2)µ anomaly [113].
background events for both the SBND and ICARUS-OA full datasets. As discussed in Sec. III C this background can be strongly reduced by selecting events with a large opening angle. Indeed, we illustrate in Fig. 4 the expected opening angle distribution of the final e + /e − pair for both background and signal events for a DM mass of M χ1 = 150 MeV at SBND. Signal events generically feature much larger opening angles than the dominant neutrino-induced background. We further show in Fig. 4 that the signal production mechanism, either meson decay or proton bremsstrahlung, plays an important role in  [56,86]. The green region is the 2σ range for the (g − 2)µ anomaly [113].
fixing the shape of the signal distribution. The narrower opening angles predicted by proton bremsstrahlung production reflect the fact that the dark photons originating from this process are typically much more energetic than their counterparts produced through meson decays. We estimate the usefulness of an opening angle selection cut by plotting in Fig. 5 the efficiency of the cut for our signal sample, sig θ , along with the expected 95% C.L. limit on the number of signal events, N lim 95%C.L. , as a function of the opening angle cut, θ cut . It is clear that a 5 • cut can significantly reduce the background (as seen from the sharp decrease of the required N lim 95%C.L. ) at a moderate cost on the signal efficiency. Larger cuts do not lead to further improvement, but do not significantly reduce the signal to background ratio either. This indicates that the proposed search strategy may be viable even for experimental setups with only mild angular sensitivity. The main caveat is that the experimentally accessible mass range will be decreased for larger angular cut, which follows from Eq. (16). We will thus present our sensitivity estimates for long-lived χ 2 semi-visible decays using a θ cut = 5 • cut.
Our results for long-lived χ 2 decay signatures are presented in Fig. 6, which shows the projected reaches for MicroBooNE, SBND, and ICARUS-OA based on their full expected datasets (see Table I).
As discussed in the previous paragraph, we have imposed a cut on the e + /e − opening angle, which helps to mitigate the neutrino-induced backgrounds. We observe that the SBN program can probe the relic density target (solid black line in Fig. 6) for DM masses in the  MeV range, where existing experimental constraints are not competitive. More broadly, the long-lived χ 2 decay signatures can test a significant range of currently unexplored parameter space for DM masses between 50 MeV and a few GeV.  [85] and BaBar [86,87]. The LSND-excluded region at ∆χ = 0.05 is extracted from [16] and recasted to Rχ = 2.3, and directly extracted from [53] for ∆χ = 0.05. The BEBC limit is extracted from [42,90] and we included the MiniBooNE exclusion [89]. The grey purple region is the CHARM bound the procedure from [53]. The orange dot-dashed line is the reach of ICARUS in off-axis configuration, assuming 7.7·10 21 PoT from NuMI. The green dashed and purple dotted line are respectively the reach for SBND and MicroBooNE, assuming 6.6·10 20 PoT. Solid black line is the relic density target.
As stressed above, although the angular cut on the electron/positron pair does not significantly modify the signal to background ratio for θ cut 5 • , it directly impacts the lower mass threshold, following Eq. (16). For instance, increasing θ cut to 10 • would move the lowest accessible mass from ∼ 55 MeV to ∼ 100 MeV in Fig. 6. The mass threshold effect also reduces the reach of the SBN experiments as ∆ χ is lowered. We show this dependence in Fig. 7, which displays the projected 95% C.L. limits at SBN as a function of the mass splitting for a fixed DM mass M χ1 = 150 MeV. Although rapidly falling for smaller splitting, the SBN program will be able to probe a significant portion of the cosmologically motivated parameter space relevant for iDM. Since the experimental sensitivity for lower masses is strongly limited by the kinematic threshold discussed in Sec. III A, a "low angular cut θ cut , large background" approach could be an interesting additional search strategy. In particular, this could potentially provide independent coverage of the parameter space probed by the re-interpretation of LSND results, as discussed earlier in Sec. II C. Note that we quote the limits derived in [53] for ∆ χ = 0.05 and from [16] in ∆ χ = 0.1. For small mass splittings and low DM masses in the few -tens of MeV range, the excited state χ 2 will be quasi-stable, such that the decay signature is no longer effective. One can instead leverage χ 1,2 -electron scattering signatures (both up-and down-scattering) to provide coverage of this lower DM mass region. Since the heavy state χ 2 is very long-lived, the signature closely mimics those of both the standard light DM elastic scattering and the familiar neutrino-electron scattering [16]. Following our discussion in Sec. III C, we require the scattered electron to be very forward with respect to the incoming χ 1,2 direction, Eq. (22), which is expected to substantially suppress the dominant neutrino-electron scattering background for NuMI produced neutrinos at the off-axis ICARUS detector. We therefore present a 3-events line for χ 1,2 -electron scattering in Fig. 6 for the ICARUS-OA setup, including the angular cuts on scattering signal given in Eq. (22).
Altogether, we find that the variety of distinctive signatures present in the iDM scenario -long-lived χ 2 decay, χ 1 -nucleon up-scattering followed by fast χ 2 decay, and 'standard' χ 1,2 up-and down-scattering off elec-  Fig. 6). The dotted green line represent the 3 χe → χe events line for SBND in off-target configuration. The grey region are excluded from NA64 [85] and BaBar [86,87]. The LSND-excluded region is extracted from [16], the BEBC one from [42,90] and we included the MiniBooNE limit [89]. The grey purple region is the CHARM exclusion following the procedure from [53]. Solid black line is the relic density target for MV = 3Mχ 1 , ∆χ = 0.05.
trons -probe complementary regions of parameter space and can significantly extend the reach beyond current experimental constraints. We finally complement our study of the SBN program reach in standard neutrino (or anti-neutrino) run mode by showing in Fig. 8 our projections for a potential off-target run configuration for both SBND and ICARUS, in which the proton beam is directed past the main beryllium target and onto the iron dump approximately 50 meters downstream. We display 3-events lines for both long-lived χ 2 decays and χ 1,2 -electron up-and down-scattering signatures. This off-target configuration was employed in a year-long run by the MiniBooNE collaboration in its dedicated light search [88,89]. Although we have only considered a 2.2·10 20 PoT dataset (analogous to the MiniBooNE off target run), the potential gain for the SBND collaboration is impressive. Indeed, such a setup could help dramatically in reducing the beam-related neutrino background rates due to the suppression of decay-in-flight neutrinos. Furthermore, the close proximity of the SBND experiment to the iron dump allows for a substantial increase in the angular acceptance of dark sector states compared to the standard Booster run mode.
Besides the SBN program, there are number of existing and proposed experiments on the horizon that can probe the iDM model. In the near future, Belle-II [115] will significantly improve on the BaBar limits, particularly with dedicated searches for displaced decays [54,56,114,116].
In the lower mass regime, NA64 will benefit from a substantial increase in statistics (5 · 10 12 electrons-on-target or more) in the coming years [117], while next generation missing-energy/momentum experiments, such as LDMX [118] or positron-based experiments [119] have the potential of extending NA64 limits by more than an order of magnitude. Along with the SBN program, other proton beam fixed target experiments, such as Dark-Quest/LongQuest [20,53], JSNS 2 [19], and SHiP [120] will have substantial sensitivity to iDM candidates in the sub-GeV range. Finally, a rich experimental program is under development for long-lived particle detectors around the LHC interaction points [32,[121][122][123][124][125][126][127][128][129][130][131]. In tandem with the main general purpose detectors, the LHC can provide powerful and complementary sensitivity to light iDM [47,51,55].

V. CONCLUSIONS AND OUTLOOK
Proton beam fixed-target experiments, including accelerator neutrino beam facilities, offer a powerful means to search for light dark sectors. In this work we have examined the prospects of the three Fermilab SBN LArTPC detectors, MicroBooNE, SBND, ICARUS, to test models of light inelastic dark matter coupled to the SM through a kinetically mixed dark photon. Substantial fluxes of iDM states can be produced with both the on-axis 8 GeV Booster proton beam as well as the higher energy 120 GeV NuMi proton beam. Production from NuMI protons is particularly important for the large ICARUS detector, which is situated six degrees off-axis from the NuMI beamline. On the other hand, SBND, which is the detector closest to the BNB target, provides the best sensitivity to production from the Booster. Considering several experimental beam-target-detector run scenarios, we find that new substantial regions of iDM parameter space, including some predicting the observed thermal DM relic abundance, can be probed in the near future by all three experiments. The simple iDM model studied here leads to a rich phenomenology at the SBN experiments, furnishing several distinctive signatures including long-lived semivisible decays of the excited dark state, DM-nucleon upscattering plus fast de-excitation of the excited state, and up-and down-scattering with electrons. These signatures arise in qualitatively different regions of parameter space and thus offer complementary probes of the iDM model.
The iDM signatures must be disentangled from the neutrino interactions the experiments were designed to detect. For the e + e − pair signature of the long-lived decay of the excited state, we have simulated the dominant beam-related neutrino backgrounds for each of the experimental scenarios. Our findings suggest that a simple cut on the electron-positron pair opening angle can significantly mitigate the main background for this signal, which arises from conversion of secondary photons produced in neutrino interactions. Searches for long-lived decays of the excited state can probe new regions of parameter space, particularly for small mass splittings and DM masses in the 100 MeV range.
For small mass splittings and low DM masses of order 10 MeV, where the excited state is effectively stable, up-and down-scattering leading to a very forward single electron can provide a relatively clean signal of the iDM model. On the other hand, for large splittings and larger DM masses of order 1 GeV, DM-nucleon up-scattering followed by prompt/displaced decays of the excited state in the SBN detectors offers a potentially striking probe. Though a full study of the potential neutrino-induced backgrounds is needed to accurately assess the SBN sensitivities, the large predicted signal event rates make this channel quite promising and worthy of further study by the collaborations.
Within the proposed on-target run of the two beams, we highlight that further development of experimental techniques, such as tagging of soft baryons and exploitation of detailed timing information could help to further reduce the backgrounds considered in this work. Such possibilities merit further study as the understanding of LArTPC detectors improves.
Although we focused primarily on the near-term prospects by estimating the experimental reach during the planned 3-year SBN run, we additionally explored the potential of a supplementary one-year run in an "off-target" configuration, in which the BNB proton beam is steered past the target onto the iron dump. While all experiments would benefit from the anticipated reduction of beam-related neutrino-backgrounds, the SBND detector was found to profit most from this setup, with almost an order of magnitude improvements in its potential reach. In this light, it would appear especially warranted to explore the technical and logistical feasibility, as well as other potential physics motivations, of such an off-target run.
The full SBN program is coming online and will reach maturity on the several-year time scale. Our work reinforces the general expectation that these experiments have an excellent near-term opportunity to search for light new physics beyond the Standard Model, including light inelastic dark matter.