Millicharged particles in neutrino experiments

We set constraints and future sensitivity projections on millicharged particles (MCPs) based on electron scattering data in numerous neutrino experiments, starting with MiniBooNE and the Liquid Scintillator Neutrino Detector (LSND). Both experiments are found to provide new (and leading) constraints in certain MCP mass windows: 5 - 35 MeV for LSND and 100 - 180 MeV for MiniBooNE. Furthermore, we provide projections for the ongoing Fermilab SBN program, the Deep Underground Neutrino Experiment (DUNE), and the proposed Search for Hidden Particles (SHiP) experiment. In the SBN program, SBND and MicroBooNE have the capacity to provide the leading bounds in the 100 - 300 MeV mass regime. DUNE and SHiP are capable of probing parameter space for MCP masses in the range of 5 MeV - 5 GeV that is significantly beyond the reach of existing bounds, including those from collider searches and, in the case of DUNE, the SLAC mQ experiment.


I. INTRODUCTION
Extensions of the Standard Model (SM) by weakly interacting particles, and their probes at the intensity frontier experiments have become an important direction of particle physics [1]. One of the simplest and most natural ways of coupling new particles to the SM is via a "kinetic mixing" or "hypercharge portal" [2,3], which at low energies may lead to millicharged particles (MCPs), that would seemingly contradict the observed quantization of electric charge in nature [4]. In recent years, a wide class of related models were studied in connection with dark matter [5][6][7] (see also [8][9][10][11][12][13][14][15][16][17]), and MCPs can be viewed as a specific limit of those theories.
It is well appreciated that both proton and electron beam dump experiments provide sensitive probes of vector portal models. In particular, production and scattering of light dark matter [9] has been studied as a function of mediator mass m A , dark sector coupling α D , dark matter mass m χ , and kinetic mixing parameter Y . Depending on the relation between these parameters, either the past electron beam dump facilities [12] or the proton fixed-target experiments with a primary goal of neutrino physics [10,13,18] provide the best sensitivity. The simplest limit of m A → 0, when the parameter space simplifies to the mass and effective charge of MCPs, {m χ , }, was analyzed only in the context of electron fixed-target experiments [19][20][21], despite earlier studies of proton fixed-target experiments' potential as a probe of MCPs [22][23][24]. Clearly, fixed target neutrino experiments, such as the existing data from MiniBooNE [25] and the Liquid Scintillator Neutrino Detector (LSND) [26], and the soon to be released data from MicroBooNE, the ongoing SBN program [27], the Deep Underground Neutrino Experiment (DUNE) [28], and the proposed Search for Hidden Particles (SHiP) [29] serve as a fertile testing ground of MeV-GeV physics due to their high statistics [10,13,[30][31][32][33]. These experiments all serve as promising avenues to probe MCPs.
The purpose of this work is twofold: First, we demonstrate that existing data from LSND provides leading bounds on MCPs (surpassing existing constraints from SLAC's mQ experiment [19]) in the low mass regime (m χ 35 MeV). Likewise, newly released data from MiniBooNE [25] represents the current leading bounds on MCPs in the mass range of 100 MeV m χ 180 MeV (offering better sensitivity than collider experiments [21,34]). Second, we predict that by optimizing search strategies at ongoing and upcoming experiments (such as MicroBooNE, SBND, DUNE, and SHiP), fixed source neutrino experiments can serve to provide leading bounds for MCPs over the full range of masses 5 MeV m χ 5 GeV. The detection signature of MCPs in these experiments is elastic scattering with electrons, and we find that detection prospects are highly sensitive to the threshold imposed on the electron's recoil energy. Therefore, significant gains in sensitivity to MCPs may be achieved by future experiments by optimizing the detection of low-energy electrons. Our setup can be easily extended to probe light dark matter coupling to a massive dark photon.
Our results have direct implications for models with late kinetic coupling of dark matter and baryons [38] that could lead to extra cooling of the baryon fluid and spin temperature at redshifts z ∼ 20, which may result in a more pronounced 21 cm absorption signal. If a fraction of dark matter is composed of MCPs, this extra cooling mechanism can be naturally realized [39,40], and fit the unexpected signal reported by Experiment to Detect the Global Epoch of Reionization Signature (EDGES) [41].  Table I. The bound on N eff [35] comes from changing the effective number of neutrinos, while the SLAC bound (below ∼ 100 MeV) and the collider bound (extends up to ∼ 30 GeV in the plot) are taken from [19] and [21,34] respectively. The projected sensitivities at milliQan are from [34,36]. Our exclusions apply independently of the existence of a dark photon, which introduces additional constraints [37].
Our analysis reveals that sensitivities from LSND, SBND, SHiP, and DUNE can explore previously unprobed regions of parameter space that are favored by the 1%-MCP fractional dark matter hypothesis [37,40,42], as well as, a proposal to use MCPs to reduce electron number (while maintaining charge neutrality) to achieve an earlier decoupling of the baryon gas from the CMB [43].

II. PRODUCTION AND DETECTION
Fixed-target neutrino experiments rely on the production of neutrinos from weak decays of charged pions. In generating a large flux of π ± these experiments necessarily also produce a similar number [i.e. O(10 20 )] of π 0 [16]. For large beam energies, other neutral mesons (e.g. η, Υ, J/ψ) are also produced. Significant branching ratios to lepton pairs necessarily imply associated decays to pairs of MCPs, resulting in a significant flux of MCPs even for extremely small charges. In the case of η and π 0 , Dalitz decays π 0 , η → γχχ dominate, while for J/ψ and Υ direct decays, J/ψ, Υ → χχ are the most important. The branching ratio for a meson, M, to MCPs is given roughly by where M is the mass of the parent meson, X denotes any additional particles, and f (m χ /M ) is a phase space factor that decreases slowly as a function of m χ /M . For the 2-body decay cases, f (m χ /M ) has an simple analytic The number of MCPs passing through the detector is a function of both the branching ratio and geometric losses which varies substantially between experiments (see Table I).
We now turn to detecting MCPs in neutrino experiments, where the predominant signature is elastic scattering with electrons. Electron scattering dominates the detection signal because of the low-Q 2 sensitivity of the scattering cross section. Explicitly, in the limit of small electron mass, we have Integrating over momentum transfers, we see that the cross section is dominated by the small-Q 2 contribution to the integral. In this limit, we have dσ eχ /dQ 2 ≈ 4πα 2 2 /Q 4 , and therefore σ eχ ≈ 4πα 2 2 /Q 2 min . In the lab frame Q min is related to the recoil energy of the electron via Q 2 = 2m e (E e − m e ) [44]. An experiment's recoil energy threshold, E (min) e , then sets the scale of the detection cross section as Consequently, sensitivity to MCPs can be enhanced by measuring low electron energy recoils (an important feature for future search strategies).

III. RESULTS
The various curves in Fig. 1 are obtained by performing a sensitivity analysis [45]: given a number of predicted background events b and data n, the number of signal events s up consistent with the observation and backgrounds at (1−α) credibility level is found by solving the equation is the upper incomplete gamma function [46]. Using a credibility interval of 1 − α = 95% we calculate the corresponding bounds implied by s up according to the formula Equations (1) and (3) [47,48]). For the SHiP and DUNE experiments, we also include J/ψ and Υ mesons as well as Drell-Yan production which are discussed in the text. We use an efficiency of E = 0.2 for Cherenkov detectors, E = 0.5 for nuclear emulsion detectors, and E = 0.8 for liquid argon time projection chambers. The data at LSND and MiniBooNE is taken from [49] and [25,47] respectively. Projections at MicroBooNE [50], SBND [27], DUNE [28] and SHiP [51] are based on expected detector performance.
by particles inside the detector. Electromagnetic decays of mesons dominate the flux for m χ m η /2, whereas Drell-Yan production (DYP) dominates for large MCP masses (only accessible at DUNE and SHiP).
To estimate how many MCPs of energy E i arrive at the detector, we model the angular and energy distributions of the mesons using empirical formulas (discussed below) that parameterize the double-differential distribution of mesons (e.g. d 2 N π /dΩdE π ) in the lab frame. Then, for each meson produced at a given angle and energy, we calculate the differential branching ratio with respoect to the energy and angle of the MCP in the lab frame. This generates a double-differential distribution of MCPs in the lab frame for each meson. Next, using this new MCP-distribution, we determine the fraction of MCPs which pass through the detector with an energy within a bin centered at E i giving us a histogram-spectrum of MCPs as a function of E χ . Repeating this procedure over all production energies and angles of the meson, and appropriately weighting the contribution by the mesondistribution, yields each meson's contribution to N χ (E i ). DYP of MCPs from a quark and anti-quark pair is treated similarly, but to generate the MCP-distribution we integrate over the full production phase-space using MSTW parton distribution functions [52], and using Heaviside functions, we select the proportion of events containing an MCP pointed towards the detector, with energy E i .
Our main qualitative (but not quantitative) results can be understood without appealing to the details above, by introducing the quantity A geo (m χ ) which is defined as the ratio between the total number of MCPs that reach the detector and the total number of parent-mesons produced [53]. Given the number of parent mesons at a given neutrino experiment and A geo , then gives a rough estimate of the total number of MCPs that pass through the detector (e.g. at LSND for m χ = 1 MeV we have A geo ≈ 2% and N π 0 = 1.3 × 10 22 and so N tot χ ≈ i N (E i ) ≈ 2.6 × 10 20 . Using the information in Table I and Eq. (4) one can approximately reproduce our results.
At LSND, the π 0 spectrum is modeled using a Burman-Smith distribution [54,55] assuming 2 years of operation on a water target and 3 years of operation on a tungsten target. Our LSND analysis is based on [49], which featured 1.7 × 10 23 protons on target (POT), a beam energy of 0.798 GeV, and a single electron background of approximately 300 events with energies ranging between 18 MeV and 52 MeV. We estimate the N e /Area in Eq. (4) to be 2.5 × 10 26 e − /cm 2 .
The meson spectrum from Fermilab's Booster Neutrino Beam (BNB) is relevant for MiniBooNE, MicroBooNE, and SBND. The BNB delivers 8.9 GeV POT and so can produce substantial numbers of both π 0 and η mesons. The former's angular and energy spectra are modeled by the Sanford-Wang distribution [16,56], and η mesons by the Feynman Scaling hypothesis [56]. These distributions are common across all three experiments. Our geometric acceptances are in reasonable O(1) agreement with [16].
At MiniBooNE we perform two analyses: First, we consider the recently updated neutrino oscillation search [25]. We combine data from both neutrino and antineutrino runs and consider a sample of 2.41 × 10 21 POT with a single electron background of 2.0 × 10 3 events and a measured signal rate of 2.4 × 10 3 . Next, we consider a parallel analysis involving electron-recoil data with 1.86×10 20 POT [47]. Backgrounds were suppressed by operating the beamline in an "off-target" (i.e. not collimating charged pions) beam-dump mode, and these are further suppressed (to 0.4 predicted events) by imposing a cut of cos θ > 0.99 on the electron's recoil angle [47,48]. In both cases we estimate an electron number density of 3.2 × 10 26 e − /cm 2 and do not include timing cuts in our analysis.
At MicroBooNE, the meson rates assume 1.32 × 10 21 POT and we estimate an electron density of 3.9 × 10 26 e − /cm 2 . The chosen recoil cuts are based on the lowest reaches achievable given the wire spacing in Mi-croBooNE's liquid argon detector [50]. The wire spacing is 3 mm and the ionization stopping power is approximately 2.5 MeV/cm, so electrons with total energy larger than at least 2 MeV produce tracks long enough to be reconstructed across two wires. Requiring that ionization signals don't shower limits our recoil window to be between 2 MeV and 40 MeV. The treatment of SBND is similar, except we assume 6.6 × 10 20 POT, which corresponds to half the run time of MicroBooNE.
At SHiP we assume 2 × 10 20 POT and a near neutrino detector 50 m from the beam stop with an electron den-sity of 2.7 × 10 26 e − /cm 2 . We include J/ψ and Υ, in addition to π 0 and η. We do not include mesons such as ρ, ω and φ, because they do not serve to significantly alter the sensitivity offered by J/ψ. At SHiP's energies, production of π 0 and η is described by the BMPT distributions [16,57]. We have compared our geometric acceptances to those obtained using [16] and found reasonable agreement, with our acceptances being smaller by a factor of four. For production of J/ψ, we use the energy production spectra described in [58]. These distributions rely on production being highly peaked in the forward direction and parameterized as dσ/dx F ∝ (1 − |x F |) 5 , where x F = 2p / √ s is the meson's longitudinal component in the COM frame of the collision. We account for geometric losses by using an empirical formula for the p T distribution provided in [59]. We assume that the production spectrum of Υ mesons are similarly given, and normalize their total cross section to the data in [60]. For J/Ψ we have reproduced the Pb rates in Table 3 of [61], while for Υ we reproduced the Pt rates in Table 1 of [62]. We estimate N J/ψ = 2.1 × 10 15 with an acceptance of A geo (100 MeV) = 8 × 10 −2 , and N Υ = 1.2 × 10 11 with A geo (100 MeV) = 7.2 × 10 −2 . For large MCP masses, DYP becomes the main production mechanism, and we calibrate our DYP calculations by reproducing the dimuon invariant mass spectrum in Fig. 11 of [63] from the FNAL-772 experiment [64].
At DUNE, our treatment of meson production is very similar to the treatment at SHiP. We model pseudoscalar meson production using the BMPT distribution, but use a beam energy of 80 GeV [28] and account for differences in the target material. We also include J/ψ and Υ mesons and treat them as described above. Our detector treatment and electron recoil cuts are motivated by MicroBooNE's liquid argon time projection chamber (LAr-TPC) detector, and its ability to measure lowenergy electron recoils. We assume 3 × 10 22 POT and a 30 tonne liquid argon detector which corresponds to 5.4 × 10 25 e − /cm 2 . We estimate N J/ψ = 3 × 10 16 with an acceptance of A geo (100 MeV) = 2.4 × 10 −3 and N Υ = 5.1 × 10 9 with A geo (100 MeV) = 3.7 × 10 −3 . It is important to point out that the attenuation of MCPs through dirt and rock only becomes important for MCPs lighter than ∼ MeV [20] with large and is neglected in our analysis. This could have small effects on the MCP fluxes in all mass ranges and deserves a seperate and dedicated study. We also neglect the multiple scatterings of MCPs inside the detectors, which could be used as an additional tool of discriminating their signature against the neutrino background.
We consider two classes of backgrounds for Table I: those coming from each experiments flux of neutrinos [i.e. νe → νe and νn → ep], and those coming from sources such as cosmics, mis-identified particles, or dirt related events.
We treat neutrino induced backgrounds by summing over the neutrino fluxes from each collaboration and accounting for the detection efficiencies E. A large background reduction is obtained by imposing the electron recoil cuts E (max) e shown in Table I. These do not significantly affect the signal (which is dominated by low-energy electron recoils), but do significantly reduce charged and neutral current backgrounds [65,66]. In the case of MiniBooNE's dark matter run based on electron recoils, the angular requirements already provide a powerful background discriminant and the maximum energy of the electron is determined by the kinematics of the event.
External sources of backgrounds are estimated by multiplying the neutrino-induced backgrounds by an overall multiplicative factor. LAr-TPC detectors can use timing information as vetoes to reduce additional sources of backgrounds; this is not possible in a nuclear emulsion chamber. Therefore, we multiply our neutrino induced backgrounds by a factor of 10 for LAr-TPC detectors (MicroBooNE, SBND, and DUNE) and a factor of 25 for nuclear emulsion detectors (SHiP); this increase in the backgrounds decreases our sensitivity to by 20 − 30%. We have chosen an overall multiplicative factor because this naturally scales the backgrounds with the overall size of each detector. We anticipate that the timing information in LAr-TPC detectors, coupled with their fantastic vertex resolution, will allow it to veto backgrounds more effectively than nuclear emulsion detectors whose tracks do not contain any timing information; we have consequently chosen the factors of 10 and 25 respectively. We emphasize that our results in Fig. 1 can be easily revised for different background assumptions according to [46].

IV. OUTLOOK
Fixed-target neutrino experiments can probe MCPs due to the large number of mesons produced with electromagnetic decay pathways. Using existing data, LSND and MiniBooNE are able to provide the leading sensitivity to MCPs for certain sub-GeV masses. Beyond serving as a test of fundamental physics, this newfound sensitivity has implications for models of physics beyond the Standard Model. In particular it further restricts the parameter space of cosmological models where a fraction of MCP dark matter results in extra cooling of baryons that modifies 21 cm physics at high redshifts.
Equally important are our projected sensitivities at MicroBooNE, SBND, DUNE, and SHiP. The success of these experiments as probes of MCPs will rely heavily on their respective collaborations' search strategies. Increasing the sensitivity to low-energy electron recoils can be enhance the signal with a scaling proportional to 1/(E e − m e ). MicroBooNE has shown preliminary work measuring electron recoils with kinetic energies as low as 300 keV [50]. If this can be achieved, it is conceivable that the combined sensitivity of LSND, SBND, MicroBooNE, DUNE and SHiP will provide the leading sensitivity to MCPs in the full range of 5 MeV m χ 5 GeV.
Further progress may come from new experimental concepts and innovations. Significant progress could come from coupling large underground neutrino detectors with purposely installed new accelerators [13,67]. Millicharged particles may also be searched by experiments in disappearance channels [68][69][70][71], where e + e − → γ + χ +χ and Z + e − → Z + e − + χ +χ production leads to anomalous missing momentum/energy from the χpair that pass through a detector without depositing energy. Because of the advantageous scaling with (second, rather than the fourth power), there are clear prospects on improving bounds on MCPs above the 100 MeV energy range.
In addition, it is also interesting to place a milliQantype detector [34,36] downstream of a proton-fixed target along the Fermilab NuMI/LBNF or CERN SPS beams. The high flux of MCPs produced on the fixed-target combined with the dedicated milliQan MCP detector could yield a much more improved sensitivity reach of MCPs with masses below 5 GeV. This idea is implemented in [72] and is named Fermilab MINIchargesearch (FerMINI) experiment.
V. ACKNOWLEDGMENT