Millicharged Cosmic Rays and Low Recoil Detectors

We consider the production of a"fast flux"of hypothetical millicharged particles (mCPs) in the interstellar medium (ISM). We consider two possible sources induced by cosmic rays: (a) $pp\rightarrow$(meson)$\rightarrow$(mCP) which adds to atmospheric production of mCPs, and (b) cosmic-ray up-scattering on a millicharged component of dark matter. We notice that the galactic magnetic fields retain mCPs for a long time, leading to an enhancement of the fast flux by many orders of magnitude. In both scenarios, we calculate the expected signal for direct dark matter detection aimed at electron recoil. We observe that in Scenario (a) neutrino detectors (ArgoNeuT and Super-Kamiokande) still provide superior sensitivity compared to dark matter detectors (XENON1T). However, in scenarios with a boosted dark matter component, the dark matter detectors perform better, given the enhancement of the upscattered flux at low velocities. Given the uncertainties, both in the flux generation model and in the actual atomic physics leading to electron recoil, it is still possible that the XENON1T-reported excess may come from a fast mCP flux, which will be decisively tested with future experiments.


I. INTRODUCTION
Cosmic rays exhibit a rich phenomenology in the Earth's local neighborhood. Highly boosted charged particles strike protons in the interstellar medium and in the upper atmosphere, to produce showers of more numerous (and less energetic) charged particles. Any charged particle which is sufficiently light and longlived will be present in cosmic ray showers, including the positron, the muon and the pion which were first discovered there [1][2][3]. In addition, charged cosmic rays interact with local magnetic fields and undergo a random walk, described well by a diffusion process [4], which leads to an enhancement in their abundance near the galactic disk. As a result, cosmic rays, particularly slow moving particles, are retained in the vicinity of the galactic disk. Like any charged particles, incoming CRs will interact with matter through Rutherford scattering, dσ/dT e ∝ 1/T 2 e , with soft collisions being most abundant. Upon reaching the dense environment of the atmosphere, the Earth itself, or a detector, CRs will loose energy rapidly. This allows for their detection with tracking and calorimetry, but also allows one to protect detectors from CR background with shielding and underground laboratories.
Millicharged particles (mCPs) are a simple extension of the Standard Model which introduce just that, a new particle with a very small charge. They can arise in simple frameworks in which a new gauge boson kinetically mixes with the photon [5]. mCPs in the MeV to several GeV range have been identified as interesting targets for experimental analyses, with limits set by SLAC's milliQ [6], milliQan [7,8], LSND [9,10], Fermilab's miniBoone [9,11] and ArgoNeuT [12,13], and proposed dedicated experiments [14][15][16][17]. Importantly, for this paper, mCPs exhibit all of the same phenomena as cosmic rays, only to a lesser degree due to their small charge.
In addition to serving as a simple extension of the SM, mCPs can also have important astrophysical and cosmological consequences. mCPs have been invoked to explain several experimental excesses/anomalies [18][19][20][21][22][23]. Earlier work has considered connections of mCPs to dark matter and its observable consequences [24][25][26][27][28][29]. More recently, there has been a resurgence of interest in mCPs, and specifically millicharged dark matter (mDM) [30][31][32][33][34][35][36][37][38][39], since the EDGES collaboration reported an anomaly in their 21 cm data [40]. The proposal is that the enhancement of mCP-electron scattering at low velocities allow for a late thermal re-coupling of baryons and dark matter, lowering the kinetic temperature of atoms, and thereby enhancing the 21 cm absorption feature. This would be in line with stronger-than expected absorption claimed to be seen by the EDGES collaboration. Because the universe is awash with electromagnetic fields, the ultimate fate of mDM is a non trivial question and the current literature is in a state of flux with no clear consensus on either the velocity or density distributions of mDM. Some authors have argued that mDM is FIG. 1. A cartoon depiction of the leaky-box/disk-diffusion model [47]. A cosmic ray collides with a proton in the ISM producing an mCP pair via an intermediate meson. This results in a fast-moving mCP that diffuses slowly about the leaky-box. The long-retention time due to magnetic diffusion results in a greatly enhanced flux of mCPs sourced from the ISM [4,48]. evacuated from the disk entirely by supernova shocks [24] so as to be unobservable in terrestrial experiments [41,42], while others have argued that the same supernova shocks, coupled with diffusive transport eventually lead to an easily observable high-velocity distribution of mDM although estimates vary by orders of magnitude c.f. [43] v.s. [44]. More recent work has suggested that there are a variety of unaccounted for collective damping effects that could severely limit the speed of mDM reaching a presentday terrestrial detector [45]. What is unambiguous, however, is that if mDM existed at the cosmic dawn and is responsible for the EDGES anomaly, then there should be a predictable (albeit model-dependent) density of mDM spread somewhere within roughly 20 kpc of our local galaxy.
Motivated by the uncertainty in the predictions for a mDM density and velocity distribution, it is natural to ask if there is a direct-detection scheme that is relatively insensitive to the details of mCP transport. We answer this question in the affirmative by considering two qualitatively different scenarios: (a) The mCPs do not constitute a component of the dark matter, but are readily produced via collisions of cosmic rays with the interstellar medium (ISM) in addition to local cosmic-rays' collisions with the upper atmosphere (as discussed in [46]).
(b) The mCPs do form a sub-component of the dark matter in the ISM and are boosted via cosmic ray collisions. We refer to this scenario as "Rutherfordboosted" dark matter (RBDM).
Scenario (a) is entirely independent of mDM modelling assumptions and can be used to set limits on mCPs in the same way a fixed-target or collider experiment could. Scenario (b) assumes an ambient population of mDM, but is relatively insensitive to the details of the number density and velocity distributions. Even if mDM is slow, cosmic-ray collisions will generate a RBDM component that is fast. Even if the mDM is evacuated from the disk, provided it lies within the cosmic-ray diffusion zone, cosmic rays will find it. Finally, even if collective damping effects are important for coherent mDM motion, the RBDM component will be incoherent having been generated by random and isotropic collisions throughout the volume of the cosmic ray diffusion zone. In both scenarios our results rely crucially on the diffusive motion of mCPs (see Fig. 1) which leads to a residency time in our local galaxy that is four orders of magnitude larger than would be expected if the mCPs were to exhibit ballistic motion. This can result in the ISM component of the flux being larger by an order of magnitude then the locally generated flux from the upper atmosphere despite the ISM being a "thin target" and the atmosphere being a "thick target".
Rates of mCP detection naturally benefit from lowthresholds, which both increase the detection cross section and allow an experiment to access low-velocity components of an mCP flux.
The XENON1T collaboration has demonstrated the ability to conduct low-background tonne-scale electron recoil experiments with a threshold in the 100-eV (S 2 ) [49] to few keV scale (S 1 -S 2 ) [50]. As we show below, such low thresholds are especially advantageous for Scenario (b). A natural question, then, is whether or not XENON1T's unique low-threshold capabilities are able to overcome the large discrepancies in fiducial volume between itself and Super-Kamiokande (1 t vs 22.5 kt), which has recently been identified as a resource for studying mCPs from cosmic rays [46].
Our main results demonstrate that in Scenario (a), single electron recoil in the Super-Kamiokande still provides superior sensitivity even compared to the most sensitive dark matter detectors, such as XENON1T. We update the sensitivity/exclusion curves by accounting for the additional channel of mCP production in collisions of CRs with the interstellar medium. In Scenario (b) we exploit the fact that the flux of RBDM has a significant infrared enhancement, and are able to derive competitive limits from the XENON1T experiments' electron recoil analysis. Given enough uncertainty in the current treatment of atomic physics, we find that the claimed few keV-scale excess of recoil electrons may still plausibly come from the fast flux of millicharged particles, but this subject requires a more detailed investigation. Departing from these two relatively conservative scenarios, we explore the sensitivity to a fast mCP flux of more exotic origin, such as e.g. dark sector decay to mCP, finding that in this case the sensitivity may extend to very small values for the millicharge.
The rest of the paper is organized as follows: in Section II we introduce the well-known, but crude, "leaky-box" model of magnetic retention that allows the ISM to contribute fluxes comparable to those coming from the upper atmosphere. Then, in Section III we outline the calculation of cosmic-ray fluxes from 1) the upper atmosphere and 2) the ISM. We also discuss the resulting RBDM flux in the case where mCPs form a fractional population of the DM. Next, in Section IV we describe how event rates at the XENON1T detector are calculated and present exclusions and projections for both Scenario (a) and Scenario (b), as well as a generalization to the fast mCP flux of unspecified origin. In Section V, attenuation due to the Earth overburden is discussed. This introduces novel dependent scaling in the case of Scenario (b). Finally in Section VI we summarize our results and comment on the future impact of low recoil detectors for mCP and mDM searches.

II. MAGNETIC RETENTION
It is well known that magnetic fields lead to a very large enhancement in the cosmic ray intensity [4,48]. The essential physics is that magnetic fields scatter particles in such a way that they undergo a random walk described by a diffusion coefficient, D, that depends on the magnetic rigidity, R = p/|q| of the particle at hand (more on this later). This statement is robust, and so the magnetic retention of mCPs can be readily obtained by re-scaling the well understood magnetic retention of protons. In what follows we do not attempt an elaborate treatment (as in e.g. GALPROP) relying instead on a simplified diffusion model.
We will quote all of our results in terms of "leaky-box" or "disk-diffusion" model wherein a particle is envisioned to diffuse about a leaky-box roughly conceptualized as a cylinder of radius ∼ 20 kpc and a height of ∼ 5 − 10 kpc; a cartoon of is shown in Fig. 1. This defines a length scale LB which conveniently parameterizes the enhancement of the mCP flux due to magnetic retention. This parameterization can be extended to account for a more sophisticated treatment of the mCP diffusion as we discuss in Appendix A.
The residency time within this container is given by τ LB ∼ 2 2 LB /D where D = D 0 β(R/R 0 ) δ is the rigidity dependent diffusion constant and the factor of two comes from the conventional definition in a diffusion equation Here, R 0 is a reference rigidity and δ is an exponent that is somewhat model dependent; a popular choice that is both theoretically well motivated and empirically successful is based on Komolgorov scaling, δ = 1/3 [47]. This can be compared to the residency time of free-passage τ FP = LB /β. Realistic values of LB would range from 1−10 kpc, while D 0 ≈ 3 × 10 28 cm 2 s −1 and R 0 = |e| × [1 GeV].
The ratio of the leaky-box retention time τ LB , to time of free-passage (i.e. ballistic motion), τ F P determine the rigidity-dependent magnetic-retention enhancement of the mCP flux This demonstrates the dramatic (four orders of magnitude) enhancement of mCPs that are sourced in the ISM relative to the naive estimate neglecting diffusive motion. We use LB to quantify the relative size of magnetic retention, and we (again) encourage the reader to interpret this a convenient parameterization of the enhanced mCP flux relative to the naive "free-passage" calculation rather than as a bonafide physical parameter. A more precise definition of R B accounting for the spatial distribution of targets, and the cosmic ray density is given in Appendix A. Foreshadowing slightly, this enhancement easily compensates for the relative paucity of scattering targets allowing mCP production in the ISM to be competitive with local production in the upper atmosphere. For instance, taking a cross section of 40 mb=4 ×10 −25 cm 2 and an ISM column density of n ⊥ ISM =(1 cm −3 )×(1 kpc)= 3 × 10 21 cm −2 , we find a probability of inelastic scattering to be ∼ 10 −3 easily enhanced to an O(1) number after accounting for magnetic retention.
Lastly, it is important to note that there are broadly two versions of mCPs in literature and magnetic retention works only for a subset of these models. In the first model, mCPs are charged directly under the SM U(1), hence the mCPs behave very similar to the SM charged particles at all length scales and the results of this section hold. Particles charged under a dark U(1) which kinetically mixes with the SM U(1) also carry effective charge when probed at energies ω m A , the mass of the dark U(1) photon. In this model, if the dark photon mass is much larger than the inverse galactic length scales, there is no magnetic retention. For this reason, the main results of our paper are applicable to the first model described above, but wherever relevant we also show results in the absence of magnetic retention to best describe limits on mCPs with a heavy dark photon mediator.  [51] (triangles) as well as the parameterization of [52] (squares).

III. FLUXES OF MCPS
Having discussed the retention of mCPs due to astrophysical magnetic fields, we now turn to their nascent production. As outlined in the introduction we consider two conceptually different scenarios and organize the following discussion accordingly.

A. Scenario (a): cosmic-ray produced millicharges
As has been recently demonstrated in [46], cosmic rays can efficiently produce mCPs in the upper atmosphere. We consider the sensitivity of XENON1T to this production channel, but also consider a secondary flux component stemming from the ISM. Much of the methodology used in this paper is outlined in detail in [46] and so we discuss only the most important details here and refer the interested reader to our earlier work.
There are two main ingredients in the production of mCPs from cosmic ray collisions: 1) the production and resulting energy distribution of electromagnetically decaying mesons, and 2) the resultant spectral shape of the millicharged flux from the subsequent meson decays in flight. The former is much more challenging than the latter, but can be well approximated if one has accurate Feynman-x distributions for the production of mesons in pp collisions [46]. Once these are known it is a straightforward numerical exercise to boost the spectra into the lab frame weighted by the meson energy spectra obtained in step 1.
In the upper atmosphere a conservative estimate of the meson yield can be obtained by considering only primary production, the probability of which is given by the ratio of the meson-production to total-inelastic cross section σ m /σ inel . In the ISM things are slightly different and this factor should be replaced by σ m n ⊥ ISM where n ⊥ ISM is the column density of the interstellar medium. This difference stems from the fact that the atmosphere is a thick target (λ Atm ∼ 1 km with λ = [σn] −1 the mean free path) whereas the ISM is a thin target [λ ISM ≈ 8 Mpc]. For the column density of the interstellar medium we take This corresponds to the average density of protons in the ISM multiplied by 1 kpc, which we take as a representative distance through which a cosmic ray will transverse the ISM.
In [46] spectra were quoted as differential with respect to γ. A simple change of variables (including a Jacobian) converts these spectra to being differential with respect to β χ γ χ . The mCP intensity can then be found (assuming isotropy) by multiplying the rate of production, by the column density of targets and the enhancement factor coming from magnetic retention, and dividing by 4π We note that while we have derived this expression within the context of an overly-simplified leaky-box model, a slightly more systematic treatment can be found in Appendix A, where n ⊥ ISM × R B can be identified as the second bracketed term in Eq. (A11).
R prod is the production rate per ISM target particle, and is given by where the factor of two accounts for the combined χ andχ flux. The conditional probabilities are determined by dσ m /dx F and the details of m → χχ decays respectively [12,46]. The flux of mCPs coming from the upper atmosphere is defined similarly but without the factor of n ⊥ ISM ×R B in Eq. (3) and with σ m being divided by σ in , the total inelastic pp cross section. Importantly the flux from the upper atmosphere scales as 2 since branching ratios involve two powers of , whereas the flux from the ISM flux scales as 2+1/3 (for Komologorov scaling) and has a slightly different spectral shape due to the dependence of R B on magnetic rigidity.
In Figure 3 we show the flux of mCPs in this scenario for = 3 × 10 −3 , and m χ = 200 MeV. The fluxes for Benchmark values of mχ = 200 MeV, and = 3 × 10 −3 have been chosen corresponding to viable parameter space to explain the EDGES anomaly as depicted in Fig. 1 of [30]. production in the upper atmosphere and in the ISM (including magnetic retention) are shown. Both fluxes peak at a few GeV, but as advertised, the ISM component is about an order of magnitude larger and slightly softer.
Finally, we would like to emphasize that in the upper atmosphere we approximated pA collisions as a sum of incoherent pp collisions, and we neglected secondary and tertiary meson production, thus our approach is conservative. In the ISM, however, collisions genuinely are pp collisions, and since the column density is much smaller than the inverse cross section, secondary production is negligible. Therefore our treatment of the ISM production is robust, with the only (but important and sizeable) uncertainties stemming from our treatment of magnetic retention.

B. Scenario (b): Rutherford boosted dark matter
Having established the ISM as a competitive source of mCPs from pp collisions, it is interesting to posit that some fraction of DM is millicharged (mDM) such that, with no additional assumptions, cosmic rays would boost this component of the DM via Rutherford scattering. Rutherford scattering is enhanced by a 1/β 2 p factor such that scattering via low-energy cosmic rays dominate the RBDM spectrum. It is therefore imperative that we use reliable data on the cosmic ray spectrum in the ISM at low (i.e. 10 MeV scale) energies. The CR flux in the ISM is known to differ substantially from measured "top of atmosphere" fluxes at low energies. Fortunately Voyager 1 has both measured the cosmic ray spectrum beyond the heliosphere [51] and provided the data in the energy domain we require to estimate the flux of RBDM. We use a simple analytic interpolation of their data for the cosmic ray spectrum in the ISM given by This matches on to more sophisticated parameterizations at higher energies [52], and reproduces the data reasonably well at lower energies as can be seen in Fig. 2. Since we anticipate that the size of the magnetic retention introduces a much greater uncertainty in the flux, Eq. (5) is adequate for our purposes. The rate of RBDM production per mDM target is given by Here, T min is a function of T χ that is found by demanding the incident proton have enough energy to kick an mCP to a kinetic energy of T χ . By solving for the maximum four-momentum transfer possible one can find The flux of RBDM is found, in direct analogy with Eq. (3), by multiplying by the number density of mDM, and the magnetic enhancement factor where we have parameterized the mDM column density as We have taken a longer column length for DM than for the ISM contribution (1 kpc vs 10 kpc) because while the local DM density is roughly equal to the local ISM density, DM is spread over a larger volume. In this paper we define the mDM fraction f χ globally with respect to the Milky Way As discussed in the introduction, mDM transport can lead to over-or under-densities of mDM in our local vicinity. We can imagine two extreme cases: One could imagine that mDM falls back onto the disk [44], and has a large enough transport cross section with baryons that the mDM matter profile closely tracks the baryonic profile in the disk. Alternatively, if mDM is ejected it could be present in the cosmic-ray diffusion volume, but be well separated from the disk (for instance ∼ 5 kpc out of the plane). We estimate that even in these extreme scenarios the re-distribution of the mDM density only leads to an O(10) enhancement of the RBDM flux arriving at Earth (see Appendix A). We can transform from T χ to β χ γ χ by including the Jacobian dT χ /dβ χ γ χ and multiplying through by the magnetic retention factor R B where Putting everything together, and allowing mDM to be some fraction, f χ , of the total DM density, we estimate the intensity of RBDM arriving at earth to be where g(β χ γ χ , m χ ) is a dimensionless function that depends on m χ and satisfies g(0, m χ ) = 1; Fig. 4 shows its behavior as a function of β χ γ χ for a few representative masses.
As an illustration of the various components of the flux arriving at a terrestrial detector (without considering effects of attenuation due to e.g. overburden or solar winds) we show the differential flux in Fig. 3 for the benchmark point where one part per million of dark matter is millicharged f χ = 10 −6 , = 3 × 10 −3 , and m χ = 200 MeV (this corresponds to an EDGES preferred point of parameter space as presented in Fig. 1 of [30]). We compare this to a Maxwellian velocity distribution cut-off at 575 km/s. RBDM adds a hard tail of fast millicahrged dark matter, as expacted.

C. Exotic Scenarios for large fluxes of millicharges
Even larger fluxes of mCPs can be motivated in more exotic scenarios than the fairly minimal Scenarios (a) and (b) considered thus far and we pause here to discuss some possibilities. The question of maximal flux of exotic particles (mCP or other variety) is quite general and deserves a separate discussion. Some exotic sources of fast moving DM include boosted DM [53,54], solar reflection [55] and SM decays [56]. For concreteness we consider here, scenarios with minimal dark sector content. The sources of exotic flux with v > 0.1 we consider are: i. dark matter decay globally and in the galaxy; ii. dark matter decay or annihilation inside astronomical objects such as the Sun and solar system planets; iii.
exotic radioactivity from the decay of normal atoms.
Scenario i. posits that throughout a history of the Universe there existed, or still exists, an unstable agent that produces pairs of mCPs, according to X → 2mCP (see [57] for more extended discussion of similar scenarios). If X constitutes the dark matter today, 5% or so of decays after the recombination until now is permissible. Notice that, for example, X can be a subdominant component of the dark matter that has fully decayed to mCPs. Then the mCP's spectrum, which is presumably (semi)-relativistic initially, can be significantly softened by the subsequent expansion. Therefore one can find a situation where only the DM experiments would be sensitive to this type of flux, while for neutrino experiments the flux would be below the energy threshold. Scenario ii. occurs when DM is effectively intercepted by large astronomical bodies, such as thee Sun, with subsequent annihilation or decay.
In this case, the maximum flux will be limited by the "supply" of DM, that is the maximum number of particles intercepted by the Sun with the subsequent annihilation to mCPs. Scenario iii. draws on recent ideas in [56,58] where Hydrogen atoms undergo decays to exotic particles, which in turn can lead to the mCPs. This scenario (H-atom decay) is somewhat less constrained that the one based on proton decays, as constraints on H-decay are less stringent than in the proton case (see [56,58] and references therein).
In each of these three cases we can estimate the maximum amount of mCP flux that can be created. Assuming a GeV-mass progenitors for the mCPs, we get: In comparison, the total virial DM flux is Φ vir Thus we see that the nonminimal options for obtaining "fast mCP flux" could be as high as ∼ 0.1% of the main DM flux. In addition to these exotic, but relatively straightforward, scenarios one can imagine slightly more elaborate mechanisms for generating a fast flux of mCs, although each of these will require further scrutiny. For example, a potentially interesting effect is related to the formation of nucleus-mCP bound states for the negatively charged particles. This possibility refers to larger mass and larger part of the parameter space.
Taking, for example, helium nucleus into consideration, one can estimate the binding energy to be ∼ (µ/(2 GeV)×( /0.05) 2 ×133 eV, where µ is the reduced mass of the mCP-Helium nucleus pair. Binding energies of this size will lead to post-BBN, but pre-recombination formation of bound states. These bound states, if surviving until today, can be accelerated at astrophysical sites almost as efficiently as other nuclei. Moreover, upon entering Earth's atmosphere, the mCPs in such bound states can be easily stripped by nuclear collisions, creating a flux of energetic mCPs. Interestingly, even if the binding of mCPs to nuclei is mediated by a massive dark photon, such that magnetic retention of RBDM is inhibited (because electromagnetism is shielded at large distances), this "nuclear delivery system" mechanism would still be viable as it relies exclusively on the non-zero electric charge of the He-mCP atom. While this mechanism is certainly interesting, it could potentially be prone to many uncertainties, and chiefly among them the abundance of bound states in specific astrophysical settings, several billion years after the Big Bang. Therefore we do not pursue it further in this paper.
At smaller values of , mCP-nucleus bound states will not form. However, one could imagine mCP-mCP bound states to be a dominant DM subcomponent provided that there are additional dark sector dynamics that bind mCPs to one another. For example a dark photon could bind two oppositely charged mCPs (possibly of two different species with different masses as considered in [59]). Binding energies could be sufficiently large that bound states would form in advance of recombination, and in this case bounds from the CMB [35] would need to be reinterpreted. A substantial fraction of DM could then be composed of neutral mCP-mCP bound states, which would later be ionized by CR scattering events leading to a larger effective value of f χ 0.004. This would contribute to (and thereby enhance) the RBDM flux. The viability of this scenario is highly model dependent, and since it would require further dark sector dynamics to facilitate mCP-mCP bound state formation we do not consider it further here.
In addition to these exotic scenarios, a larger-thanestimated enhancement from magnetic retention is easily envisioned. For example if the diffusion constant were to rise more steeply at low velocities than we have modelled (e.g. D ∼ Dβ η (R/1 GeV) δ with η ≤ 1 as suggested in e.g. [60][61][62][63]) then the RBDM flux would be even further enhanced in the infrared. A similar effect could be achieved by considering Kraichnan (δ = 1/2) rather than Komologorov scaling. Also, as discussed in Appendix A, the distribution of mDM itself can lead to a large magnetic retention. If the mDM is constrained to lie in the galactic disk rather than in an NFW-like profile then the overlap of mDM targets with cosmic rays is greater, and the local flux is enhanced because all of the RBDM sources are closer to Earth.
It is also possible that low-energy cosmic rays have a much larger intensity in substantial regions of the Milky Way than the measurements of the Voyager-I mission would suggest. Indeed, authors have speculated in the past that the cosmic ray spectrum could be several orders of magnitude larger at low energies, and that such a spectrum could explain unexpectedly large concentrations of H + 3 [64]. Since it is the low-energy part of the cosmic-ray proton spectrum that dominates the RBDM flux, if the low energy (sub GeV) cosmic ray intensity were to be this large anywhere in the galaxy with a significant mDM population, then it lead to substantial increase in the RBDM flux arriving at Earth relative to the estimates presented in this paper.

IV. RESULTS
Having established momentum dependent fluxes for both Scenario (a) and Scenario (b), we now turn to setting limits from the resultant rate in the electron ionization analyses in XENON1T as well as projections for the future. The differential cross section for ionization involves atomic physics form-factors which are described in Appendix B. We estimate the rate of electron recoils  [49] are shown in black. Also shown is the parameter space that explains the excess events in the 1-7 keV bins in S1-S2 data [50] in red. In gray we show a compilation of current terrestrial limits on MCPs [46] and also include the recent milliQan-demonstrator results [8] .
in an energy bin as where eff(E t ) is the electron efficiency as a function of energy transfer E t and the differential cross section dσ dEt is provided in Appendix B. The factor of 2π comes from the fact that we only consider zenith angles above the horizon. The minimum velocity is set by a combination of the experiment's electron recoil energy threshold and overburden attenuation. This latter effect is dependent on 2 /m χ and introduces some non-trivial dependence of the rate on for Scenario (b) (RBDM) as shown in Fig. 10. We therefore defer a discussion of the overburden until after presenting our results. We consider two analyses from XENON1T in this work. First, we consider the S 2 only analysis from [49] with 22 tonne day exposure. Since the S 2 only analysis has poorly understood systematic uncertainties we only use it for limit setting purposes. Following [65], we take 22 events in the (0.2 − 0.5) keV bin and 5 events in the (0.5 − 1) keV bin as limiting rates (all parameter space leading to larger event counts is excluded). The efficiency eff is assumed to be unity. The second analysis considered is the recently reported excess in the S 1 -S 2 analysis from [50] with 0.65 tonne year exposure. The efficiency eff, is taken from Fig.2 of [50]. The collaboration observed 53 excess events in the 1−7 keV bins with the excess peaked at the (2 − 3) and (3 − 4) keV bins.
A. Scenario (a): Cosmic-ray produced mCPs As seen in Fig. 3, the flux of mCPs in Scenario (a) is peaked at (semi-)relativistic βγ. This allows almost the entire spectrum to penetrate the overburden of 3650 mwe, for the values of we consider here. As we will see in Sceanrio (b), at lower velocities the overburden plays an extremely important role in determining the experimental sensitivity.
First, we compare limits from the different XENON1T analyses in Fig. 5. We show limits obtained from the two S 2 bins in black. The red curve corresponds to the 53 event excess in 1 − 7 keV in the S 1 -S 2 analysis. We find that for parameters capable of reproducing the 53 event excess in the 1−7 keV bins, the expected signal in the 0.5-1 keV bin of the S 2 only analysis is roughly 6−7 events (as against the 5 event limit). Important to our predictions is a treatment of the atomic physics governing the detector response that differs from commonly advocated approximations in the literature. We discuss the details of our approach in Appendix B and point out that a plane-wave treatment with a Sommerfeld/Fermi function overestimates the event rate, and dramatically so at low recoil energies.
In addition to (slightly) overpredicting the S 2 only event yield in the 0.5−1 keV bin, the preferred parameter space of the S 1 -S 2 analysis is firmly ruled out by existing constraints on mCPs. This includes the non observation of cosmic-ray produced mCPs from the upper atmosphere in Super-Kamiokande [46], beam-induced mCPs in ArgoNeuT [13], and recent results from the milliQan demonstrator [8] (shown together in shaded gray in Fig. 5). Figure 6 (Left) compares the exclusions with and without magnetic retention for the (0.5 − 1) keV bin (S 2 only). The modest gain in flux due to retention (relative to the component produced in the upper atmospheere) exhibited in Fig. 3 translates to stronger exclusions. Limits on mCP models with a light, but not massless, dark photon are equivalent to the line without magnetic retention, because, although laboratory scattering events are unaffected, on large scales the dark photon mass screens astrophysical magnetic fields nullifying the effects of magnetic retention.
Finally, in Fig. 6 (Right), we show projections for models that predict a higher flux, or for future analyses with higher experimental sensitivity. We find that if the experimental sensitivity can be improved by a factor of ten then new constraints can be set on mCPs independent of any assumptions about dark matter. We also note that this improved sensitivity could stem not only from improvements to the detector, but also if the flux were to be larger than estimated here because of e.g. a larger diffusion coefficient.   [49] are shown in black. Also shown is the parameter space that explains the excess events in the 1-7 keV bins in S1-S2 data [50] in red. In gray we show a compilation of current terrestrial limits on MCPs as well as DD limits and its ceiling [66]. (Right): Comparison of excess events for the parameter choice mχ = 100 MeV, = 10 −3 , and fχ = 0.004, a sample point from the left panel with the data taken from [50] .
greatly exceed that of Super-Kamiokande because of the huge gain in flux offered by the lower kinematical threshold on β χ γ χ . Second, the effect of overburden attenuation is more important here than in Scenario (a). At low velocities the incident mCPs have less energy (being more easily stopped in rock) and interact more strongly with matter (due to the 1/β 2 scaling of the Rutherford cross section).
All of our constraints depend on the fraction of DM that carries a millicharge, f χ . Constraints from the CMB power spectrum require that mDM compose only a small (sub-percent) fraction of DM [35]. A convenient benchmark parameter choice that is "safe" from the perspective of CMB constraints is f χ = 0.004 and we adopt this as a benchmark value of f χ throughout this section. This same parameter choice also ensures that the mDM we consider is not ruled out by constraints on self-interactions in galaxy cluster collisions [67].
We start by comparing the different XENON1T ionization analyses in Fig. 7 (Left) for f χ = 0.004. Much like in Scenario (a), we find that the parameter space capable of producing a 53 event excess in the S 1 -S 2 analysis (plotted in red) results in a slight overpopulation of the 0.5 − 1 keV S 2 only bin (relative to the 5 events) and a slight under-population of the 0.2 − 0.5 keV bin (relative to 22 events). Thus, the S 2 only and S 1 - Here we illustrate the dramatic consequences of magnetic retention by plotting the parameter space that yields an expected rate of 5 events in the 0.5 − 1 keV bin of the S2 only analysis. We show the parameter space accounting for magnetic retention (solid) relative to the parameter space that would be expected without magnetic retention (dashed). A RBDM flux that remains unenhanced could occur if scattering is mediated by a light, but not massless, gauge boson such that astrophysicial magnetic fields are effectively screened. We choose an unrealistically large concentration of mDM (fχ × exp = 4 with exp denoting experimental sensitivity) so that the unenhanced parameter space is visible to the reader.
S 2 rates are not in strong tension with one another, however a proper comparison will require a detailed distorted wave calculation that is beyond the scope of this work. Unlike Scenario (a), the parameter space that is being probed by the S 2 only and S 1 -S 2 analyses is hitherto unconstrained parameter space. Also shown in gray are terrestrial limits (a combination of constraints from Super-Kamiokande [46], LSND [9], ArgoNeuT [13]), recent results from the milliQan demonstrator [8], and the ceiling for virial DM in surface experiments [66]. While the tension with the S 2 only analysis is weak, we find that the parameter space capable of producing a 53 event excess [50] in the S 1 -S 2 analysis does not reproduce the lowest energy bin well as shown in Fig. 7. This inability to reproduce the lowest energy bin's event rate persists for all choices of mDM mass m χ and charge. However, a terrestrial build up of thermalized DM as considered in [68] could produce monochromatic binding energy release, when binding to Xenon thus explaining this excess.
As we have already emphasized in Section III, the flux of RBDM is enhanced by four to five orders of magnitude via magnetic retention relative to what one would expect from conventional cosmic ray boosted dark matter [69][70][71][72][73][74]. We demonstrate the consequences of this enhancements explicitly in Fig. 8 for f χ = 0.004 and a hypothetical experimental sensitivity that is 1000 times as large as the current XENON1T data release with the same constraint of 5 events in this 0.5-1 keV bin. Compared to Scenario (a) we find that the effect of retention is much more dramatic in the m χ − plane. Naively the scaling with respect to should be identical between Scenario (a) and Scenario (b) since the production cross section scales as 2 and both fluxes experience the same rigidity dependent magnetic retention factor, however this turns out not to be the case. The reason for the unexpected scaling of signal with respect to stems from the IR-enhanced RBDM flux and effects related to the overburden which we discuss in Section V. As noted above, we emphasize that the "No Magnetic Retention" line can naturally arise in models with a light, but not massless, dark photon mediator whose mass shields the galactic magnetic fields, turning off retention. Detection in XENON1T would be phenomenologically identical to an mCP, but the resultant flux would be roughly five orders of magnitude smaller.
This unexpectedly large sensitivity to the millicharge, , means that XENON1T's current sensitivity, just beyond the limits from Super-Kamiokande and ArgoNeuT, can be pushed well into regions of unexplored parameter space with relatively modest improvements.
To illustrate this point we show limits for mCPs with fixed f χ = 0.004 but with f χ ×(exp) increasing Fig. 11 with exp standing for experimental sensitivity. A large value of f χ × exp could occur either from a larger exposure (assuming fixed backgrounds), a large value of LB (implying larger magnetic than assumed here), or if some other exotic physics were to source a fast flux of mCPs. The advantageous scaling with respect to , with signal scaling as 1.86 [c.f. Eq. (20)], means that even modest improvements in experimental sensitivity, and/or modelling of magnetic retention will allow XENON1T to probe large swathes of as-of-yet unexplored parameter space in the near future.

C. XENON1T and a mono-energetic flux of mCPs
Although physically unrealistic, a monochromatic flux of mCPs is a useful example to work through because it fixes the velocity distribution and so is less sensitive to overburden attenuation. Motivated by the recent proposal of a fast sub-component of mDM as an explanation of low-energy excesses at a variety of experiments we study a benchmark velocity of β χ γ χ = 0.1. We update bounds from XENON1T (those presented in Fig. 6 of [75] assume a virial velocity and rely on nuclear recoils), and include the ceiling from overburden attenuation for XENON1T (3650 mwe).
Our results show that any parameter space capable of explaining the excess events observed at e.g. EDELWEISS via energy deposition into a collective plasmon mode, would necessarily overpopulate the S 2only and S 1 -S 2 bins from XENON1T. This statement is dependent on β χ , and for small enough values of β χ need not hold true. The electron recoil threshold ϵ

Existing Constrains
Plasmon Signal Ceiling Sensitivity of XENON1T to a monoenergetic flux of mCPs moving with a velocity of βχ = 0.1. We find that XENON1T robustly excludes an mCP-induced plasmon-excess (as proposed in [75]) because βχ = 0.1 is sufficiently fast to penetrate the overburden at Gran Sasso and kick electrons above the XENON1T thresholds. For lower velocities β ≤ 0.02, the electron recoil drops below the thresholds for XENON1T's S2-only analysis. Since a plasmon excitation requires βχ 0.00675, it is possible for a fast flux of mCPs to explain the plasmon-excess claimed in [75] without being in tension with XENON1T provided 0.00675 ≤ βχ ≤ 0.02.
at XENON1T naively sets a minimum velocity for the incident mCP of β χ 0.02. Excitation of a bulk plasmon does not require quite as fast an mCP with a velocity threshold of roughly β χ 0.007 [75]. Therefore, for an electron energy loss spectroscopy (EELS) like mechanism to be at the source of the low-energy excess in EDELWEISS, the incident mCP must have a velocity satisfying 0.007 β χ 0.02. This neglects the possibility of a coherent atomic recoil which requires a more sophisticated treatment of the atomic wavefunctions. Searches for these events could in principle further tighten the window of acceptable velocities, but would require a more sophisticated treatment that is beyond the scope of this work.

V. OVERBURDEN ATTENUATION
Since XENON1T has low recoil thresholds, it is possible for non-relativsitic mCPs to leave detectable signatures. The slowest possible velocity in the S 2 -only analysis [49] (with T e ∼ 200 eV) is β χ γ χ ∼ 0.02. At such slow velocities, however, the overburden can stop mCPs for 2 /m χ 10 −12 MeV −1 . This is especially important in the context of RBDM whose flux is dramatically enhanced as β χ γ χ → 0.
To account for this effect we make use of range tables for protons passing through water [76], and use XENON1T's overburden thickness quoted in meters water equivalent. Since dE/dx is a function of βγ alone, we re-express this quantity as We then define an "effective range" via the continuous slowing down approximation as where βγ thr is the threshold value of the experiment under consideration (for instance at XENON1T βγ thr ranges from 0.02 for the S 2 -only analysis to 0.1 for the S 1 -S 2 analysis). The true range, ∆x (in water), is then given by Inverting this expression gives a minimum initial values of β χ γ χ as a function of 2 /m χ , and we cut-off all of our spectra below this value. Our calculation agrees well with the similar βγ min analysis in [44], as well as the overburden attenuation for the specific case of virial DM (βγ ≈ 10 −3 ) in [66].
To a very good approximation the overburden serves to cut-off the flux at βγ min (we have checked by explicitly migrating the flux using the continuous slowing down approximation). Naively this truncation reduces the rate by an independent pre-factor, however this intuition fails for IR-peaked differential fluxes like the one predicted by RBDM. Instead the overburden attenuation changes the scaling of the signal with respect to .
The naive signal scaling is 2 × 1/3 × 2 with factors coming from the production cross section, magnetic retention, and detection cross sections respectively. Notice, however, that βγ min is a function of 2 /m χ , and that the event rate scales as βmin dβ I RB χ (β)/β 2 where the 1/β 2 comes from the non-relativistic limit of Coulomb scattering. Since the RBDM intensity scales as 1/β 3.33 this gives a signal that scales as (1/βγ min ) 4.33 .
In the region where 10 −11 2 /m χ 10 −7 we find that βγ min ∼ ( 2 /m χ ) 0.285 , such that the signal scales as This scaling relation of signal with respect to eventually breaks as βγ min asymptotes to βγ thr as shown in Fig. 10.
This unintuitive scaling of signal with respect to paints a highly optimistic future for the discovery (or exclusion) of RBDM using underground detectors. Rather than the rapidly falling sensitivity of 4.33 that one would naively expect, the 1.86 scaling derived above is extremely advantageous. Decreases in overburden have a non-negligible impact on signal strength. We estimate that 300 mwe of overburden would increase the signal by a factor of 50-100 relative to the 3650 mwe at Gran Sasso, however it is unclear if this is a net gain in sensitivity given the increased cosmogenic activity in the detector. Even without adjusting overburden, we expect that XENON1T can significantly extend its sensitivity. Larger iterations of large noble gas detectors (i.e. future XENON-n-ton experiments) will have the advantage of both increasing exposure, and benefiting from enhanced self-shielding. Conventional mDM searches assume a virial velocity distribution, and their sensitivity is therefore limited in many cases by an overburden ceiling at at ∼ 10 −4 . In contrast RBDM is sufficiently fast to both penetrate the overburden at Gran Sasso and to overcome XENON1T's thresholds, such that there may be room for XENON1T to probe significant regions of the as-of-yet unexplored parameter space shown in Fig. 11.

VI. CONCLUSIONS
The advent of sub-keV direct detection thresholds in large scale detectors (like XENON1T) brings qualitatively new physics to the forefront of mCP direct detection. Of fundamental importance is the realization that for m χ m e , the electron recoil threshold sets a limit on the velocity of mCPs. For recoil thresholds in the sub-keV regime, mCPs can be non-relativistic and this has three important consequences. First, the detector response is sensitive to the atomic physics of the target material and care must be taken to properly address this issue (see Appendix B). Second, the confluence of the IR-enhanced Rutherford cross section σ ∼ α 2 /β 2 and an IR-dominated flux can lead to rates that are  11. Scenario (b). Projections are shown for different fχ × exp with exp standing for experimental sensitivity.
Increased experimental sensitivity could be achieved either via detector improvements (e.g. longer exposure) or because of an enhanced flux of RBDM relative to our conservative estimate. For example a local over-concentration of mDM in the galactic disk can enhance the flux by a factor of 20 relative to our estimates (see Table I). Limits correspond to XENON1T S2 only 0.5 − 1 keV bin.
incredibly sensitive to the threshold energy (much more so than the naive 1/E thr based on the Rutherford cross section would suggest). Finally, this sensitivity to the low-velocity components of an incident mCP flux can alter the parametric dependence of the rate with respect to . The overburden determines the minimum velocity of the mCP spectrum that ultimately reaches the detector. This minimum velocity is, in turn, determined by the charge-squared to mass ratio 2 /m and therefore induces a non-trivial dependence in the event rate. In the case of RBDM studied here, this leads to a rate that scales as 1.86 rather than the naive 4.33 that one would expect after accounting for magnetic retention. This leads naturally to the second major conclusion of our study, which is that diffusive motion of mCPs can dramatically increase the flux of mCPs arriving to Earth, to the point where the ISM can outshine the upper atmosphere. The mechanism of producing RBDM outlined here is robust, and is only mildly sensitive to various mDM transport hypotheses, and the RBDM subcomponent is best conceptualized as an incoherent gas, and is therefore insensitive to collective damping mechanisms; this differentiates RBDM from mDM that has been Fermi accelerated. Cosmic ray collisions with mDM will take place provided the mDM is located anywhere within roughly 15-20 kpc of the local disk, and the intensity of RBDM arriving at Earth only varies by a factor of a few for different mDM configurations. The largest achievable flux from redistributing mDM occurs when the mDM occupies the volume of the disk, and this results in an enhancement by a factor of roughly ten as shown in Table I.
We find that XENON1T can realistically compete with large-scale neutrino detectors such as Super-Kamiokande in searches for mCPs produced via inelastic pp collisions in both the upper atmosphere and the ISM. This minimal "Scenario (a)" offers a direct probe of mCPs independent of any dark matter assumptions. Our modelling of the magnetic transport of mDM is crude, and this subject deserves a more sophisticated treatment given that XENON1T and Super-Kamiokande could both benefit from a more precise prediction of the flux from the ISM (the flux from the upper atmosphere can already be reliably calculated).
On the mDM front, our study shows that a RBDM subcomponent exists in any mDM model, and that this flux can be estimated, modeled, and its resultant sensitivity estimated. There are a number of uncertain ingredients in our formalism, however our estimates are relatively conservative: we make use of empirical measurements of the cosmic ray spectrum at low energies, we use a leaky box size of 10 kpc (rather than the 20-200 kpc suggested by Table I), and we have included reasonable treatments of the overburden attenuation and low-energy atomic physics governing the detector response. We find that XENON1T specifically is capable of probing RBDM, and that with increased exposure it can likely probe large swaths of parameter space out of reach of conventional direct detection experiments. This sensitivity reaches down "from above", i.e. for ≥(the direct detection ceiling) in contrast to conventional direct detection which is limited by an overburden ceiling; searches for RBDM may then close the gap between accelerator and fixed target searches for mCPs and virialized mDM searches in underground experiments. The scaling of signal with respect to is especially advantageous because of the IR-enhanced nature of the RBDM flux. As one probes smaller values of , larger fractions of the slow (and consequently large) RBDM flux are able to penetrate the overburden.
In conclusion, we have shown that mCP production in the ISM, magnetic retention, low detection thresholds, non-relativistic enhancements of Rutherford scattering, and non-trivial atomic physics all play a role in determining the sensitivity of the XENON1T experiment to mCPs. These conclusions apply especially to RBDM because it is the slowest available particles that dominate the rate.
In this scenario overburden attenuation introduces a non-trivial modification of the signal scaling, and a proper account of the low-energy atomic physics is essential; this issue demands more attention. Our treatment of magnetic retention can be substantially refined, and we hope to inspire future work in the cosmic ray community in this direction.
A more precise prediction of the RBDM spectrum would allow XENON1T unprecedented access to the lingering window of viable parameter space for mDM and push the experiment to the forefront in the search for mDM.
where the factor of 4π has disappeared because the intensity is per-steradian. Having calculated the rate of production per unit volume, we now turn our attention to the resultant steady-state intensity distribution.
Since I ∝ f 0 the intensity also satisfies a convective transport equation. The steady-state solution is defined via ∂ t I 0 = 0, such that the steady-state boosted DM intensity is given by Then, integrating along the line of sight defined byv we find If one takes I CR to be independent of position then this can be written equivalently as in agreement with [69], where n ⊥ DM is the line-of-sight integrated column density.
The presence of strong magnetic fields in our galaxy means that mCPs move diffusively rather than ballistically 1 and therefore dI/dt = ∂ t I + 1 2 D∇ 2 I. Setting ∂ t I = 0 to find the steady-state solution, we arrive at a similar equation to Eq. (A7), but with the term on the left-hand side arising from diffusive transport We note here, that we are implicitly assuming a slow distribution of initial targets, however the generalization to a finite-speed distribution is obvious from the above discussion. If we impose boundary conditions at infinity in analogy with electrostatics then we can make use of the "point-charge" Green's function 2 . Let us multiply through both side by |v| (the velocity of the RBDM) and |k | 2 such that the intensity of RBDM at the location of the earth is given by where in the second line we have assumed that I CR (y, p) = N CR (y)I voy (p) with the subscript "voy" referring to the spectrum measured by Voyager I outside the heliosphere [51]. Notice that the second quantity in curly braces has the same units as a column density and can be interpreted as an "effective column desnity" that accounts for the enhancement of the density for a system undergoing diffusive transport relative to ballistic transport.
where we have implicitly assumed a magnetic rigidity of 1 Only for q/m ≤ 10 −12 GeV −1 do mCPs undergo ballistic rather than diffusive transport [38]. 2 For the Greens function appropriate for boundary conditions on a cylindrical boundary (as is often assumed in the disk diffusion model) see §III of [78] .
1 GeV, and used GeV (with R GeV a shorthand for R 1 GeV ). The magnetic retention factor is defined as the ratio of the flux in the diffusive regime relative to the naive ballistic-transport estimate involving a line-of-sight column density Comparing to Eq. (1), we can define the leaky-box size, in a specific model of our target and cosmic-ray distribution as If we consider the ISM production in Scenario (a), then we can express the target number density as n T = N T ×(1 cm −3 ), and re-express the pre-factor in units of kpc such that where N CR and N T are dimensionless functions defined as the ratio of the local cosmic ray intensity to the voyager measurement, and the ratio of the local target density. Notice that, as defined in Eq. (A15), LB can depend on the target under consideration. For instance, if mDM has a very different spatial distribution than the interstellar gas, then a different choice of LB will be required to describe the magnetic enhancement of the ISM and mDM components of the flux as illustrated explicitly in Table I. As discussed in Section III C, the flux of mDM arriving can be sourced by mechanisms other than Rutherford scattering, however the spatial distribution of mDM alone can also lead to an enhanced flux. For example, if there are overlapping regions where both the CR intensity and the mDM density are simultaneously overdense, then this can lead to a larger LB as defined in Eq. (A15). Similarly, if the low-energy CR intensity is larger than what Voyager has measured (as was expected in some models), then this could also enhance the low-β χ part of the RBDM flux (which dominates the event rate estimates at XENON1T).
To make things concrete we calculate the integral in Section III C for three different profiles of N T in Eq. (A15) with N CR held fixed. These choices are meant to be illustrative rather than realistic, and should only be interpreted as a useful exercise for gaining intuition. Nevertheless they give some sense of a credible interval of values that LB can assume.
1. The mDM's number distribution follows a spherically NFW profile.
2. The mDM is coupled to baryons and follows the baryonic density exactly; all of mDM is located inside the disk.
3. The mDM has been evacuated into two diskshaped pancakes each containing half the mDM population, with both pancakes displaced some ∼ 5 kpc perpendicular to the disk.
In each of these scenarios we can calculate the intensity of mDM arriving at Earth as predicted by Eq. (A11) for a rigidity of R = 1 GeV, and then infer the equivalent leaky-box length used to parameterize the results in the main text.
For the cosmic ray density profile we have produced a (very rough) interpolation of Fig. 11 of [63]. The explicit functional form is given by which has been written in cylindrical coordinates with r and z expressed in kpc. We take the Earth to be located at r ⊕ = 8 kpc, and z ⊕ = 0. We perform the same exercise for Scenario (a) using a disk-like exponential profile for the ISM density (A17) The normalization is chosen such that the ISM density at Earth's location is 1 cm −3 . Using this profile, we estimate that the effective column density is 1.14 × 10 5 kpc/cm 3 , corresponding to a leaky-boxy length of 19 kpc.
Appendix B: Adjusting the free electron cross section to describe atomic ionization In this Appendix we describe our approach to the atomic ionization by mCPs. The Born approximation for our problem remains valid at all times due to the smallness of , or more precisely due to α/β χ 1 condition. The most precise prescription for treating atomic ionization must result from atomic theory calculations and/or experimental measurements using ordinary Standard Model projectiles. While the former is computationally challenging for the problem at hand, the latter often exists, albeit in domains of energy/momentum transfers not directly applicable to our case. When the energy transfer to the electron much exceeds its binding energy, the scattering is quasi-free. Nevertheless, one can consider approximate treatments of the atomic physics for the needs of direct detection problem (see e.g. [79]). These include i) the replacement of the final state electron wave function by plane waves exp(ikr), as well as ii) the prescription where matrix elements are calculated using final state plane waves, but the cross section is additionally augmented by the so-called Fermi (or Sommerfeld) factor: where for the inner shells α → Z eff α substitution is Sandwich profile 9.7 × 10 6 16.
TABLE I. Effective column density for a mCP with a rigidity of 1 GeV for the three scenarios outlined in the text. The sandwich profile is composed of two disk profiles, each with half of the mDM, located 5 kpc away from the galactic disk (out of plane). The NFW profile is normalized such that the local mass density of mDM equals fχ × ρDM. The disk profile assumes nmDM(r, z) ∝ exp[−|z|/(0.7 kpc − r/(3 kpc)], and is normalized such that the mass of the disk is equal to fχ times the mass of the Milk Way's dark matter population.
used. The reliability of these approximate methods for the xenon atom is questionable. In addition, if the plane waves are used in the final state that are not orthogonal to the bound electron wave function, an adhoc substitution in the transition operator, exp(iqr) → exp(iqr) − 1, is often used to restore the correct qualitative behavior of the transitional matrix element in the small momentum transfer q limit. In order to gain insight into this problem, we conduct the comparison of these methods in a case where there exists closed analytic formulae that takes into account exact outgoing wave functions, i.e. the hydrogen atom. The full expression for the ionization rate can be found in e.g. the textbook [80], section 148, problem 4. Figure 12 shows the results for the hydrogen atom, as a function of the momentum of the ionization electron, using atomic units with a B = 1 αme = 1; m e α 2 = 1, so that the energy of the final state electron is k 2 /2 and the total energy absorbed by hydrogen is (1 + k 2 )/2. These cross sections, dσ/dk, are integrated over all possible momentum transfers q and are given in units of σ 0 = 8π(m e α) −2 α vχ 2 , where v χ is the velocity of passing mCP particles. We show the full cross section (first order Born approximation with respect to the incoming particle, with exact accounting of the ionization electron final state wave functions), as well as the results of the approximations i and ii. (In both approximations, we use the already-mentioned "subtract 1" prescription.) We see that in both cases the approximate cross section represents a significant overestimation in the k ∼ O(1) region. In particular, the prescription ii overestimates the peak cross section by over one order of magnitude. At the same time, we also plot the result for the cross section on free and static electrons that has dσ/dk ∝ 1/k 3 scaling, which gives a good fit everywhere for k > 1. In order to remove what would be an unphysical in atomic setting 1/k 3 infrared tail, we replace the free cross section with the constant value at 0 < k < 1. The resulting curve deviates from the full answer by at most a factor of 2, which we deem to be acceptable on the scale of other uncertainties involved elsewhere in our calculations. Therefore, we resort to employing the approximate treatment of the ionization, and for that purpose modify the cross section of mCP scattering on free electrons using the following prescription. Consider the ionization of an electron with binding energy E b , resulting in the energy transfer to the atom E t . When the kinetic energy of the final state electron becomes equal or larger than the initial binding energy, we are using scattering on free electrons. For the kinetic energy of final state electron 0 ≤ E t − E b ≤ E b we approximate the cross sections by a constant, and for E t < E b the cross section is assumed to be zero: These formulae are written in the limit of m χ m e and m χ β 2 χ /2 E t . The total ionization cross section at a given energy transfer E t is obtained by summing over all shells: where N b = 2J b + 1 is a given shell electron multiplicity. The values of energies for xenon shells are taken from Ref. [81]. When m χ β 2 χ /2 E t is not assumed, we need to include an additional multiplicative factor 1 − β 2 in the third line of Appendix B, and 1 − β 2 in the second line for continuity, where E max is the maximum energy that kinematically can be passed to a static electron by a passing particle of mass m χ and velocity β χ .
Given the log-log nature of the mCP plots explored in this paper, we believe that the accuracy provided by the above treatment is sufficient. If, however, an mCP flux, or any other light DM flux, is entertained to be behind recent XENON1T excess, a far more thorough atomic physics treatment is needed (possibly with multielectron correlations taken into account), as the exact status of the excess vs exclusion limits may indeed be quite sensitive to atomic details.