Millicharged Particles from the Heavens: Single- and Multiple-Scattering Signatures

For nearly a century, studying cosmic-ray air showers has driven progress in our understanding of elementary particle physics. In this work, we revisit the production of millicharged particles in these atmospheric showers and provide new constraints for XENON1T and Super-Kamiokande and new sensitivity estimates of current and future detectors, such as JUNO. We discuss distinct search strategies, specifically studies of single-energy-deposition events, where one electron in the detector receives a relatively large energy transfer, as well as multiple-scattering events consisting of (at least) two relatively small energy depositions. We demonstrate that these atmospheric search strategies -- especially this new, multiple-scattering signature -- provide significant room for improvement in the next decade, in a way that is complementary to anthropogenic, beam-based searches for MeV-GeV millicharged particles. Finally, we also discuss the implementation of a Monte Carlo simulation for millicharged particle detection in large-volume neutrino detectors, such as IceCube.

For nearly a century, studying cosmic-ray air showers has driven progress in our understanding of elementary particle physics. In this work, we revisit the production of millicharged particles in these atmospheric showers and provide new constraints for XENON1T and Super-Kamiokande and new sensitivity estimates of current and future detectors, such as JUNO. We discuss distinct search strategies, specifically studies of single-energy-deposition events, where one electron in the detector receives a relatively large energy transfer, as well as multiple-scattering events consisting of (at least) two relatively small energy depositions. We demonstrate that these atmospheric search strategiesespecially this new, multiple-scattering signature -provide significant room for improvement in the next decade, in a way that is complementary to anthropogenic, beam-based searches for MeV-GeV millicharged particles. Finally, we also discuss the implementation of a Monte Carlo simulation for millicharged particle detection in large-volume neutrino detectors, such as IceCube. FIG. 1. Artistic rendition of this work. Millicharged particles, χ, are produced in cosmic-ray showers from neutral meson decay. They lose energy as they travel through the Earth, with only the highest energy ones reaching the detector preferentially from directions of small overburden. These produce electromagnetic signatures in underground neutrino and dark matter detectors. In this work, we study both single-and double-hit signatures and relate their sensitivities.

CONTENTS
proposing further searches that can be done with atmospheric MCP. Motivated by Ref. [21], we explore multiplehit signatures in which a given MCP traversing a detector can scatter off two or more electrons, leaving a faint track. This search is highly advantageous in detectors that can identify low-energy electrons, such as the upcoming multi-kiloton-scale unsegmented liquid-scintillator neutrino experiment JUNO [33,34]. Fig. 1 offers an artistic rendering of this work's main ideas. High-energy SM particles, like protons, bombard the atmosphere, producing rich particle showers. If MCP exist, then they can emerge in these showers and travel through the Earth. Large-volume detectors provide excellent targets for these MCP, which can scatter once (right track), potentially imparting enough energy on the target electron to be a strong signal in these detectors. If the MCP scatters multiple times (left track), the faint track it provides is difficult for background processes to mimic. Some flux of MCP could also travel through the Earth (bottom track), if its mean free path through the Earth is long enough, and come upwards through the detector.
The remainder of this work is organized as follows: in Section II we revisit previous simulations of atmospheric MCP production and discuss the tools we use for our simulation. Section III discusses MCP propagation through the Earth, including energy loss and possible absorption, leading to an attenuation of the upward going MCP flux passing through a detector. In Section IV we provide the details of our experimental simulations for both current and upcoming experiments and construct sensitivity estimates for both the single-scattering and multiple-scattering searches. Section V discusses some Monte Carlo techniques for searches in even larger detectors, such as IceCube and the upcoming IceCube Upgrade. Finally, in Section VI we offer some concluding remarks.

II. MILLICHARGED PARTICLE PRODUCTION
A careful consideration of the production of millicharged particles in the upper atmosphere is crucial to this proposed search strategy. We provide the details of this approach with respect to its formalism in Section II A. Several sources of uncertainty are relevant as well, and these uncertainties, unfortunately, can plague any search for millicharged particles that relies on their production from the decays of neutral pseudoscalar/vector mesons. We discuss these uncertainties in Section II B.
We make a minimal assumption regarding the nature of this new millicharged particle -that it has some mass and small coupling to the Standard Model photon via a small electric charge. It is possible that such an MeV/GeVscale particle constitutes some fraction of the Dark Matter in the universe. If this is the case, a number of additional constraints apply. This has been explored as a potential solution to the EDGES anomaly [35] in, for instance, Refs. [36][37][38]. Various effects of millicharged particles as dark matter, leading to stringent constraints in the MeV-GeV mass range, are discussed in Refs. [31,32,[39][40][41][42][43]. If a millicharged particle is discovered via the atmospheric-production and scattering-detection approach we propose, then it is at most a very small fraction of the relic dark matter in the universe. After such a discovery, a great deal of scrutiny must be applied to determine a consistent picture of this new particle with cosmological and astrophysical observations.

A. Formalism
Given that cosmic rays are composed mostly of protons, their collision with the upper layers of Earth's atmosphere mimics the setup encountered in a proton beam dump experiment, with nuclei in the air playing the role of the target. An extensive cascade of radiation, ionized particles, and hadrons is generated with energies ranging from a few GeV up to dozens of EeV [44]. Among the mesons produced in the cascade, it is possible to find pseudoscalar mesonssuch as π 0 , η -and vector mesons -such as ρ, ω, φ and J/Ψ. If a millicharged particle exists with a mass below half of any given meson mass, then it can be produced in two-or three-body decays, replacing the final-state electrons in relatively common processes such as π 0 → γe + e − and J/Ψ → e + e − .
All of these mesons are unstable and have very short lifetimes and although none of them reaches the surface of the Earth, they can decay via a photon-mediated process to millicharged particles. The MCP are assumed to be stable particles and can reach the surface of the Earth and propagate through it 1 in such a way that they could be detected in underground detectors such as neutrino and dark matter experiments.
In this work, we will adopt a minimal model were the MCP (χ) is described by a stable particle coupled to the photon with strength ε × e and mass m χ . Taking this into account, the production profile of millicharged particles generated in air-showers from a parent meson m, can be described by the cascade equation [45] dΦ where ρ(X) is the atmospheric density at column depth X, λ m = γ m β m cτ m is the decay length of the parent meson m, and dΦm dEm d cos θ is the production rate of the meson with energy E m at zenith angle θ. Here, dn dEχ is the energy distribution of the millicharged particle in the decay, which can be written in terms of the branching fraction and the decay rate distribution as The branching fraction Br(m → χ) and decay width Γ(m → χ) must be specified accordingly for the three-body decays of the pseudoscalar mesons and the two-body decay of the vector mesons. The expressions for these quantities, as well as the kinematic integration limits in Eq. (1), are specified in Appendix A. The production rate of the mesons in the atmospheric cascade can be solved numerically. We have used the Matrix Cascade Equation (MCEq) software package [46,47], which includes several models for the cosmic ray spectrum, hadronic interactions, and atmospheric density profiles. In this work, for our benchmark results we have used the SYBILL-2.3 hadronic interaction model [48], the Hillas-Gaisser cosmic-ray model H3a [49], and the NRLMSISE-00 atmospheric model [50] to obtain the production rate of mesons. Given that millicharged particles are stable, a direct integration of Eq. (1) in the column depth X can be performed to find the differential energy spectra for a fixed zenith angle θ. In Fig. 2, we show the energy and angular dependence of the MCP flux at the surface of the Earth, for each one of the parent mesons considered in this work.

B. Uncertainties
The uncertainties in the production of MCPs are mostly due to the cosmic-ray model (CRM) chosen to generate the primary spectrum, as well as the hadronic interaction model (HIM) used to simulate the production rate of mesons in the atmospheric shower, the latter providing the dominant source of uncertainties, which grow as the energy increases. Ref. [51] recently explored these uncertainties in a fashion similar to how we have done here.  TABLE I. Comparison of meson production rate using different hadronic interaction models. Numerical values of ∆, the relative integrated production rate of a given meson, for the various cosmic ray models (top three rows) and hadronic interaction models (bottom three rows) considered in this work.
The composition and energy spectra of the primary cosmic radiation are characterized by a CRM that describes the primary spectrum with a power law that is ultimately fitted to air-shower data. The power law spectrum may break or not depending on the origin of the cosmic ray and the energy range observed. Additionally, the spectrum is typically characterized by the steepening that occurs for proton energies between 10 6 and 10 7 GeV, the so-called "knee," and an extra feature around 10 9 GeV called the "ankle." The origin of these features remains unclear, and it is an active research topic (see, for instance, Refs. [52,53]).
As mentioned before, the benchmark CRM we use to compute the MCP flux is the H3a model, which is widely used for calculations of atmospheric lepton fluxes [54,55]. However, there are other realistic CRMs which could be considered and that would yield to a more optimistic/pessimistic estimate of the MCP flux. To illustrate this point, we have considered the production rate of the π 0 meson in other CRMs, such as the Thunman-Ingelman-Gondolo model [45], the Gaisser-Stanev-Tilav model [56], and the poly-gonato model [57]. 2 The comparison of the differential production rate for these other CRMs as well as the ratio with respect to our benchmark model H3a as a function of the energy, can be seen in the right panel of Fig. 3.
On the other hand, to model the interactions of a primary cosmic ray with the atmosphere, we need a suitable model of hadron-hadron, hadron-nucleus, and nucleus-nucleus collisions. Yet, hadronic collisions at very high energies involve the production of particles with low transverse momenta, where present theoretical tools such as QCD are not enough to understand this feature. To address this problem, phenomenological models with Monte Carlo implementations are used. Since a hadronic interaction model provides the interaction coefficients of the coupled cascade equation used to describe the production rate of mesons, we can directly prove its impact by estimating the production rate, for a fixed CRM. In this work, we choose SYBILL-2.3c as benchmark model, and the QGSJET-II-02 [58], DPMJET-III [59], and EPOS-LHC [60] HIMs for comparison. The left panel of Fig. 3 displays the production rate of π 0 for all of these models as well as the ratio to our benchmark HIM.
To estimate the impact of the benchmark models used in this work, we evaluate the ratio of total production with respect to a different CRM/HIM, defined by where m is the meson of interest, j the index used to denote the model that is being compared with the benchmark model BM, E min the minimum energy available in MCEq, and Λ = 10 3 GeV is the upper energy cut that we use in order to obtain the MCP flux from a parent meson. Table I shows the ∆ coefficients for the different CRMs and HIMs for π 0 , η, and φ. A similar result can be found for the other mesons, except for J/Ψ, in which case there are not enough statistics on the meson production rate to evaluate the uncertainty properly. As can be seen from this result, the biggest source of uncertainty comes from the hadronic interaction models which would induce a difference in the MCP flux from about ∼ 16% up to ∼ 68%, depending on the parent meson. Even though all of the event generators considered here are updated with LHC data, the cross-section measurements for hadron production have rather large uncertainties [61]. On the other hand, we must also take into account the differences encountered in the phenomenological models that arise from the treatment of inelastic hadronic collisions within the framework of Reggeon Field Theory [62]. We stress that HIM uncertainties are present in most if not all modern MCP searches, as any beam-based 3 search needs to simulate hadronic interactions at some level. The bottom panel depicts the production ratio relative to SIBYLL.2.3.c, which we use for our production estimates. Right: Comparison of production rates for different cosmic ray models: TIG (orange), GST (pink), H3a (blue), and PG (green). The bottom panel shows the ratio relative to H3a, which we use for our production estimates. See text for further discussion of CRMs and HIMs.

III. PROPAGATION THROUGH EARTH AND ENERGY LOSS
As we will show in Section IV, underground neutrino detectors provide an excellent opportunity to search for MCPs produced in cosmic-ray air showers. However, the detection of MCPs requires detailed knowledge about charged particle propagation in the medium, since these experiments uses the Earth's crust to shield from atmospheric muons, which usually constitutes the main source of background. Because of this, it is expected that the MCPs that arrive at the detector lose part of the energy they had when reaching the Earth's surface. Just as with any other charged particle, the MCP would lose energy by ionizing the medium through which they propagate and by interacting with nuclei. The average energy loss along the MCP trajectory can be parametrized by where dX is the column density traversed, ε is the MCP coupling, and a x and b x are various categories of energy loss with (potentially) different scaling with the MCP energy E and ε. We simplify this expression, adopting the righthand side of Eq. (4) in our simulations, where a and b are the energy loss parameters given in units of [GeV/mwe] and [mwe] respectively ("mwe" being "meters of water equivalent"). On the other hand, we can estimate the overburden length traversed by an MCP approaching a detector located at depth d and along a trajectory with angle θ by with R ⊕ being the Earth's radius.
The probability that a millicharged particle with energy at the surface E i arrives at the detector with an energy E f will depend on the coupling α em ε and the incident direction cos θ. This probability is given by:  where the average distance R(E i , ε, E f ) can be obtained from Eq. (4), and is given by We have taken the energy loss parameters for a standard rock shielding with a = 0.223 GeV/mwe and b = 4.64 × 10 −4 mwe as reported in Ref. [64]. The left panel of Fig. 4 shows the probability that an MCP arrives at a detector at a depth of 1.5 km with E χ > 1 GeV for two different values of the coupling ε, with a variety of different initial energies as depicted in the label. The flux at the detector can be obtained by convolving the flux at the surface with the survival probability P. The right panel of Fig. 4 shows the angular distribution of the flux at detector for an MCP with a mass of 10 MeV arriving at the detector with energies of 1.0 and 5.0 GeV.
We note here that for relatively large ε 2 = 10 −2 , the upward-going flux, −1 ≤ cos θ ≤ 0 is highly attenuated. On the other hand, for smaller ε 2 = 10 −4 , because energy losses are suppressed by two orders of magnitude, the flux is mostly independent of cos θ. For experiments that can be sensitive to such small millicharges, the effective area of the sky that the detector can search will nearly double. Moreover, searches for upward-going events can assist in reduction of background events, for instance from atmospheric muons entering the detector from above.

IV. CURRENT AND FUTURE EXPERIMENTAL SEARCHES
To date, the most stringent searches for millicharged particles in the ∼MeV-GeV mass range have used antropogenic beams, where the MCP are produced in either the collision of two beams (collider searches) or when a beam impacts a target (beam-dump searches). Colliders place the strongest constraint on MCPs with masses above ∼1 GeV, primarily through a combination of direct searches at LEP and measurements of the invisible width of the Z boson [65]. The LHC can probe even larger, 30 GeV, masses [66]. Below 100 MeV, the strongest constraint is from the SLAC mQ electron fixed-target experiment [63]. Future dedicated experiments have been proposed in the context of both collider [22,25] 4 and fixed-target [24,26,27] environments, reaching sensitivity to ε 2 10 −6 and a wide range of masses, up to nearly 5 GeV.
Recently, searches for MCP in neutrino experiments have garnered attention. Sensitivities for accelerator production have been estimated in Refs. [20,21] and carried out by the ArgoNeuT collaboration in Ref. [29]. The latter has set some of the strongest constraints over a wide MCP range. Atmospheric production has been considered in a similar fashion to this work for the Super-Kamiokande detector and its future successor Hyper-Kamiokande [69] in Ref. [30], as have limits from particles produced in the interstellar medium and from Earth relics recently discussed in Refs. [31] and [32], respectively.
In the remainder of this section, we will discuss and derive constraints for MCP masses between 10 MeV and 1.5 GeV. Our results for current experiments, particularly Super-Kamiokande, demonstrate an improvement on current constraints for some masses, in agreement with the results of Ref. [30]. We divide the search strategy into two different types of analyses, those relying on single-hit signals (Section IV A) and those that rely on multiple energy depositions from a single MCP particle traversing the detector (Section IV B). The former is advantageous for largevolume, high-energy-threshold experiments (e.g., water Cherenkov detectors like Super-Kamiokande), whereas the latter is strongest in scintillator experiments with precision timing and low-energy thresholds (for instance, JUNO). We find that these multiple-scattering searches offer great potential for discovering millicharged particles in the 10 MeV to 1.5 GeV mass range and warrant additional focus in the coming decade.

A. Single-Scattering Searches
Via the same coupling ε that allows for the production of MCP in the upper atmosphere, the particles that reach a detector are capable of scattering (via a t-channel photon exchange) off detector materials. We will consider scattering off electrons for the remainder of this work. This massless-mediator scattering yields a differential cross section that peaks for small electron recoil energy, and so detectors with capabilities of identifying and reconstructing low-energy electrons will be advantageous in this endeavor. However, as we consider electrons of lower and lower energy, more and more backgrounds become relevant. In the following subsections, we will discuss the characteristics of the signal events, the various backgrounds, and the experimental limits we are able to derive in this single-hit analysis including the statistical techniques employed.

Signal
Once the flux of MCPs arrives at the detector, the millicharged particles can interact with the detector medium by scattering off electrons. We will be interested in low-energy recoils, as those are more frequent due to the shape of the differential cross section. The differential cross section between MCPs and electrons is given by [20] where E χ is the energy of the incoming MCP, E r is the recoil energy (assuming an initial stationary electron in the laboratory frame), α EM is the electromagnetic coupling constant, and m χ and m e are masses of the MCP and the electron, respectively. From the equation above, it is easy to verify that in the ultrarelativistic limit E χ E r , m χ , m e , the differential cross section scales like E −2 r , and therefore as the threshold of a given experiment becomes lower, the signal becomes much stronger.
The event rate from MCPs can be obtained by a convolution of the millicharged particle flux with the differential cross section multiplied by the number of targets and the detection efficiency. This is given by where n e is the total number of electrons in the detector and (E r ) is the (detector-dependent) reconstruction efficiency of these single-electron events. The event rate distribution for a 1 kt effective mass detector is shown in Fig. 5, as a function of the incoming direction from zenith angle and the recoil energy, for different values of ε 2 and an MCP mass of 10 MeV. Notice that for the angular distribution an integration in the recoil energy is needed, and the values in the label indicate different detection thresholds. The MCP event rate shown in this figure is a general result, which can be applied to any underground particle detector located at a depth of ∼ 1.5 km, while the mass dependence scales in a similar way as the MCP flux, with higher event rates at lower masses, where the millicharged particle receives contributions from most of the mesons.

Backgrounds
Here we discuss the different sources of background events that contribute to such single-scattering millicharged particle searches in water Cherenkov or liquid scintillator detectors. The expected background rates (summed over all contributions) in these two environments, along with a signal expectation, are shown in Fig. 6 for Super-Kamiokande phase II (left) and JUNO (right). More details on the Super-Kamiokande event distributions are given in Appendix B.
a. Penetrating muons: Our atmosphere is filled with muons generated from the decays of charged mesons that are dominated by pions at the lowest energies and kaons at higher energies [70]. These muons are highly boosted, and even though they have a short lifetime (∼ 2 × 10 −6 s), they can easily reach the surface, propagate through the Earth and leave a signal in the detector. The most efficient way to suppress this and other possible cosmogenic backgrounds is to have a sufficient amount of overburden over the location of the detector. This background may be rejected by measuring the opening angle of the Cherenkov cone of the candidate electron in the event [71].
b. Neutral-current neutrino events: Atmospheric neutrinos induce a source of background due to neutral-current interactions with nuclei in the medium. In the case of liquid scintillator detectors, the background is induced by interactions with 12 C whose de-excitation produces radiation. In Cherenkov detectors, elastic, neutral-current scattering of neutrinos off the target nuclei can lead to small energy depositions with a similar spectrum to the millicharged particle signature. However, this is a relatively small component of the background in Cherenkov detectors like Super-Kamiokande.
c. Charged-current neutrino events: Specifically in water or ice Cherenkov detectors (i.e., Super-Kamiokande, Hyper-Kamiokande, or IceCube), both ν e /ν e and ν µ /ν µ charged-current scattering events from the atmospheric neutrino flux may mimic our desired signal of a low-energy electron. The former directly creates low-energy electrons indistinguishable from the signal as long as the hadronic particles produced in the interaction are below the Cherenkov threshold. Despite the fact that they cannot be separated on an event-by-event basis, the recoil energy distribution of the electrons has a spectral shape different than the MCP signal, which allows them to be distinguished statistically. Similarly, ν µ events can mimic the background if the outgoing µ ± is below the Cherenkov threshold, but the Michel e ± coming from the µ ± decay are visible. As in the previous case, the Michel electron energy distribution is distinct and can be disentangled statistically; see Ref. [71] for further details.
d. Radioactive and anthropogenic backgrounds: Radioactive backgrounds are produced by the materials in and around the detector. These can produce electrons with energies similar to those struck by MCPs, e.g., in radioactive beta decays of nuclei, and depend on the experimental setup. For example, the dominant radioactive backgrounds in JUNO, in the energy range relevant for this analysis, are 8 B, 10 C, 11 C, and 11 Be, which we reproduce from Ref. [33] and include as background in the right right of Fig. 6. We note that this rate is significantly larger than currently allowed MCP signatures for E e 12 MeV and will be the limiting factor for single-scattering searches in detectors like JUNO. Searches for signals with an associated neutron in the final state, most prominently scattering of ν e from the diffuse supernova neutrino background, can reject these backgrounds by requiring coincidence with the signal of neutron capture. Since MCP signatures do not produce an outgoing neutron, but only the electron recoil, this method cannot be used to reduce the background. However, as we will discuss in Sec. IV B, multiple-hit signatures can be used to reduce radioactive backgrounds. The dominant anthropogenic backgrounds in the search of MCPs are produced by either nearby nuclear reactors or by accelerator facilities. In JUNO, the antineutrinos produced by nuclear reactors are the source used to study neutrino oscillations and can be separated from the MCP signature by the presence of the coincident neutron capture mentioned above. In Super-Kamiokande, accelerator neutrino events produced in the Tokai accelerator facility can be removed by searching for MCPs during off-beam times.

Experimental setups considered
In this section, we will briefly describe the experiments considered in this work and the constraints that can be established from the lack of significant MCP signal. We will derive our constraints by performing a forward-folding binned-likelihood analysis, where the data and expectations are organized in equally sized bins of electron recoil energy. To construct our test statistic, we compute, for each one of the experiments considered, the number of signal events expected in a given bin of electron recoil energy. The number of signal events in reconstructed electron recoil energy is given by where the symbols inside the inner most integral are defined in Eq. (9), T is the time exposure, E b is the bin size, E i is the bin center, and (E i ) is the average detection efficiency in a given bin.
To compute current experimental constraints, we construct a background-agnostic test statistic, comparing the observed data with that expected from an MCP with mass m χ and mixing ε, ignoring any potential background contribution. This procedure will always result in a relatively conservative upper limit on the parameter space of MCP and cannot allow for a potential preferred region. We adopt this strategy for the current Super-Kamiokanade and XENON1T experimental results for different reasons, which we will detail in their respective paragraphs to follow.
The background-agnostic test statistic is calculated using the bin-by-bin likelihood function [72] where d i and µ s i are the data and expected signal (given m χ and ε) in bin i, respectively, and P (d i |µ s i ) is the Poissonian likelihood of observing d i given expectation µ s i . This form of the likelihood is both background-agnostic and one-sided, guaranteeing the setting of a constraint instead of a preferred region. Our test statistic then is where the denominator is the signal-free likelihood function and will return 1 when using the background-agnostic, one-sided likelihood of Eq. (11). When deriving a constraint, we will assume that Wilks' theorem holds and set limits assuming we have two degrees of freedom, m χ and ε. In reality, the signal event distributions are all of the same shape (peaking at low electron recoil energy), where ε scales with the rate and m χ determines which mesons can contribute to the MCP flux, also impacting the overall normalization. This relation implies that these two parameters could, in principle, be viewed as a single degree of freedom. In light of this, our use of two degrees of freedom, combined with the above background-agnostic, one-sided likelihood function approach, should be viewed as conservative.
a. Super-K: Our analysis uses the event selection in Ref. [71], which was designed to search for diffuse supernova background neutrinos. This event selection was previously discussed in Ref. [30], where they estimate event rates of MCPs produced in the atmosphere. We use the data from SK-I, SK-II, and SK-III and perform a joint likelihood analysis combining the three phases. For our analysis, we include the signal efficiencies provided in Ref. [71]. The background expectations estimated by the collaboration can be extracted from the same reference. We demonstrate the expected event rate from the collaboration's background as well as the observed data during SK-II operation in Fig. 6 (left) as a function of electron recoil energy. However, it must be stated that this background expectation does not include the potential signal from diffuse supernova neutrino background (DSNB) events, the signature of interest in Ref. [71]. If we perform a likelihood analysis 5 comparing the observed data with the expected background (with no DSNB contribution) plus the MCP signal, we obtain a moderate (∼90% CL) preference for the existence of MCP; however, this signal is degenerate with the DSNB one. So, for this reason, we choose to adopt the conservative, background-agnostic approach discussed above to derive an upper limit on ε 2 as a function of m χ at 90% CL.
b. XENON1T: Dark matter direct-detection experiments can measure electron recoil energies down to a few keV. We use the recent XENON1T experiment [73] result to search for signatures of millicharged particles produced in the Earth's atmosphere. The background models presented in Ref. [73] are notably insufficient to explain the observed rate of electron recoil events. However, one potential background to explain the excess is the presence of an unconstrained tritium beta-decay background. Because of this possibility, we choose to adopt the background-agnostic approach as we did for Super-Kamiokande. If we calculate the likelihood using the given B 0 background model of Ref. [73], we observe a preferred region of parameter space at over 95% CL. This preferred region of parameter space is nearly excluded completely by other constraints at 2σ confidence and will be robustly tested by JUNO. See Appendix C for more details, including the preferred region of parameter space that we obtain. When using Eq. (9), we have assumed that all of the electrons in a detector are unbound -this is not the case in XENON1T, especially when comparing the observed recoil energy distribution with the nuclear binding energies of interest. In that context, our signal rates predicted in XENON1T will be overpredicted by a factor of several -see, for instance, Refs. [31,74] for further discussion of this effect. Given the scaling of our signal rate with ε 4 and the fact that XENON1T is not the most powerful experiment considered in this work, we disregard this effect in estimating our constraints.
c. JUNO: Given that JUNO offers a future search for single-scattering events, we calculate our test statistic assuming that the collected data in each bin is consistent with the background expectation and perform a comparison Single-hit constraints on millicharged particles. Constraints on millicharged particle parameter space from SLAC mQ [63], ArgoNeuT [29], the milliQan Demonstrator [68], and colliders [65], compared with our 90% CL single-hit constraints from Super-Kamiokande (purple), JUNO (orange), and XENON1T (grey).
with the expected signal-plus-background distributions. See the right panel of Fig. 6 for these backgrounds from Ref. [33] as well as a signal expectation. We then calculate our bin-by-bin likelihood using where µ b i (µ s i ) is the expected background (signal) in bin i. With this form of the binned likelihood, we then use the same expression for the test statistic T S in Eq. (12) above and again, assume two degrees of freedom when projecting JUNO sensitivity at 90% CL. For our analysis, we assume an exposure of 170 kt − years, which corresponds to approximately ten years of data collection. Fig. 7 shows our new constraints at 90% CL, compared with existing upper limits from the SLAC-mQ experiment [63], ArgoNeuT [29], the milliQan Demonstrator [68], and colliders [65] 6 . Additional constraints, potentially subject to large background uncertainties and systematics associated with hand-scanning, can also be found in Ref. [28]. Our Super-Kamiokande analysis, shown by the purple line and associated shaded region, yields results comparable results to those derived in Ref. [30] for the bulk of parameter space. The small discrepancies between these two results arise from the following differences. The analysis of Ref. [30] set constraints assuming an exposure of 22.5 kt − years at Super-Kamiokande, while here we have used the full 2853 day (or 7.81 × 22.5 kt − years) exposure of SK I-III and the statistical analysis discussed above to set this constraint. While the details of backgrounds are complicated, Hyper-Kamiokande should be able to improve on such constraints in the coming decades -we direct the reader to Ref. [30] for more on millicharged particle searches at Hyper-Kamiokande and Refs. [69,75] for a detailed discussion of searches for the DSNB at Hyper-Kamiokande and associated challenges.

Single-scattering constraints and sensitivity projections
As shown in Fig. 7, the strongest current limit for 100 MeV m χ 500 MeV comes from this Super-Kamiokande analysis. However, a similar duration of a JUNO single-hit analysis (orange line) will be able to improve on this, due to JUNO's ability to reach very low electron recoil energy. Although XENON1T (grey) can reach keV-scale electron recoils in this analysis, its small volume-exposure limits its capabilities relative to SK or JUNO.

B. Multiple Scattering Searches
Throughout Section IV A, we focused on scenarios where the atmospheric-produced MCP scatter inside a detector, providing enough energy to an electron to leave a distinct signature in the experiment. When discussing the signal of single-scattering events, we pointed out that the rate of signal events grows with smaller electron recoil energy E r . However, in the case of both Super-Kamiokande and JUNO, the background do as well. Especially for JUNO, the background rate (dominated by radioactive emissions) became prohibitive for E r 12 MeV, preventing a singlescattering analysis from reaching lower recoil energies.
One means for reducing (or even eliminating) these low-recoil-energy backgrounds is to require multiple scatterings within one small time window -the radioactive backgrounds are each associated with some relatively large half-life and so the possibility of two such emissions in a ∼10 −6 s window is very rare. Consequently, this type of analysis requires that a given MCP deposits energy at least twice as it traverses the detector, which will scale with even more powers of the millicharge ε than our single-scattering analysis. However, because this strategy allows us to search for even smaller recoil energies E r where the scattering cross section grows, we will see that large event rates are still possible. This principle was exploited in Ref. [21] for beam-produced millicharged particle searches in ArgoNeuT/DUNE.
In the remainder of this subsection, we will utilize this same approach for our atmospheric searches, focusing on the JUNO experiment. We will demonstrate that this approach can even outperform the single-scattering analysis presented in Section IV A due to the low-energy capabilities of the JUNO detector's liquid scintillator.

Signal characteristics
We will determine the rate of multiple-scattering signal events in a given experimental setup by means of comparison with the single-scattering analyses discussed in Section IV A. For simplicity, we will consider total event rates instead of the event distribution as a function of electron recoil energy E r . The important quantities to consider will be T H , the minimum recoil energy for a hard-scattering event (i.e., those considered in the single-scattering analyses) and T S , the minimum recoil energy capable of detecting events with multiple soft scatters in a small time window.
We perform the comparison between single-and multiple-scattering by computing the mean free path that a given MCP travels before either scattering in a "hard" manner (yielding an electron energetic enough to appear in singlescattering analyses) or in a "soft" manner (with electrons of too low an energy to be useful in single-hit searches, but energetic enough to be detected and used in a multiple-scattering search).
Assuming that a given MCP travels through a region with electron density n e and length L det. , we can approximate the probability that the MCP interacts using its hard-scattering mean free path λ H = 1/(n e σ H ), where σ H is the hard-scattering cross section, i.e., the cross section where we require the outgoing electron to be energetic enough to be distinct from radioactive backgrounds. In JUNO, this requirement is E r 12 MeV. Ref. [20] approximates σ H as The soft-scattering cross section and the corresponding mean free path between soft scatters are given by σ S and λ S , respectively. The cross section σ S is identical to Eq. (15) with T S replacing T H . In principle, T S can be significantly lower than T H , depending on the properties of the detector. For JUNO, we estimate that T S ≈ 0.01 − 1 MeV due to the 10 4 photons/MeV and O(10 3 ) photo-electrons/MeV of the detector's liquid scintillator [34,76]. The ratio of mean-free-paths λ H /λ S then is proportional to T H /T S .
In addition to the single-scattering probability P 1 in Eq. (14), we can also consider the probability that an MCP scatters softly n times inside the detector. If we divide the total detector length into "segments" using L det. = N seg. L seg. , this probability is Assuming we can take the N seg. → ∞ (or L seg. → 0) limit for an unsegmented detector, this probability approaches The probability of two or more hits, ∞ n=2 P (n) can be written as The ratio of multiple soft-scatter events to the single hard-scatter events is proportional to P n≥2 /P 1 . This quantity depends on ε 2 , the soft-and hard-scattering minimum recoil energies T S and T H , and the total detector length L det. . For two choices of T S and T H , we present the ratio of these probabilities as a function of L det. and ε 2 in Fig. 8. In each panel, above the solid black lines, the multiple-hit rate exceeds the single-hit one, implying that searches for these multiple-scattering events can be at least as powerful as the single-scattering searches. If JUNO is capable of searching for 0.1 MeV soft-scatters, this implies that these multiple scattering events can be favorable for ε 2 ≈ 10 −5 (assuming path lengths on the order of 10 m).
One final feature of the multiple-scattering signature that is not very useful in single-scattering searches is directionality. Because the single-scattering searches are seeking soft electrons, obtaining the direction of the incident millicharged particle is difficult. With a multiple-scattering analysis, we can obtain the angular distribution of these events, which should match the flux prediction, including attenuation through Earth, as discussed in Section III. This can be further used to statistically separate our signal from background events. time frame of an event for an experiment. The large background rates present at low recoil energy for JUNO ( 8 B, 10 C, 11 C, and 11 Be radioactive decays, specifically) will be suppressed by requiring multiple scatterings in this small window. For JUNO, the emission timescale of the liquid scintillator is at most ∼200 ns [33] -if we can restrict the time window to be O(10 −6 s), then these backgrounds are reduced to being O(1) for the ten years of data collection we assume for JUNO.

Projected sensitivity with multiple scattering
Combining the above information, we project the sensitivity using this multiple-scattering search in JUNO here in Fig. 9. In order to determine the "detector length" of JUNO that we discussed above, we simply average the path length of incoming MCP over the 35-m diameter spherical vessel. The average path length through a sphere is 2D/3 = 23.3 m for JUNO.
We assume that the backgrounds in JUNO are small enough that 10 signal events are statistically significant in the 170 kt-yr exposure -we do not use any spectral information about the energy of the electron events, as we expect that resolution at such sub-MeV energies will be difficult. We take three different assumptions about the minimum threshold energy of electrons that JUNO can detect -1 MeV (dot-dashed), 100 keV (dashed), and 10 keV (solid). If this 10-keV threshold is attainable (which could be possible given the 10 4 photons/MeV of energy deposited in the liquid scintillator), we see that this multiple-hit strategy will far exceed single-scattering searches for atmospheric MCP. Even with MeV thresholds, this is a complementary approach, especially at large m χ . The realistic 100-keV threshold seems particularly promising as a target for JUNO.

V. MILLICHARGED PARTICLES AT NEUTRINO TELESCOPES
The search for millicharged particles from natural sources is a statistically-limited problem, whose optimal detector would be a low-threshold, small-background, large-mass underground device. In this work, we have so far studied kiloton-mass-scale low-threshold and small-background detectors, such as XENON1T, and tens of kilotons mass lowthreshold detectors, such as JUNO, in this section we will study the sensitivity of megaton and gigaton detectors, known as neutrino telescopes. These detectors use natural transparent media to detect Cherenkov light produced by neutrino interactions and are designed to detect faint high-energy astrophysical neutrino sources. However, their large detector volume allows them to also have unique capabilities to observe low-energy neutrinos produced in supernovae.
Since MCPs are long-lived and have small energy losses their signature in neutrino telescopes will be a long, faint track. The main backgrounds for this search are coincident or dim atmospheric muons that can accidentally mimic the signature. Searches for long-lived faint-tracks have been performed by the IceCube collaboration by looking for an isotropic flux of fractional charged particles with charges motivated by the quark charges [77,78]. This analysis yielded a sensitivity approximately ten times stronger than previous constraints from Kamiokande and MACRO, however the assumption of an isotropic flux of fractional charge particles is not consistent with the propagation of these particles through the Earth, as we discussed in Sec. III. Despite this, Refs. [77,78] motivate further study of the sensitivity to these signals.
The yield of Cherenkov photons produced per energy per unit path length of the distance travelled by a particle of charge εe is given by the Frank-Tamm equation [17] where θ c is the Cherenkov angle in the medium, which we have normalized to the Cherenkov angle in ice, θ c ≈ 41 • [79]. The relevant Cherenkov photon wavelength for IceCube is around 400 nm, which results from the convolution of the IceCube PMT quantum efficiency and the glass-housing transmission probability [80,81]. At this wavelength the IceCube digital optical module (DOM) acceptance is approximately 0.15, and decreases by approximately half by reducing the wavelength to 350 nm or increasing it to 550 nm. Considering this, the yield of relevant Cherenkov photons is approximately 50ε 2 cm −1 . This implies that an IceCube-through-going MCP will produce, on average, 5 × 10 6 × ε 2 relevant Cherenkov photons. Next we need to estimate how many of these Cherenkov photons could reach an IceCube module, this implies taking into account the absorption of the Antartic ice [79] and the single-photoelectron efficiency [82]. The number of photons from a single-point light source, which would be a given MCP energy loss, is reduced by a factor of 10 −3 at a 10 m distance and by 10 −6 at a 100 m [83]. Thus, if an MCP passes at a 10 m distance from an IceCube DOM the Cherenkov photon yield will be approximately 10ε 2 photons. This estimation implies that the relevant MCP parameter space that can produce enough Cherenkov photons in IceCube is above approximately 5 GeV, cf. Fig. 7.
In the above discussion, we have only taken into account the Cherenkov light, but there are other sources of light production from charged particles, such as ionization. To study the sensitivity of neutrino telescopes to millicharged particles, while taking into account these other losses, we have developed a dedicated Monte Carlo to estimate the trigger rate in detectors such as the IceCube Neutrino Observatory in the South Pole. This Monte Carlo can be found in our GitHub repository at this https URL.
In order to estimate the trigger rate we need to know the precise locations of the energy depositions of the MCPs and their distance to the detection units. Our Monte Carlo consists of the following stages: 1. We produce MCPs according to a power-law energy distribution and spread them uniformly across the surface of the Earth.
2. We propagate the MCPs from the Earth surface to a cylinder that contains the detector instrumented volume.
To perform this we use PROPOSAL, which is a Monte Carlo package that simulates the energy losses of charged leptons in different media. The energy losses implemented in PROPOSAL are equivalent to the ones discussed in Sec. III, however, unlike our previous discussion which focused on the mean energy loss, PROPOSAL provides a detailed simulation of continuous energy losses and stochastic ones. Thus, in this step, for each MCP particle in our Monte Carlo we obtain the location, type, and amount of energy loss along the particle trajectory.
3. We convert the energy losses produced in the detector vicinity to the number of Cherenkov photons produced in ice. In order to do this, we use the publicly available PPC photon propagation code, which uses the parameterizations given in [64] to convert each type of loss (ionization, bremsstrahlung, photo-hadronic, and pair-production) into an ensemble of Cherenkov photons.
4. The in-ice photons are then propagated by PPC through the detector instrumented volume. Within this volume, we specify the detection units in PPC by providing the coordinate of each sensor. With this information PPC returns the number of detection units hit by photons.
Using this Monte Carlo approach one can estimate the trigger rate in a given detector. However, this is beyond the scope of this article and instead we provide this Monte Carlo and MCP fluxes computed in Sec. II in order for experiments to estimate the MCP yield in their specific setup.

VI. DISCUSSION AND CONCLUSIONS
Whether particle charges are quantized and whether any charged particles exist beyond the Standard Model are two questions that have been asked for generations. Searches for new particles with fractional charges work to address both of these questions, and significant progress has been made in these searches in recent years, particularly in the MeV to GeV mass range.
Focus in this mass range has been divided between two general categories -searches for millicharged particles produced by collider and fixed-target experiments, and searches for millicharged particles naturally produced in the atmosphere. We have focused on the latter approach. In this work, we have revisited current constraints from the Super-Kamiokande neutrino experiment, qualitatively confirming the existing literature on searches of this type. We have also analyzed existing data from the XENON1T dark matter direct-detection experiment and projected future capabilities of the JUNO reactor neutrino oscillation experiment in this parameter space, demonstrating paths for improvement in the next decade.
Going beyond this, we have combined a number of millicharged particle search strategies by proposing the search for multiple-scattering events in a liquid scintillator detector (specifically JUNO), allowing for sensitivity to significantly smaller millicharges than conventional, single-scattering searches. This is because the multiple-scattering searches allow for analyses to probe even smaller energies where backgrounds dominate the conventional search.
We have focused predominantly on searches for these particles in tens-of-kiloton-scale detectors. However, even larger detectors, such as the IceCube Neutrino Observatory (and forthcoming IceCube Upgrade) can offer an interesting, complementary means of searching for millicharged particles. We demonstrated that detectors with longer path lengths (of traversing millicharged particles) are well-suited for multiple-scattering searches. IceCube is one such detector, and it can potentially perform a search for these "faint track" signatures in the coming years. We have developed a Monte Carlo package to simulate the propagation and energy deposition of millicharged particles through a detector such as IceCube.
As we look forward to the next decade of searches for new, fractionally charged particles, it is imperative that we combine as many search strategies as possible. This maximizes the chances of discovery, and, in the hopeful event of one discovery, a combined approach is our best way to interpret and understand such a momentous result. Atmospheric searches, particularly those looking for multiple-scattering events, offer a powerful means to search for these particles, complementary to current and upcoming collider and fixed-target based searches.
where Br(m → e + e − ) is the experimentally-measured branching ratio of the meson into electron/positron pairs and I (2) (x, y) is a dimensionless quantity relating these two processes, The energy distribution of the millicharged particles in the parent meson rest frame is flat. The lab-frame distribution can be obtained by transforming between frames using the boost of the parent meson. The allowed energy of the millicharged particle in the lab frame E χ can be determined by requiring where λ(x, y, z) is the Källén function [84] λ(x, y, z) = x 2 + y 2 + z 2 + 2xy + 2xz + 2yz. (A4) Three-Body Decays: For the pseduoscalar meson (π 0 and η) decays, the process of interest is m → γχχ. The branching ratio for this decay can be related to the branching ratio of m → γγ using where I (3) (x) is, similar to I (2) (x, y), a dimensionless function [24], To obtain the lab-frame distribution of the millicharged particle energy, we use the quantity z ≡ E χ /γ m , where γ m is the meson boost. This distribution can be expressed as For a three-body decay we can determine the allowed range of E χ using with E max ≡ When discussing expected signal and background event distributions in Super-Kamiokande (cf. Fig. 6 left), we showed the expected distributions with Super-Kamiokande II for simplicity. For completeness in Fig. 10, we provide the analogous distributions for all three stages of Super-Kamiokande data collection that enter our analysis. Each panel displays the expected background events (blue), signal assuming m χ = 300 MeV and ε 2 = 5 × 10 −5 (orange), and their sum (green), compared against the observed data (black crosses).
In practice, we do not use the expected background distributions (blue) in our analysis. This is because they do not include the possibility of any diffuse supernova neutrino background events, which would contribute at low recoil energy. For this reason, we adopt the background-agnostic, one-sided likelihood approach discussed in Section IV A 3.  10. Distributions for our millicharged particle search for the three stages of Super-Kamiokande. Data is shown as crosses and expected background is shown as a blue histogram. An example millicharge particle distribution is shown, for the same parameters, as an orange histogram.

Appendix C: XENON1T Preferred Parameter Space
In Section IV A we discussed constraints on millicharged particle parameter space as a function of m χ and ε 2 coming from single-hit searches, including Super-Kamiokande, XENON1T, and a future search in JUNO. When considering current data from Super-Kamiokande and XENON1T, we took a background-agnostic approach with a one-sided likelihood function to derive a conservative upper limit on ε 2 as a function of m χ . For XENON1T, this was done in part due to the much-discussed potential tritium background, absent in the nominal B 0 background model of Ref. [73].
In this appendix, we perform an alternate test of the XENON1T data where we assume that the B 0 model is robust and calculate the likelihood comparing µ s i (the expected number of signal events in bin i) plus B 0,i (the expected background in bin i) with the observed data d i . This process, due to the excess electron-like events in the lowest energy bins, will yield a preferred region of millicharged parameter space instead of a constraint.
The potential signals proposed in Ref. [73] include neutrino magnetic moments and solar axions, which improve the fit to data over the background-only explanation by 3.2 − 3.4σ. We find that the millicharged particle signature improves our test statistic over the background-only hypothesis by 10.43 units. When we (conservatively) consider two degrees of freedom, this implies a preference of 2.8σ over the null hypothesis -if we had considered the highlycorrelated (m χ , ε 2 ) to account for a single degree of freedom, this would imply a preference of 3.2σ, as preferable as the hypotheses presented by the XENON1T Collaboration. Fig. 11 presents our 1σ (yellow) and 2σ (green) preferred region of parameter space from XENON1T, compared against all other existing limits (including our Super-Kamiokande analysis) in grey. We see that the bulk of the ±2σ preferred parameter space is already excluded by one or more other constraint, however, some regions survive near m χ ≈ 350 MeV. These regions will be tested by JUNO's single-hit analysis and thoroughly explored by JUNO's double-hit analysis.
We note here, as in the main text, that our XENON1T analysis and signal-event prediction rate from Eq. (9) do not include the fact that the electrons in the XENON1T detector are bound [31,74]. This is relevant because their binding energies are non-negligible compared to the recoil energy observed in the detector, and this implies that our signal rate predictions are optimistic. For this reason, in addition to the tritium background discussion above, we caution the reader from inferring that the MCP solution to the XENON1T excess is a likely explanation. . Preferred parameter space from XENON1T. Yellow/green bands depict the ±1σ/±2σ preferred parameter space of our XENON1T analysis when we assume that the B0 background of Ref. [73] accurately represents all backgrounds. All existing constraints are shown in grey.