Earth-bound Milli-charge Relics

Dark sector particles with small electric charge, or millicharge, (mCPs) may lead to a variety of diverse phenomena in particle physics, astrophysics and cosmology. Assuming their possible existence, we investigate the accumulation and propagation of mCPs in matter, specifically inside the Earth. Even small values of millicharge lead to sizeable scattering cross sections on atoms, resulting in complete thermalization, and as a consequence, considerable build-up of number densities of mCPs, especially for the values of masses of GeV and higher when the evaporation becomes inhibited. Enhancement of mCP densities compared to their galactic abundance, that can be as big as $10^{14}$, leads to the possibility of new experimental probes for this model. The annihilation of pairs of mCPs will result in new signatures for the large volume detectors (such as Super-Kamiokande). Formation of bound states of negatively charged mCPs with nuclei can be observed by direct dark matter detection experiments. A unique probe of mCP can be developed using underground electrostatic accelerators that can directly accelerate mCPs above the experimental thresholds of direct dark matter detection experiments.


I. INTRODUCTION
Charge quantization is a century old mystery. While explanations for quantization exist, the resultant predictions of magnetic monopoles and/or manifestation of grand unification (GUT) have not been observed despite systematic efforts. This has led to the more open-minded approach to charge quantization, and exploration of the possible existence of non-quantized charges also referred to as milli-charge particles (mCPs). In recent years mCPs have received further theoretical and experimental scrutiny (see e.g. a selection of papers on theoretical and experimental efforts: [1][2][3][4][5][6][7][8][9][10][11]).
On the theoretical side, models with pure mCPs as well as models where smallness of effective electric charge is achieved via photon mixing with a new nearly massless gauge boson have been considered [12]. Since their stability is guaranteed by their U (1) charge, a non-trivial relic abundance surviving from the Big Bang can be expected. Depending on their mass and charge, they could explain all or part of the observed dark matter, called milli-charge dark matter or mCDM with their abundance set by the freeze-out or freeze-in. (Freeze-out refers to the self-depletion through annihilation from the initially fully thermally excited abundance, while the freeze-in is a sub-Hubble-rate-induced population corresponding to smaller couplings.) Regardless of cosmological abundance of mCPs, there exists a smaller yet irreducible abundance arising from the interaction of cosmic-rays with intervening matter [7,10].
Owing to the enhancement of mCP scattering cross sections at low momentum transfer, they have been invoked recently as an explanation of certain low-energy * hramani@stanford.edu anomalies, such as enhanced absorption of CMB by 21cm absorbers [13][14][15], and excess of the keV scale ionization in the Xenon 1T experimental results [10,16,17].
Regardless of possible anomalous results explained by mCPs, there have been a plethora of efforts looking for mCPs in collider and beam-dump experiments, that should be viewed in a broader context of exploring the dark sectors [18]. mCP relics depending on their speed could also be detected in dark matter direct detection and neutrino experiments. In addition, there are strong limits on mass vs coupling parameter space arising from cosmology [19][20][21] and galactic astrophysics [22,23] as well as from stellar energy losses [24,25].
Despite these efforts, there is a tantalizing window of parameter space that current and future experimental efforts cannot access. This window corresponds to m Q ≈ 10 MeV and heavier, where BBN bounds do not apply (notice, however that BBN can still limit such models through the excess abundance of dark photons in corresponding models, see e.g. [26]). Defining the mCP charge as e, with e the electric charge, 0.1 are not directly limited by collider and beam dump experiments, for m Q ≈ 1 GeV and heavier. If these mCPs make up a fraction f Q of the DM, then for large enough charge, the atmospheric or rock overburden is enough to slow them down to small values of kinetic energies and making them inaccessible to current direct detection (DD) experiments.
In this paper, the main point to be investigated and exploited for possible novel signatures is mCPs slowing down inside the Earth, resulting in a dramatic increase of their number densities at the locations of underground laboratories. This paves the way to novel methods of searching for mCP that we also explore in this paper. A direct consequence of the mCP's precipitous slow-down is that this mCP thermalizes with the atmosphere (earth) and for large enough m Q , does not possess a large enough velocity to escape the planet subsequently. Barring subsequent evaporation, this builds up through the age of the Earth t ⊕ ∼ O(4 × 10 9 ) years leading to terrestrial densities of mCPs several orders of magnitude larger than the virial density of weakly-interacting DM. If the incoming mCP flux makes up a fraction f Q of the incoming DM flux, then terrestrial densities as high as n terr Q ≈ f Q 10 14 cm 3 can be obtained. Depending on the precise value of the mass, this tremendous density increase may be concentrated inside the Earth's core, or be spread out through the whole Earth's volume. Even in the case of heavy masses, the constant vertical downward drift of thermalized mCPs is slow, leading to the "traffic jam" effect that we have discussed earlier [27][28][29].
Previous literature has explored the build up of this large density of DM that has large cross-section with the SM nuclei [30], as well as its consequences on Earth, stars [31], comets [32] and even exoplanets [33]. However to our knowledge, the build up of milli-charges specifically and its consequences has missed scrutiny. In this work, we explore various sources of mCP flux on Earth and the subsequent build-up terrestrially. It was shown previously that masses above a GeV sink to the center of the Earth and below a GeV evaporate away leading to a narrow window of mass where this terrestrial accumulation is relevant for experiments near the surface [30]. However, due to the massless-mediator nature of mCP-SM interactions, the mCP slow-down increases this crosssection leading to larger abundances expanding the mCP masses that accumulate appreciably near the Earth's surface.
The large enhancement in local density compared to the virial density by extremely large factors (e.g. for some parts of the parameter space by as much as 10 14 ), opens up novel detection strategies. In this paper, we discuss only a small subset of possible new phenomena and strategies for mCP detection: -Terrestrial mCPs with relatively large charge can bind with SM nuclei which can be looked for in exotic isotope searches. If the binding energy is large enough, there is a chance of seeing negative mCP -atomic nucleus "recombination" inside a detector (cf. Refs. [34][35][36] for analogous ideas).
-For small enough , where binding is not allowed, annihilations of mCPs with anti-particles into SM can be looked for in the large volume neutrino detectors, and specifically in super-Kamiokande.
-Milli-charged particles that have accumulated inside electrostatic accelerators can be "accidentally" accelerated (for large enough Q) and gain energy. The subsequent scattering in the low-threshold direct detection experiments is capable of providing very strong sensitivity to the mCPs. This is especially relevant in light of new efforts to install MeVscale accelerators (for the studies of rare nuclear reactions) in the underground laboratories [37].
The last idea on this list is somewhat reminiscent of the proposal [38] that seeks to perturb the flow of mCPs by EM fields, with subsequent detection of this perturbation in the adjacent spatial region. The proposal of Ref. [38] is aimed at smaller mass and smaller compared to those explored in this paper.
The rest of this work is organized as follows. In Section. II, we explore the capture of virial mCDM, subsequent evaporation and the resultant density near the surface. In Section. III, sources of high rigidity mCP fluxes are explored and their resultant density near the surface is calculated. This is followed by Section. IV which deals with bound state formation between mCPs and SM particles. Section. V deals with the detectable consequences of mCPs annihilating inside volumes of neutrino detectors. Section. VI provides a novel proposal to detect mCPs in an electrostatic accelerator. We present concluding remarks in Section. VII. The solid line corresponds to the volume averaged equilibrium number density n Q . Due the Earth's gravity, larger masses sink faster into the earth and only a small fraction of n Q is accessible near the surface. The dashed line corresponds to the local number density n loc for large enough millicharge .

II. TOP-DOWN ACCUMULATION
The goal of this subsection is to consider, in broad strokes, the accumulation of virial mCDM with large enough charge such that thermalization on Earth is rapid. We only consider asymmetric DM in this subsection, such that annihilations can be ignored. The symmetric case where annihilations are relevant for resultant terrestrial densities, is considered in Section. V.
In the literature, two types of mCPs have been considered. The first type is minimal mCPs, truly millicharged under the SM U (1), without photon-dark-photon mixing. These particles have properties identical to SM charge at all length scales. This property results in highly complex dynamics for the mCP propagation through the galaxy, through the solar system and on Earth. This is because these mCPs interact with the galactic magnetic field, with the solar wind and finally with the electric field between the ionosphere and ground. We leave this problem for future work.
Here, we instead consider particles Q that are charged under a dark U(1) with charge g Q , with the dark photon kinetically mixing with the SM photon with mixing parameter κ. At energy scales ω m A , the mass of this dark photon, these particles Q act effectively as mCPs with charge e = κg Q , where e is the electron charge. A large enough m A can be chosen so as to turn off the long range effects described above and yet keep the milli-charge properties which lead to testable consequences that we highlight below. It is likely that such a set-up would lead to tensions with cosmology via the increase of N ef f [26], which perhaps could be circumvented at the cost of adding additional ingredients to the evolution of the primordial universe.

A. Capture and evaporation
For the couplings we are interested in, all of the mCDM gets captured, or in other words, it is in the multiple scattering regime, σ × (n atom ) 1, and will lose its initial kinetic energy down to characteristic thermal energy. The Earth volume averaged number density of dark matter captured is Here f Q is the fraction of virialized DM in milli-charges, defined as where ρ DM is the local dark matter density, ρ DM 0.3 GeV/cm 3 . v vir refers to an average velocity of galactic mCPs. Eqn (1) neglects gravitational focussing, and in the case of the Earth's capture it is well-justified.
For smaller masses, dark matter that thermalizes with atmosphere (water, rock etc) has a thermal velocity v th > v esc . Therefore, there exists a "last scattering surface" somewhere in the atmosphere, from which the most velocitized mCPs can freely escape, i.e. evaporate. Adopting earlier results, see e.g. [30], the evaporation rate per one mCP particle can be estimated as, The equilibrium density on Earth is given by, This is plotted as the solid curve in Fig. 1. One can see that above 1 GeV the evaporation is no longer a factor, and the captured number density displays the familiar m −1 Q scaling. It is evident from this plot that the the accumulated density through the lifetime of the Earth can be up to fifteen orders of magnitude larger than the galactic density of the mCPs.
At the same time, the regime of light mCPs experiencing strong evaporation is more difficult to analyze. In particular, slow-down in the upper atmosphere and subsequent diffusion and evaporation can be altered by many effects including the macroscopic mass transport. Precise analysis of this regime (e.g. m Q 10 − 500 MeV) goes beyond the scope if this work and is excluded from Fig. 1.

B. Density near the surface
While the number density averaged over the Earth volume is given in Eqn. 4, the equilibrium density profile as a function of depth depends on the mass of the mCDM. The presence of gravity as well as pressure and temperature gradients results in rearrangement with a density profile that is mass dependent. This profile was evaluated in [30] for the resultant stable population of strongly interacting particles. The main conclusions of [30] were that, for relevant cross-sections, the number density at the surface n jeans ≈ n Q at m Q ≤ 1 GeV, while there was diminishing number density of dark matter near the surface for m Q ≥ 1 GeV owing to sinking to a greater depth.
However this sinking is not immediate. Diffusion rates and terminal velocities determine the net sinking of heavier dark matter to lower altitudes. To estimate these rates we need the transfer cross section in terrestrial medium (which we will call "rock"). The transfer crosssection σ T for thermalized dark matter with atoms is estimated in Appendix A. In the perturbative regime, one can get good estimates with simple models of the charge distribution inside an atom. We find that to a very good approximation for both attractive and repulsive interactions, In this formula, µ 2 rock,Q stands for the reduced mass of an atom-mCP system, Z is the atomic number and v th is the typical thermal velocity of a particle with mass µ .
While all of the dark matter is captured in the atmosphere or close to the solid Earth's surface, the random walk due to thermal motion will cause DM to diffuse deeper into rock. The time taken to diffuse to depth h is given by, assuming simple Brownian motion. Here n rock is a typical number density of atoms. If t diff ≤ t ⊕ , the age of the Earth, this causes the DM to spread out all over the Earth volume and reach the average densities given in Eqn. 4. However, in the presence of gravity there is a mechanism for vertical sinking of DM, through the gravitational pull. Under the gravitational interactions mCPs can acquire terminal velocity that can be estimated as [39], While this causes suppression of DM densities near the surface relative to the volume averaged number density in Eqn. 4, there is still a traffic jam effect on the way down as pointed out in [27] in the limit v term v vir . The enhancement in number density occurs because of flux conservation, as smaller velocity at a depth h translate to larger number densities: Note that this enhancement should never exceed the volume average in Eqn. 4, after all, diffusion will always be more efficient. This enhancement factor is the most relevant for the heavy dark matter that has its equilibrium position much closer to the center of the Earth than to its surface. However, for the GeV-scale dark matter, this type of "traffic jam" enhancement can be still small compared to equilibrium Jeans-type contribution. Putting these two enhanced populations together, finally, the number density for m Q ≥ 1 GeV at an underground laboratory can be estimated as, while for m Q < 1 GeV it is simply given by n jeans calculated in [30]. The dashed line in Fig. 1 illustrates this total enhanced density as a function of DM mass m Q . This is applicable only for the millicharge large enough such that it stops in the corresponding overburden. Thus, this line is applicable only for For such large , the transfer cross-section is given by the second term in Eqn. 5 which is independent of and hence this line is applicable for all in this range. One notices three distinct regimes. For m Q 4 GeV, the density is dominated by n jeans . For larger masses, there is sinking, and the number density is instead given by n tj . The kink at around 50 GeV, corresponds to shift from m Q ≤ m rock to m Q > m rock in Eqn. 7.

III. BOTTOM-UP ACCUMULATION
The previous section dealt with top-down accumulation; mCPs rapidly thermalizing in the overburden followed by diffusion/gravity populating lower altitudes. However, mCPs with large enough rigidity ( momentum charge ) could penetrate the overburden and get deep into the Earth before thermalizing. They then diffuse through rock and the atmosphere before finally evaporating. This diffusion time effectively acts as the time of accumulation of these mCPs leading to moderate local density. This large rigidity could arise either due to mCP possessing large momenta or small charge. An irreducible source of a fast flux occurs due to cosmic ray produced mesons which decay into mCPs. The flux for such mCPs was treated in detail in [7,10], and we do not repeat it here. An alternative source of the fast flux could be the cosmic ray collisions with mCDM that accelerates mCDM particles to higher velocity via Rutherford scattering [10]. Finally, virial mCPs with small enough charge could also penetrate the overburden and diffuse subsequently.

A. Fast Flux
We next estimate the accumulation of mCPs due to the atmospheric fast flux. In the absence of evaporation, and in the assumption that all mCPs generated in the atmosphere are captured and retained, we would have, where dΦ d(βγ) is the incoming mCP flux per interval of βγ. An mCP with mass m Q , charge and boost factor βγ penetrates a distance d pen (m Q , , βγ) in the rock, that we estimate using the Bethe-Bloch formula. We cut off the integral approximately at βγ max , above which particles cannot be stopped by the entire column density of the Earth. βγ max is given by equating the penetration depth to the Earth's radius: d pen (m Q , , βγ max ) = R ⊕ . We take the estimate for the atmospheric flux from [7]. This quantity n loc is plotted in Fig. 2 (left panel) as a function of and m Q . Thus, this plot gives an expected average density created by cosmic rays, if evaporation can be neglected. Also shown are a compilation of existing constraints on mCPs from cosmic rays [7], beam dumps and colliders [11]. N eff constraints from [20] are displayed as well. We also show the region that is accessible to surface and deep underground direct detection  Left: number densities neglecting evaporation; Right: realistic number densities upon accounting for evaporation. Also shown are a compilation of existing constraints on mCPs from cosmic rays [7], beam dumps and colliders [11]. N eff constraints from [20] are displayed as well. We also show the region that is accessible to surface and deep underground direct detection (DD) [40]  Contours of n Q /f Q arising from accumulation due to virial mCDM density are plotted.
(DD) [40]. We find that for 5 × 10 −6 , there is negligible terrestrial accumulation since the mCP interacts feebly enough to penetrate the entire Earth without thermalization. As is increased, there is also a larger flux due to preferential meson decays resulting in larger accumulation. We find that densities up to n Q ≈ 1cm −3 can be achieved barring evaporation.
It is clear that this density will be diminished due to evaporation, and the total local density will depend sen-sitively on the retention time. This can be thought as the time taken for the mCP to diffuse out to the surface (with subsequent evaporation determined by m Q ) is given by the diffusion time t diff (d pen ) given in Eqn. 6. We approximate the total number of mCPs collected in the infinitesimal shell with depth d pen to have been distributed with linearly decreasing density in the shell of thickness d pen . Thus we have for the local density, This quantity is plotted in Fig. 2 (right panel). The effect of evaporation is severe for lighter masses, due to their su-perior thermal velocities which leads to shorter diffusion times. Above a GeV, evaporation is negligible and the left and right panels present near identical densities. In the region currently allowed by terrestrial bounds, densities upto n Q ≈ 10 −4 cm 3 can be achieved. While this is several orders of magnitude smaller than the densities found in Section II for mCDM, it is important to note that this is an irreducible density with no assumptions regarding the relic density of these mCPs.

B. Virial mCDM with Small
Alternatively, virial mCPs with lower charge could also reach significant depths before thermalizing. This thermalized mCP then diffuses outwards before eventually evaporating. Unlike the model variation presented in Section II, these mCPs can reach surface and underground detectors without significant slow-down, leading to strong limits from existing DD experiments shown in Fig. 3, left panel, which we rescaled from the f Q = 1 constraints shown in Ref. [40]. However, as seen in the figure, the limits relax for subcomponent mCDM, with no existing limits below f Q ≈ 10 −8 .
We next calculate the local thermalized density to explore a complementary probe of this parameter space.
The local density can be calculated using Eqn. 12 with the flux given by, where the non-relativistic approximation βγ ≈ v Q is used, and g(v Q ) is the Maxwell Boltzmann distribution boosted to the Earth frame given in Appendix B of [41]. We provide contours of the resultant mCDM density at 1km depth in rock in Fig. 3 (Right panel). For small enough the entire Earth is not enough to stop virial DM and hence there is no accumulation. Above the black dashed line, mCDM stops in the overburden and we leave the estimation of the number densities for top-down accumulation to future work.
We note that for part of this parameter space (e.g. for ∼ 10 −6 , m Q ∼ 100 MeV − 1 GeV and f Q ∼ 10 −8 concentrations of mCPs in the underground laboratories can reach 10 2 cm −3 , which may still provide some basis for future detection, despite the smallness of , as discussed in Section VI.

IV. BOUND STATES
The thermalized dark matter now has kinetic energy on the order of kT room ∼ 0.025 eV. The negatively charged mCDM can now form bound states with atoms. At a sufficiently large mass, the lowest orbit of such bound states is inside the atomic K-shell if, Here µ Q,N is the reduced mass of the mCDM-nuclear system, a 0 is the Bohr radius, and Z is the atomic number of the nucleus. If this condition is satisfied, the binding energy will be on the order of E B Z 2 2 µ Q,N /2 13.6 eV × Z 2 2 (µ Q,N /m e ).
One may worry that if the bound states form, the existing atomic electron gains some positive energy due to effective screening of the atomic nucleus, Z → Z − . The total binding energy of an atom scales as 16eVZ 7 3 . This is obtained by observing that to a good approximation, the electrons Z e = Z in number at a distance a 0 Z Requiring E B > 0 gives This is always weaker than Eqn. (14) for Z > 1, and therefore we take (14) as the main criterion for the bound state formation by negatively charged mCPs. These exotic bound objects have net charge . The nucleus could attract more negative charged mCPs but it is unlikely because there are fewer mCPs than nuclei. They could form bound states with positively charged mCPs as well.
The positive mCPs could form bound states with free electrons when the bound state energy exceeds their kinetic energy due to thermal equilibrium with the room temperature, or could form bound states with negatively charged mCDM / atom-negative charge hybrids. This bound state has energy, Notice, however, that the bound states of Q + Q − can occur due to A exchange, and can be significantly deeper than electrostatic bound states of the same particles. The existence of the Q + Q − binding can have profound consequences for cosmological abundances and late annihilation of mCP, as is discussed e.g. in Refs. [42][43][44][45].
In both Eqs. (17) and (18), we have required that the depth of these bound states exceed typical thermal energy, because otherwise they can be easily broken up by thermal collisions. An approximate position of the critical dividing lines for the bound states is shown in Fig. 4. The region above the lines correspond allow stable bound states at 300 K. As is easy to understand, the electrostatic attraction between Q − and a nucleus is strongest, and the solid black line is plotted for a typical nucleus with Z ∼ 30.
We briefly outline observable physics effects that can occur due to a formation of (Q − N + ) bound states. The cross section leading to these bound states is not necessarily small. A typical formation of this bound state will occur with an Auger-type ejection of an electron, and subsequent cascade of the bound state down to its ground state, The rate for such a process can be large, and we make a crude estimate of the cross section by accounting for the relatively small size of the nucleus-Q − bound state, ∼ π(a 0 /Z) 2 , the probability of an outer electron to be within that distance from the nucleus ∼ Z −2 . The cross section will contain 1/v Q , the inverse velocity of the incoming particle, which will be made dimensionless by the typical velocity of an electron inside the K-shell, ∼ Zα. This way, one get the following estimate for the cross section of bound state formation, This size of the cross section will ensure relatively rapid capture of Q − , if the bound state formation is possible. The refinement of this estimate along the line of computations performed in [34,35] is possible. The capture of Q − by the nuclei of light elements may lead to exotic concentrations of (HQ), (CQ) and (OQ). However, the search techniques that involve ionization and mass spectrometry [46,47], as well as "alternative" chemical history for the milli-charged bound states poses certain difficulties in applying such bounds.
A less uncertain approach to search for an mCP "recombination" with an atom would consist in searching for heat/ionization provided by the process (19). For example, just below the (Q − N + ) boundary in Fig. 4 the recombination with light elements is not possible but recombination with atoms such as Xe or I happens readily. An ideal setup for such a probe would be the DM-Ice experiment [48] that utilizes NaI crystals shielded by ∼2.5 km of ice. Assuming the range of parameters that does not allow the formation of bound states with H, O and elements in the atmosphere, one could still expect -for a right range of {m Q , }, the exothermic reaction of Q − association with iodine atoms. Taking into account that the capture cross sections can be significant (20), all (or nearly all) of the negative mCP incident on NaI crystal may undergo the capture process. In this case one should expect that the counting rate is If the binding energy is between few keV to a 100 keV, one should compare it with the counting rates observed by these experiments that do not exceed O(10) kg −1 day −1 keV −1 , which for 10 kg crystals and two decades in energy does not exceed the total counting rate of 0.1 Hz. Therefore, for 100 GeV particles it translates to sensitivity to n Q at the level of 10 −7 cm −3 , and given results of Fig. 1, to f Q as small as 10 −19 . Conversely, one can achieve some sensitivity to cosmic ray generated mCP flux. Further gains in sensitivity can be achieved by exploiting that one and the same amount of energy is released in the formation of the bound states with a given atom, which would then show as an "unidentified" line in the spectrum taken by a detector. We finish this subsection by acknowledging the fact that the full exploration of the sensitivity to bound states throughout {m Q , } parameter space is difficult, as the binding to lighter elements will change patterns of mCP accumulation and distribution with depth. We leave more detailed exploration of the bound state related observables to forthcoming work.

V. ANNIHILATION INSIDE LARGE VOLUME DETECTORS
In this section we explore the possibility of mCDM being in equal amounts of particle anti-particle pairs. The presence of a large number of positive and negative mCDM terrestrially can lead to annihilations. Notice that the negative charge can be bound deep inside atoms, and the probability of annihilation may get significantly reduced due to the electrostatic repulsion. (A positively charged mCP would not be able to approach the orbit of bound negatively charged mCP). On the other hand, the A -induced attraction between the two mCP particles may actually overcome the Coulomb repulsion, and the annihilation may proceed even if the negative mCP is locked inside an atomic bound state. Unfortunately, reliably predicting the abundance of mCP at locations of underground laboratories when negative mCPs are intercepted by atoms is extremely challenging. For this reason, we concentrate on the region of parameter space where both the positive and negative mCDM are unbound by atoms. Annihilation of mCPs can occur via a variety of different channels. The largest cross section presumably occurs due to annihilation to two dark photons, QQ → A A . The cross-section for this process, including Sommerfeld enhancement for small velocities is [42], Here, α D = g 2 D 4π and the superscript D denotes dark final state to differentiate from annihilations to SM which we will look at shortly. v Q is the velocity of each particle which we will take to be the thermal velocity v th . This large annihilation rate can result in depletion of the large densities calculated in Section. II. The Jean's density n jeans calculated in [30] is a result of accumulation through the age of the Earth, and we find that with annihilations turned on, it is Here → 1 as expected. In the opposite limit, i.e. for large cross-sections, tanh(ζ) ζ → 1 ζ . The traffic jam contribution takes on average a time duration, L ob v sink t odot before reaching the detector and hence the suppression is smaller. Here L ob ≈ 1 km, is the length of the overburden. We conservatively estimate the traffic jam density to be i.e. 0.2 of the density that corresponds to the DM column density required to annihilate all of the incoming thermalized DM. Finally, we use n ann loc = Max n ann jeans , Min n ann tj , n tj , n Q Thus, the local density of mCPs in the lab is critically dependent on the dark fine structure α D , with larger α D leading to smaller local abundances. We next turn to observables produced by Q + Q − annihilation. While annihilation to two A is dominant, and subsequent conversion/oscillation of A into A can occur, it is severely suppressed at small mass of A . Hence we look at annihilation into SM in this work.
Visible annihilation channels via an s-channel virtual photon, QQ → A * → SM, will occur at a rate suppressed by but will result in immediate release of energy in the form of SM charged particles.
The annihilation cross-section of QQ into e + e − is given by, Next, let us estimate the number of annihilation events inside a volume that we will take to correspond to the fiducial volume V of the Super-Kamiokande experiment. The event rate is given by, (n ann loc ) 2 σ SM ann v th V = 400 10 −5 which potentially may result in a very strong sensitivity to × n Q m −1 Q . The quantity between the {} provides the additional Sommerfeld enhancement. Sensitivity to α D enters through this enhancement factor as well as through n ann loc in Eqn. 25. Fig. 5 illustrates this through limit contours for the different input values of f Q . We take O(200) events per year to be roughly a limiting count rate. The left panel corresponds to κ = 0.1 (small α D ) and the right panel corresponds to κ = 10 −4 (large α D ). For a fixed small , larger α D results in larger Sommerfeld enhancement and thus stricter limits on f Q . However, sensitivity to large is lost for large α D because this results in high annihilation rates as well, and hence extremely small local number densities n ann loc . Nonetheless new parameter space is probed for both the small α D and large α D regimes.

VI. ELECTROSTATIC ACCELERATORS
Self-annihilation, enhanced by terrestrial accumulation, provides non-trivial limits on the parameter space of the mCDM model. At the same time, quadratic scaling with abundance does not allow to probe very small f Q . The smallest f Q where experiments like Super-K will have sensitivity is for f Q ∼ 10 −7 . In this section we propose a novel strategy to test even smaller densities.
As alluded to in the Introduction, the local density of mCPs could be accelerated in a large electric field and the accelerated mCPs could then be detected. Owing to disparate charge to mass ratio compared to SM particles, oscillating field accelerators will not be suitable for mCPs. Instead electrostatic accelerators such as Van de Graaf generators and Cockcroft-Walton accelerators would be suitable. Modern accelerators with potential difference (∆V ) in the megavolt range are used in nuclear physics experiments. Examples among these are LUNA (∆V = 3.5 MV) [49], JUNA (∆V = 0.4 MV) [50] and CASPAR (∆V = 1.1 MV) [51].
Since we are considering mCPs with a dark photon, we need to first determine the plasma masses so as to ensure that the electric field is not shielded by mCPs. The plasma mass of the dark photon A in the presence of a number density n Q of mCPs Q is given by, Here g Q → 1 is set, to be conservative. The range of the electric field is thus, In other words, the screening length is larger than 1 meter for n Q 10 8 cm 3 . If the concentrations of mCPs exceed this level, it is likely that even the acceleration of "normal" protons will get compromised.
We consider the accelerator field to be turned on, but with the proton source inside the accelerator being "off". While the mCPs outside the region with the electric field receive no net acceleration, mCP particles on the inside may get accelerated. Given a E thr required for detection, it is clear that in order to have sensitivity one should require where ∆V is the accelerating voltage.
To create a more realistic description, we model the accelerator tube as a r = 1 mm radius, L = 1 meter long tube similar to the LUNA setup [49] at 1 km depth. This concept is illustrated schematically in Fig. 6. The flux of mCPs seeping into the pipe that gets accelerated to detectable energies is given by This expression is derived in detail in Appendix. B. In Fig. 7, we plot the rate for accelerated mCPs to come out of this setup. The accelerated flux of the mCPs is relatively collimated, and could be detected with relatively compact dark matter detectors. To translate the rate of accelerated mCPs to the counting rate we need to estimate the probability of generating one signal event by a single accelerated mCP particle.
We assume a dark matter detector that completely covers the geometric parameters of the mCP beam, but with an opening in the shielding that allows mCPs to enter. The probability for mCP to create a signal event can be estimated as where is taken to be ∼ 5 cm is a realistic size for a small-to-medium size DM detector, n atoms is the number density of active targets in the detector. We take n atoms 4.5 × 10 22 cm −3 to correspond to the number density of germanium atoms in Ge crystal. Finally, σ Q,atom is the cross section leading to atomic recoil or FIG. 6: Scheme of electrostatic accelerator concept. Positive mCPs between the two plates held at a potential difference ∆V get accelerated towards the more negatively charged plate and enter the low threshold detector to be detected. ionization by the accelerated mCPs. Depending on the mass of the mCP, and the type of detector, the relevant scattering is on the atomic electrons or elastic scattering on a whole atom, with subsequent ionization created in the inter-atomic collisions. In the latter case, the resulting ionization energy is typically quenched by a factor of ∼ 0.1. The transfer of energy in elastic scattering on an atom is the most efficient if m Q ∼ m atom . Using a non-pertubative estimate (Appendix A) for such a cross section, σ ∝ 4π × (µv Q ) −2 ∝ 10 −23 cm 2 for a typical mass of mCP in the 100 GeV range, and its kinetic energy in the ∼ keV range, it is easy to see that for this size of the cross section, the probability of scattering within the 5 cm detector is O(1). In addition, in order to maximize the signal from such detectors in terms of m Q , it would be advantageous to use devices with a wide range of atomic masses, such as CaWO 4 of CRESST [52].
The lowest thresholds for detection, and sensitivity to the lowest m Q can be achieved through scattering on electrons/atomic ionization. Despite the fact that even after the acceleration, the mCPs are relatively slow, one could apply perturbation theory for estimating the cross sections leading to ionization. Taking the outer shell atomic electron to be localized within space region ∼ a 0 , its interaction with an incoming mCP to scale as U ∼ α a −1 0 , perturbation theory is valid as long as U × ∆t ∼ α v −1 Q < 1 [53]. Taking m Q below a few GeV, in the 10 −6 − 10 −3 range and v Q ∼ (| e∆V |/m Q ) 1/2 ∼ 1/2 × (10 −2 − 10 −1 ), we see that the perturbativity condition is satisfied, and the cross section can be estimated to scale as πa 2 0 × (α v −1 Q ) 2 . Again, this estimate exceeds 1/(n atom ) making P ∼ O(1). Thus we conclude that the rate of accelerated mCPs in Fig. 7 can indeed translate to a similar counting rate in a small-to-medium size detector intercepting the "beam" of mCPs. Q becomes greater than 1, and the mCP-electron interaction becomes adiabatic. This regime would require additional analysis of ionization efficiency, and would also benefit from detectors that are sensitive to energy release, e.g. phonons, spread between many atoms.) In Fig.  7 (Right) we illustrate the rate sensitivity and energy thresholds required to probe mCDM with masses below 1 GeV that have accumulated via bottom-up mechanism. Given rapid advance in dark matter detectors, some part of the parameter space can be probed with existing technology. It is also clear that if at some point in the future, the detection thresholds for DM-induced recoil can be brought to a sub-eV level (see e.g. [54]), even 's as small as 10 −8 can be probed via accelerating mCPs in the underground MV voltage electrostatic accelerators.
Also of notice is potential sensitivity to mCPs via their generation by cosmic rays which results in an irreducible population on Earth. This irreducible number densities shown in Fig. 2 translate to acceleration rates plotted in Fig. 8. We find that a maximum obtainable rate can be as high as 10 3 Hz, while a lot of new parameter space can be explored with counting rates reaching down to 1 Hz and below.
In addition to novel probes discussed in this paper, we have also examined a number of other ways of constraining { , n Q } parameter space, summarized in Appendix C. None of them carry as much promise/sensitivity as the three pathways outlined in this work.

VII. CONCLUSION
We have shown that the direct probes of milli-charged particles can be advanced using their accumulation inside the Earth. Owing to a relatively large cross-sections and short free path inside dense media, the number densities of mCPs can indeed be many order of magnitude larger than the cosmological abundances. In this paper we have analyzed the main mechanisms of their accumulation, finding that when evaporation is impossible, strong enhancements of number densities are expected at the locations of the underground laboratories. The enhancement factors can approach 10 15 , and therefore even very subdominant fractions of the cosmological mCP DM can be probed. Cosmic ray induced production of mCPs creates far less abundant concentration, which nevertheless could reach up to 10 −4 cm 3 .
We pointed out three different methods that could lead to the most precise probes of mCP properties to date. Concentrations upward of 10 7 cm −3 can be probed via the annihilations of mCP particles to the SM states inside large volume detectors, such as super-Kamiokande. A minimal adjustment of already existing searches will be able to refine sensitivity in 1-to-10 GeV mass range (Fig.  5). mCPs created by cosmic rays cannot be probed this way, as the predicted densities are too small.
Another method for potentially probing even very small abundances of dark matter is via the formation of bound states of negatively charged mCPs and atomic nuclei. If {m Q , } parameters are right (i.e. just below the solid line on Fig. 4), the binding to heavy elements may be possible, while the binding to light elements is not. In this case, the dark matter experiments that use heavier elements (Xe or I) may be quite sensitive to the energy release accompanying the formation of bound states. In particular, the DM-Ice experiment that utilizes NaI shielded by Antarctic ice is a good example of a device that is extremely sensitive to such a scenario. We argue that n Q as small as 10 −7 cm −3 can be probed, that due to accumulation enhancement leads to the sensitivity to tiny f Q . (Detailed analysis of {m Q , } parameter space is left for future work, as formation of bound states throughout column density may significantly alter the predictions for the density profiles as function of depth.) We also note that for a fixed mCP mass and charge, the amount of binding energy to a given type of an atom is fixed, and therefore the formation of bound states leads to a mono-energetic energy deposition. It is well known that the recently observed Xenon 1T excess events [16] is consistent with a monochromatic signal (see e.g. [55][56][57][58]). It is then possible to speculate that the excess is maybe coming from Q − -Xe nucleus bound state formation. However, more work needs to be done to understand whether the required densities of the mCPs can be realistically expected for the environment (Gran Sasso lab) where the experiment is operating.
Finally, perhaps the most direct way of testing the mCP particles is their acceleration in underground accelerators (that see their primary application for measuring astrophysically relevant nuclear reaction cross sections). Even an "accidental" acceleration of mCPs may result in their kinetic energy going up from thermal to accelerated energies ∼ ∆V . For MV-type electrostatic accelerators, and for a generous range of , 10 −5 − 10 −1 , the resulting gained kinetic energy can be far above the thresholds of direct detection experiments at 10 eV (and possibly lower in the near future). Therefore, combining underground accelerations with the specially placed dark matter detectors along the mCP accelerated trajectory can bring significant new sensitivity, and indeed test local concentrations of mCPs down to unprecedented low values. This way, many "physics targets" can be covered. A dedicated effort in this direction has a potential to explore mCP densities created by cosmic rays, and access the region of parameter space consistent with the explanation of the EDGES anomaly.
We find that a very good approximation for both attractive and repulsive interactions, These cross sections form the basis for our code that calculates the effective slow-down, sinking velocity, diffusion coefficients, and ultimately n Q (h) in the main body of the text.

Appendix B: Accelerator geometry
The dark matter is thermal and hence has velocity The differential angular flux coming in per infinitesimal pipe length is Here θ ∈ {0, π 2 }, is the angle between the incoming particle velocity and the beam axis, while ϕ ∈ {0, π 2 } subtended by the velocity on the radial direction.
The time spent by the particles inside the pipe is, τ (θ, ϕ) = 2r v th sin θ cos ϕ (B3) The maximum time the particles can spend (because of acceleration along the beam axis) is, If a particle spends time τ = Min[τ max , τ (θ, ϕ)], inside, it is accelerated to, now, Given a threshold for subsequent detection E thr , there is an l min , l min = E thr e∆V L (B7) There are several recent experimental limits which set constraints on a local population of thermalized DM. We consider them in turn and show that none of them set constraints on terrestrial mCPs.
Cryogens Anomalous heating of cryogens was considered in [29,30,60]. However all of these experiments employ shielding to reduce black-body radiation and hence will cool the mCDM to the shield temperature before they hit the cryogen if the mean-free-path corresponding to the transfer cross-section is smaller than the thickness of the shield. This can be evaluated using λ MFP = 1 n shield σ T ≈ 10 −6 cm µ GeV (C1) Here we approximate n shield = 10 22 cm 3 and take the crosssection in Eqn. 5 to be σ T ≈ 4π . This MFP is much smaller than any realistic shielding. As a result mCDM will not be constrained by cryogens.
LHC beam lifetime In [30], limits were placed on accumulated DM with contact interactions from LHC beam lifetime, and these arguments can be generalized to other particle storage rings and accelerators. Rutherford scattering on mCPs can also take particles from the beam. This will come on top of particle loss due to the residual scattering on atomic constituents of the residual gas molecules inside the beam pipe. The scattering on mCP (that can be treated perturbatively in this problem) is at least 2 times smaller than scattering on electrons and nuclei. This way, one can estimate the "best case" sensitivity via comparing the beam losses on atoms vs mCPs n Q ∼ n atoms −2 ∼ 10 10 cm −3 × (10 −2 / ) 2 , where we took a realistic residual gas density at 10 6 particles per cm 3 . This may provide some additional sensitivity to mCP if f Q is large, but it is inferior to other probes discussed in this paper. Stability of nuclear isomers: 180 m Ta.
In [27,28] the non-observation of the decay of isomeric Tantalum was used to set limit on terrestrial DM. The form-factor to scatter with 180 m Ta naturally picks mCP masses in the hundreds of GeV and above. For positively charged mCDM, there is Coulomb repulsion with nuclei. The probability of overcoming this is given by the Gamow factor. where E g = 2µ (παZ Ta ) 2 and T is the ambient temperature. This factor evaluates to tiny values rendering the cross-section too small for [28] or any future projections of a Tantalum experiment to set relevant limits. Negatively charged mCDM of this mass in the heavy mass range can induce 180m → 180 g.s. transition. If the parameter range is such that the formation of bound states with nuclei is possible, de-excitation of 180m isotope will not provide competitive sensitivity to e.g. nucleus-mCP recombination process discussed in section IV. The main reason for that is relatively small isotopic abundance of 180m that would not be competitive with many kilograms of I or Xe employed by direct detection experiments. If is very small (and the mass is in ∼ 100 GeV range to ensure efficient de-excitation), then surface direct detection experiments will be more sensitive to Q − than Tantalum 180m.
Anomalous heat transport The heat conductivity of Earth can be modified by the presence of an Earthbound exotic species [30]. However, the heat conductivity is proportional to the mean-free-path which is extremely small in rock as explained above. As a result, there are no useful limits that we can derive from this effect.