# Earth-bound Milli-charge Relics

Maxim Pospelov  School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  William I. Fine Theoretical Physics Institute,
School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA  Harikrishnan Ramani<EMAIL_ADDRESS>Stanford Institute for
Theoretical Physics, Stanford University, Stanford, CA 94305, USA

###### Abstract

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: Haas:2014dda ;
Izaguirre:2015eya ; Magill:2018tbb ; Kelly:2018brz ; Harnik:2019zee ;
Acciarri:2019jly ; Plestid:2020kdm ; Bloch:2020uzh ; Foroughi-Abari:2020qar ;
Harnik:2020ugb ; Marocco:2020dqu ).

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 Holdom:1985ag . 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
Plestid:2020kdm ; Harnik:2020ugb .

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 anomalies, such as enhanced absorption of CMB by 21cm absorbers
Barkana:2018qrx ; Boddy:2018wzy ; Liu:2019knx , and excess of the keV scale
ionization in the Xenon 1T experimental results Aprile:2020tmw ;
Farzan:2020dds ; Harnik:2020ugb .

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
Lanfranchi:2020crw . 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 Kovetz:2018zan ; Creque-Sarbinowski:2019mcm ; Slatyer:2018aqg
and galactic astrophysics Stebbins:2019xjr ; Kachelriess:2020ams  as well as
from stellar energy losses Davidson:2000hf ; Chang:2018rso .

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}\approx 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. Vogel:2013raa ). Defining the
mCP charge as $\epsilon e$, with $e$ the electric charge, $\epsilon\lesssim
0.1$ are not directly limited by collider and beam dump experiments, for
$m_{Q}\approx 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_{\oplus}\sim O(4\times 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_{Q}^{\rm terr}\approx f_{Q}\frac{10^{14}}{\textrm{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 Pospelov:2019vuf ; Lehnert:2019tuw
; Rajendran:2020tmw .

Previous literature has explored the build up of this large density of DM that
has large cross-section with the SM nuclei Neufeld:2018slx , as well as its
consequences on Earth, stars Busoni:2017mhe , comets DeLuca:2018mzn  and even
exoplanets Leane:2020wob . 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 Neufeld:2018slx . However, due to the massless-mediator nature of
mCP-SM interactions, the mCP slow-down increases this cross-section 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. Pospelov:2008qx ; An:2012bs ;
Fornal:2020bzz  for analogous ideas).

  * -

For small enough $\epsilon$, 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 MeV-scale
accelerators (for the studies of rare nuclear reactions) in the underground
laboratories Broggini:2017cez .

The last idea on this list is somewhat reminiscent of the proposal
Berlin:2019uco  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. Berlin:2019uco  is aimed at smaller mass and smaller
$\epsilon$ 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.

##  II Top-Down Accumulation

Figure 1: Density of millicharge particles $m_{Q}$ on Earth. The solid line
corresponds to the volume averaged equilibrium number density $\langle
n_{Q}\rangle$. Due the Earth’s gravity, larger masses sink faster into the
earth and only a small fraction of $\langle n_{Q}\rangle$ is accessible near
the surface. The dashed line corresponds to the local number density $n_{\rm
loc}$ for large enough millicharge $\epsilon$.

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 milli-charged 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 $\kappa$. At energy scales $\omega\gg
m_{A^{\prime}}$, the mass of this dark photon, these particles $Q$ act
effectively as mCPs with charge $\epsilon e=\kappa g_{Q}$, where $e$ is the
electron charge. A large enough $m_{A^{\prime}}$ 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_{eff}$ Vogel:2013raa , which perhaps could be circumvented at
the cost of adding additional ingredients to the evolution of the primordial
universe.

###  II.1 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, $\sigma\times(n_{\rm
atom}\ell)\gg 1$, and will lose its initial kinetic energy down to
characteristic thermal energy. The Earth volume averaged number density of
dark matter captured is

$\displaystyle\langle n^{\rm cap}_{Q}\rangle$ $\displaystyle=\frac{\pi
R_{\oplus}^{2}v_{\rm vir}t_{\oplus}}{4/3\pi
R_{\oplus}^{3}}f_{Q}\frac{\rho_{\rm DM}}{m_{Q}}$
$\displaystyle\approx\frac{3\times
10^{15}}{\textrm{cm}^{3}}\frac{t_{\oplus}}{10^{10}\textrm{year}}f_{Q}\frac{\rm
GeV}{m_{Q}}$ (1)

Here $f_{Q}$ is the fraction of virialized DM in milli-charges, defined as

$f_{Q}=\frac{m_{Q}n_{Q}}{\rho_{\rm DM}},$ (2)

where $\rho_{\rm DM}$ is the local dark matter density, $\rho_{\rm DM}\simeq
0.3$ GeV/cm3.

$v_{\rm vir}$ refers to an average velocity of galactic mCPs. Eqn (II.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_{\rm th}>v_{\rm 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. Neufeld:2018slx , the evaporation rate per one mCP particle can be
estimated as,

$\Gamma_{\rm loss}\approx\frac{3v_{\rm
th}}{2\pi^{\frac{1}{2}}R_{\oplus}}\left(1+\frac{v_{\rm es}^{2}}{v_{\rm
th}^{2}}\right)\exp\left(-\frac{v_{\rm es}^{2}}{v_{\rm th}^{2}}\right)$ (3)

The equilibrium density on Earth is given by,

$\langle n_{Q}\rangle=\langle n^{\rm
cap}_{Q}\rangle\frac{1-\exp\left(\Gamma_{\rm
loss}t_{\oplus}\right)}{\Gamma_{\rm loss}t_{\oplus}}$ (4)

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_{Q}^{-1}$ 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}\simeq 10-500\,{\rm MeV}$) goes beyond the scope if this work and
is excluded from Fig. 1.

###  II.2 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 Neufeld:2018slx  for the resultant
stable population of strongly interacting particles. The main conclusions of
Neufeld:2018slx  were that, for relevant cross-sections, the number density at
the surface $n_{\rm jeans}\approx\langle n_{Q}\rangle$ at $m_{Q}\leq 1$ GeV,
while there was diminishing number density of dark matter near the surface for
$m_{Q}\geq 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 cross-section $\sigma_{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,

$\displaystyle\sigma_{T}$ $\displaystyle\sim\textrm{Min}\left(\frac{2\pi
Z^{2}\alpha^{2}\epsilon^{2}}{\mu_{\rm rock,Q}^{2}v_{\rm
th}^{4}},\frac{4\pi}{\mu_{\rm rock,Q}^{2}v_{\rm th}^{2}}\right)$ (5)

In this formula, $\mu_{\rm rock,Q}^{2}$ stands for the reduced mass of an
atom-mCP system, $Z$ is the atomic number and $v_{\rm th}$ is the typical
thermal velocity of a particle with mass $\mu$ .

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,

$t_{\rm diff}(h)\sim\frac{h^{2}n_{\rm rock}\sigma_{T}}{v_{\rm th}},$ (6)

assuming simple Brownian motion. Here $n_{\rm rock}$ is a typical number
density of atoms. If $t_{\rm diff}\leq t_{\oplus}$, 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 landau1981course ,

$\displaystyle v_{\rm term}$ $\displaystyle=\frac{3m_{Q}gT}{m_{\rm
rock}^{2}n_{\rm rock}\langle\sigma_{T}v_{\rm th}^{3}\rangle}\quad\quad$
$\displaystyle m_{Q}>m_{\rm rock}$ $\displaystyle=\frac{m_{Q}g}{3n_{\rm
rock}T}\left\langle\frac{v_{\rm th}}{\sigma_{T}}\right\rangle\quad\quad$
$\displaystyle m_{Q}<m_{\rm rock}$ (7)

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 Pospelov:2019vuf  in the limit $v_{\rm
term}\ll v_{\rm vir}$. The enhancement in number density occurs because of
flux conservation, as smaller velocity at a depth $h$ translate to larger
number densities:

$\eta_{\rm tj}\equiv\frac{n_{\rm tj}}{n_{\rm vir}}=\frac{v_{\rm vir}}{v_{\rm
term}}$ (8)

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}\geq 1$
GeV at an underground laboratory can be estimated as,

$n_{\rm loc}=\textrm{Max}\left(n_{\rm jeans},\textrm{Min}\left(n_{\rm
tj},\langle n_{Q}\rangle\right)\right)$ (9)

while for $m_{Q}<1$ GeV it is simply given by $n_{\rm jeans}$ calculated in
Neufeld:2018slx . 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 $\epsilon$ large enough such that it stops in the corresponding
overburden. Thus, this line is applicable only for

$\displaystyle\epsilon$ $\displaystyle\gtrsim 2\times
10^{-4}\sqrt{\frac{m_{Q}}{\rm GeV}}\quad$ surface $\displaystyle\gtrsim
3\times 10^{-6}\sqrt{\frac{m_{Q}}{\rm GeV}}\quad$ 1 km mine (10)

For such large $\epsilon$, the transfer cross-section is given by the second
term in Eqn. 5 which is independent of $\epsilon$ and hence this line is
applicable for all $\epsilon$ in this range. One notices three distinct
regimes. For $m_{Q}\lesssim 4$ GeV, the density is dominated by $n_{\rm
jeans}$. For larger masses, there is sinking, and the number density is
instead given by $n_{\rm tj}$. The kink at around 50 GeV, corresponds to shift
from $m_{Q}\leq m_{\rm rock}$ to $m_{Q}>m_{\rm rock}$ in Eqn. 7.

##  III Bottom-Up Accumulation

Figure 2: Accumulated terrestrial density of mCPs arising from decay of mesons
produced by cosmic rays in the atmosphere. 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 Plestid:2020kdm , beam dumps and colliders Marocco:2020dqu .
$N_{\rm eff}$ constraints from Creque-Sarbinowski:2019mcm  are displayed as
well. We also show the region that is accessible to surface and deep
underground direct detection (DD) Emken:2019tni

Figure 3: Left: Existing DD limits on mCDM parameter space rescaled from the
$f_{Q}=1$ constraints of Emken:2019tni , Right: Contours of $n_{Q}/f_{Q}$
arising from accumulation due to virial mCDM density are plotted.

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 ($\frac{\rm momentum}{\rm
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
Plestid:2020kdm ; Harnik:2020ugb , 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 Harnik:2020ugb . Finally, virial mCPs with small enough charge
$\epsilon$ could also penetrate the overburden and diffuse subsequently.

###  III.1 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,

$n_{\rm loc}=\int^{\beta\gamma_{\rm
max}}d(\beta\gamma)\frac{d\Phi}{d(\beta\gamma)}\frac{\pi
R_{\oplus}^{2}}{\frac{4}{3}\pi R_{\oplus}^{3}}t_{\oplus}$ (11)

where $\frac{d\Phi}{d(\beta\gamma)}$ is the incoming mCP flux per interval of
$\beta\gamma$. An mCP with mass $m_{Q}$, charge $\epsilon$ and boost factor
$\beta\gamma$ penetrates a distance $d_{\rm pen}(m_{Q},\epsilon,\beta\gamma)$
in the rock, that we estimate using the Bethe-Bloch formula. We cut off the
integral approximately at $\beta\gamma_{\rm max}$, above which particles
cannot be stopped by the entire column density of the Earth. $\beta\gamma_{\rm
max}$ is given by equating the penetration depth to the Earth’s radius:
$d_{\rm pen}(m_{Q},\epsilon,\beta\gamma_{\rm max})=R_{\oplus}$. We take the
estimate for the atmospheric flux from Plestid:2020kdm . This quantity $n_{\rm
loc}$ is plotted in Fig. 2 (left panel) as a function of $\epsilon$ 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 Plestid:2020kdm , beam dumps and
colliders Marocco:2020dqu . $N_{\rm eff}$ constraints from Creque-
Sarbinowski:2019mcm  are displayed as well. We also show the region that is
accessible to surface and deep underground direct detection (DD) Emken:2019tni
. We find that for $\epsilon\lesssim 5\times 10^{-6}$, there is negligible
terrestrial accumulation since the mCP interacts feebly enough to penetrate
the entire Earth without thermalization. As $\epsilon$ 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}\approx 1\textrm{cm}^{-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 sensitively 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_{\rm diff}(d_{\rm pen})$ given in Eqn. 6. We approximate the total number
of mCPs collected in the infinitesimal shell with depth $d_{\rm pen}$ to have
been distributed with linearly decreasing density in the shell of thickness
$d_{\rm pen}$. Thus we have for the local density,

$n_{\rm loc}(h)\approx\int^{\beta\gamma_{\rm
max}}d(\beta\gamma)\frac{d\Phi}{d(\beta\gamma)}\frac{\pi R_{\oplus}^{2}t_{\rm
diff}(d_{\rm pen})}{\frac{4}{3}\pi\left(R_{\oplus}^{3}-(R_{\oplus}-d_{\rm
pen})^{3}\right)}\frac{h}{d_{\rm pen}}\approx\int^{\beta\gamma_{\rm
max}}d(\beta\gamma)\frac{d\Phi}{d(\beta\gamma)}\frac{h}{v_{\rm th}\lambda}$
(12)

This quantity is plotted in Fig. 2 (right panel). The effect of evaporation is
severe for lighter masses, due to their superior 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}\approx\frac{10^{-4}}{\textrm{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.

###  III.2 Virial mCDM with Small $\epsilon$

Alternatively, virial mCPs with lower charge $\epsilon$ 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. Emken:2019tni . However, as seen in the figure, the limits relax
for subcomponent mCDM, with no existing limits below $f_{Q}\approx 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,

$\frac{d\Phi}{d(\beta\gamma)}\approx\frac{d\Phi}{dv_{Q}}=v_{Q}\frac{f_{Q}\rho_{\rm
DM}}{m_{Q}}g(v_{Q})$ (13)

where the non-relativistic approximation $\beta\gamma\approx v_{Q}$ is used,
and $g(v_{Q})$ is the Maxwell Boltzmann distribution boosted to the Earth
frame given in Appendix B of Essig:2015cda .

We provide contours of the resultant mCDM density at 1km depth in rock in Fig.
3 (Right panel). For small enough $\epsilon$ 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 $\epsilon\sim
10^{-6}$, $m_{Q}\sim 100\,{\rm MeV}-1\,{\rm GeV}$ and $f_{Q}\sim 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 $\epsilon$, as discussed in Section VI.

##  IV Bound states

The thermalized dark matter now has kinetic energy on the order of $kT_{\rm
room}\sim 0.025~{}\textrm{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,

$\displaystyle r$
$\displaystyle=\frac{a_{0}}{Z\epsilon}\frac{m_{e}}{\mu_{Q,N}}<\frac{a_{0}}{Z}$
$\displaystyle\epsilon$ $\displaystyle>\frac{m_{e}}{\mu_{Q,N}}$ (14)

Here $\mu_{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}\simeq Z^{2}\epsilon^{2}\mu_{Q,N}/2\simeq 13.6\,{\rm eV}\times
Z^{2}\epsilon^{2}(\mu_{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\to Z-\epsilon$. The total binding energy of an atom scales as
$16\textrm{eV}Z^{\frac{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_{N}^{-\frac{1}{3}}$ from the nucleus with effectively un-screened
charge $Z_{N}\simeq Z$. Substituting $Z_{N}\rightarrow Z_{N}-\epsilon$, the
net binding energy is

$\displaystyle\Delta
E_{B}=\left(13.6\epsilon^{2}Z^{2}\frac{\mu_{Q,N}}{m_{e}}-21.3\epsilon
Z^{\frac{4}{3}}\right)\textrm{eV}$ (15)

Requiring $E_{B}>0$ gives

$\displaystyle\epsilon>\frac{8\times 10^{-4}}{Z^{\frac{2}{3}}}\frac{\rm
GeV}{\mu_{Q,N}}.$ (16)

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 $\epsilon$. 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,

$\displaystyle E_{B}$
$\displaystyle=13.6~{}\textrm{eV}\epsilon^{2}>0.025\textrm{eV}$
$\displaystyle\implies\epsilon$ $\displaystyle>0.042.$ (17)

or could form bound states with negatively charged mCDM / atom-negative charge
hybrids. This bound state has energy,

$\displaystyle E$
$\displaystyle=13.6~{}\textrm{eV}\epsilon^{4}\frac{m_{Q}}{m_{e}}>0.025\textrm{eV}$
$\displaystyle\implies\epsilon$
$\displaystyle>0.2\left(\frac{m_{e}}{m_{Q}}\right)^{\frac{1}{4}}$ (18)

Notice, however, that the bound states of $Q^{+}Q^{-}$ can occur due to
$A^{\prime}$ 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. ArkaniHamed:2008qn ;
Pospelov:2008jd ; An:2016gad ; Cirelli:2016rnw .

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\sim 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,

${\rm Atom}+Q^{-}\to({\rm Atom}^{+n}Q^{-})+ne+(\gamma).$ (19)

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, $\sim\pi(a_{0}/Z)^{2}$, the probability of an
outer electron to be within that distance from the nucleus $\sim 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, $\sim Z\alpha$. This way, one get the following
estimate for the cross section of bound state formation,

$\sigma_{\rm capture}v_{Q}\sim\frac{\pi\alpha a_{0}^{2}\times c}{Z^{3}}\propto
10^{-23}\,{\rm cm}^{2}\times c\times(30/Z)^{3}.$ (20)

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 Pospelov:2008qx ;
An:2012bs  is possible.

The capture of $Q^{-}$ by the nuclei of light elements may lead to exotic
concentrations of $({\rm H}Q)$, $({\rm C}Q)$ and $({\rm O}Q)$. However, the
search techniques that involve ionization and mass spectrometry Smith:1982qu ;
Hemmick:1989ns , 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 deSouza:2016fxg  that utilizes NaI crystals shielded by
$\sim$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},\epsilon\\}$, 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

$\displaystyle\frac{\rm Events}{\rm time}$ $\displaystyle\sim n_{Q}(2.5\,{\rm
km})\times v_{\rm th}\times{\rm Area}$ (21) $\displaystyle\propto 10^{6}{\rm
Hz}\times\frac{n_{Q}}{1\,{\rm cm^{3}}}\times\frac{\rm Area}{100\,\rm
cm^{2}}\times\left(\frac{100\,\rm GeV}{m_{Q}}\right)^{1/2}.$

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)\,{\rm kg}^{-1}{\rm day}^{-1}{\rm 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}\,{\rm 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},\epsilon\\}$ 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.

Figure 4: Parameter space that allows for the bound states of mCPs with nuclei
(Solid), electrons (Dashed) and between themselves (Dotted) at $300$ K. The
line corresponds to the bound state energy $E_{B}=300$ K and the region above
corresponds to deeper bound states.

##  V Annihilation inside large volume detectors

Figure 5: Limits corresponding to more than 200 events in 22 kton year
exposure at Super-K for $Q\bar{Q}\rightarrow e^{+}e^{-}$ for different DM
fractions $f_{Q}$ for two different choices of the kinetic mixing parameter
$\kappa$, Left: $\kappa=0.1$ and Right: $\kappa=10^{-4}$ are shown. For the
same millicharge $\epsilon$, smaller $\kappa$ corresponds to larger
$\alpha_{D}$ and hence larger Sommerfeld enhancement for annihilation. This
leads to limits shifting downward for smaller $\kappa$.

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^{\prime}$-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, $Q\bar{Q}\to A^{\prime}A^{\prime}$. The cross-section for this
process, including Sommerfeld enhancement for small velocities is
ArkaniHamed:2008qn ,

$\sigma_{\rm
ann}^{D}v_{Q}=\frac{\pi\alpha_{D}^{2}}{m_{Q}^{2}}\frac{\pi\alpha_{D}}{v_{Q}\left(1-e^{-\pi\alpha_{D}/v_{Q}}\right)}$
(22)

Here, $\alpha_{D}=\frac{g_{D}^{2}}{4\pi}$ 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_{\rm th}$. This large annihilation rate can result in
depletion of the large densities calculated in Section. II. The Jean’s density
$n_{\rm jeans}$ calculated in Neufeld:2018slx  is a result of accumulation
through the age of the Earth, and we find that with annihilations turned on,
it is

$\displaystyle n_{\rm jeans}^{\rm ann}=\frac{\tanh(\zeta)}{\zeta}n_{\rm
jeans}$ (23)

Here $\zeta^{2}=\frac{3n_{\rm vir}v_{\rm vir}t_{\oplus}^{2}\sigma_{\rm
ann}^{D}v_{\rm th}}{4R_{E}}$. Note that for $\sigma_{\rm ann}^{D}\rightarrow
0$, $\frac{\tanh(\zeta)}{\zeta}\rightarrow 1$ as expected. In the opposite
limit, i.e. for large cross-sections,
$\frac{\tanh(\zeta)}{\zeta}\rightarrow\frac{1}{\zeta}$ .

The traffic jam contribution takes on average a time duration, $\frac{L_{\rm
ob}}{v_{\rm sink}}\ll t_{odot}$ before reaching the detector and hence the
suppression is smaller. Here $L_{\rm ob}\approx 1~{}\textrm{km}$, is the
length of the overburden. We conservatively estimate the traffic jam density
to be

$\displaystyle n_{\rm tj}^{\rm ann}=0.2\frac{v_{\rm term}}{L_{\rm
ob}\sigma_{\rm ann}^{D}v_{\rm th}}$ (24)

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_{\rm loc}^{\rm ann}=\textrm{Max}\left(n_{\rm jeans}^{\rm
ann},\textrm{Min}\left(n_{\rm tj}^{ann},n_{\rm tj},\langle
n_{Q}\rangle\right)\right)$ (25)

Thus, the local density of mCPs in the lab is critically dependent on the dark
fine structure $\alpha_{D}$, with larger $\alpha_{D}$ leading to smaller local
abundances. We next turn to observables produced by $Q^{+}Q^{-}$ annihilation.
While annihilation to two $A^{\prime}$ is dominant, and subsequent
conversion/oscillation of $A^{\prime}$ into $A$ can occur, it is severely
suppressed at small mass of $A^{\prime}$. Hence we look at annihilation into
SM in this work.

Visible annihilation channels via an $s$-channel virtual photon, $Q\bar{Q}\to
A^{*}\to{\rm SM}$, will occur at a rate suppressed by $\epsilon$ but will
result in immediate release of energy in the form of SM charged particles.

The annihilation cross-section of $Q\bar{Q}$ into $e^{+}e^{-}$ is given by,

$\displaystyle\sigma_{\rm ann}^{\rm SM}v$
$\displaystyle=\frac{\pi\epsilon^{2}\alpha^{2}}{3m_{Q}^{2}}\frac{\pi\alpha_{D}}{v_{Q}\left(1-e^{-\pi\alpha_{D}/v_{Q}}\right)}$
(26)

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,

$\displaystyle(n_{\rm loc}^{\rm ann})^{2}\sigma_{\rm ann}^{\rm SM}v_{\rm
th}V=400\left(\frac{\epsilon}{10^{-5}}\right)^{2}\left(\frac{n_{Q}}{10^{8}/\textrm{cm}^{3}}\right)^{2}$
(27) $\displaystyle\times\frac{V}{{22000\,\rm m}^{3}}\left(\frac{\rm
GeV}{m_{Q}}\right)^{2}\frac{1}{\textrm{year}}\times\\{4\times
10^{5}\alpha_{D}\sqrt{m_{Q}/\textrm{GeV}}\\}.$

which potentially may result in a very strong sensitivity to $\epsilon\times
n_{Q}m_{Q}^{-1}$. The quantity between the $\\{\\}$ provides the additional
Sommerfeld enhancement. Sensitivity to $\alpha_{D}$ enters through this
enhancement factor as well as through $n_{\rm loc}^{\rm ann}$ in Eqn. 25.

Fig. 5 illustrates this through limit contours for the different input values
of $f_{Q}$. We take $\mathcal{O}(200)$ events per year to be roughly a
limiting count rate. The left panel corresponds to $\kappa=0.1$ (small
$\alpha_{D}$) and the right panel corresponds to $\kappa=10^{-4}$ (large
$\alpha_{D}$). For a fixed small $\epsilon$, larger $\alpha_{D}$ results in
larger Sommerfeld enhancement and thus stricter limits on $f_{Q}$. However,
sensitivity to large $\epsilon$ is lost for large $\alpha_{D}$ because this
results in high annihilation rates as well, and hence extremely small local
number densities $n_{\rm loc}^{\rm ann}$. Nonetheless new parameter space is
probed for both the small $\alpha_{D}$ and large $\alpha_{D}$ regimes.

##  VI Electrostatic accelerators

Figure 6: Scheme of electrostatic accelerator concept. Positive mCPs between
the two plates held at a potential difference $\Delta V$ get accelerated
towards the more negatively charged plate and enter the low threshold detector
to be detected.

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}\sim 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 ($\Delta V$) in the megavolt range are used in nuclear physics
experiments. Examples among these are LUNA ($\Delta V=3.5$ MV) sen2019high ,
JUNA ($\Delta V=0.4$ MV) liu2016progress  and CASPAR ($\Delta V=1.1$ MV)
Robertson:2016llv .

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^{\prime}$ in the presence of a
number density $n_{Q}$ of mCPs $Q$ is given by,

$\Pi=g_{Q}^{2}\frac{n_{Q}}{m_{Q}}\approx(3\times
10^{-5}\textrm{eV})^{2}\frac{n_{Q}}{10^{14}/\textrm{cm}^{3}}\frac{\textrm{GeV}}{m_{Q}}$
(28)

Here $g_{Q}\rightarrow 1$ is set, to be conservative. The range of the
electric field is thus,

$\lambda\approx
7\textrm{mm}\times\sqrt{\frac{10^{14}/\textrm{cm}^{3}}{n_{Q}}}\sqrt{\frac{m_{Q}}{\textrm{GeV}}}$
(29)

In other words, the screening length is larger than 1 meter for
$n_{Q}\lesssim\frac{10^{8}}{\textrm{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_{\rm thr}$ required for detection, it is clear
that in order to have sensitivity one should require

$\epsilon e\Delta V>E_{\rm thr},$ (30)

where $\Delta 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 sen2019high
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

$\displaystyle\Phi[E>E_{\rm thr}]$ $\displaystyle=2\pi rL\left(1-\frac{E_{\rm
thr}}{\epsilon e\Delta V}\right)\frac{f_{Q}\eta\rho_{Q}}{m_{Q}}v_{\rm th}$
$\displaystyle\times\textrm{Min}[1,\frac{r\epsilon e\Delta V}{L\sqrt{TE_{\rm
thr}}}]$ $\displaystyle=2\times 10^{21}\textrm{Hz}f_{Q}\left(\frac{\rm
GeV}{m_{Q}}\right)^{\frac{5}{2}}\left(1-\frac{E_{\rm thr}}{\epsilon e\Delta
V}\right)$ $\displaystyle\times\textrm{Min}[1,\frac{r\epsilon e\Delta
V}{L\sqrt{TE_{\rm thr}}}]$ (31)

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.

Figure 7: The rates of accelerated mCDM for $f_{Q}=10^{-8}$ in an
electrostatic accelerator at $1$ km depth. In the Left panel, the rate
corresponding to $m_{Q}\gtrsim 1$ GeV, where top-down accumulation is dominant
is shown with the requirement for the final energy to exceed 10 eV, 100 eV, 1
keV. In the Right panel, the rate corresponding to $m_{Q}\lesssim 1$ GeV,
where bottom-up accumulation is dominant is shown with the requirement for the
final energy to exceed 100 meV, 1 eV, 10 eV.

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

$P\sim\sigma_{Q,\rm atom}\times n_{\rm atoms}\ell,~{}~{}{\rm if\,}P<1.$ (32)

where $\ell$ is taken to be $\sim 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_{\rm atoms}\simeq 4.5\times 10^{22}$ cm-3 to
correspond to the number density of germanium atoms in Ge crystal. Finally,
$\sigma_{Q,\rm atom}$ is the cross section leading to atomic recoil or
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 $\sim 0.1$. The transfer of energy in
elastic scattering on an atom is the most efficient if $m_{Q}\sim m_{\rm
atom}$. Using a non-pertubative estimate (Appendix A) for such a cross
section, $\sigma\propto 4\pi\times(\mu v_{Q})^{-2}\propto 10^{-23}{\rm
cm}^{2}$ for a typical mass of mCP in the 100 GeV range, and its kinetic
energy in the $\sim$ 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 CaWO4 of CRESST Abdelhameed:2019hmk .

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 $\sim a_{0}$, its interaction with an incoming mCP to scale as
$U\sim\alpha\epsilon a_{0}^{-1}$, perturbation theory is valid as long as
$U\times\Delta t\sim\alpha\epsilon v_{Q}^{-1}<1$ Landau1981Quantum . Taking
$m_{Q}$ below a few GeV, $\epsilon$ in the $10^{-6}-10^{-3}$ range and
$v_{Q}\sim(|\epsilon e\Delta
V|/m_{Q})^{1/2}\sim\epsilon^{1/2}\times(10^{-2}-10^{-1})$, we see that the
perturbativity condition is satisfied, and the cross section can be estimated
to scale as $\pi a_{0}^{2}\times(\alpha\epsilon v_{Q}^{-1})^{2}$. Again, this
estimate exceeds $1/(n_{\rm atom}\ell)$ making $P\sim 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.

Figure 8: The rates of accelerated mCPs produced in meson decays in the
atmosphere in an electrostatic accelerator are shown with the requirement for
the final energy to exceed 10 eV, 100 eV, 1 keV

Fig. 7 summarizes possible counting rates of accelerated mCPs for
$f_{Q}=10^{-8}$. Since the probability of detection can indeed be $P\sim
O(1)$, these rates translate to possible counting rates inside a detector
placed along the path of the accelerated particles. Given that it is realistic
to find detectors with background counts as low as $10^{-3}$ Hz, it is clear
that a dedicated search along the lines suggested in this section could probe
$f_{Q}$ down to extremely small values. For $m_{Q}\gtrsim 1$ GeV, Fig. 7
(Left), achieving sensitivity to $f_{Q}$ as small as $10^{-20}$ looks
realistic, with $f_{Q}\approx 10^{-8}$ achievable all the way to very large
$m_{Q}$. (At large $m_{Q}\gg 100$ GeV one should be cognizant of the fact that
the recoil energy of atoms drops as $m_{\rm atom}/m_{Q}$, and ionization of
atoms may be suppressed if $\alpha\epsilon v_{Q}^{-1}$ 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. Knapen:2017xzo ), even $\epsilon$’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 $\\{\epsilon,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}$cm3.

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},\epsilon\\}$ 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},\epsilon\\}$ 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 Aprile:2020tmw  is consistent with a monochromatic
signal (see e.g. Alonso-Alvarez:2020cdv ; An:2020bxd ; Baryakhtar:2020rwy ;
An:2020tcg ). 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 $\sim\epsilon\Delta V$. For MV-type
electrostatic accelerators, and for a generous range of $\epsilon$,
$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.

##  VIII Acknowledgements

We are indebted to Drs. A. Berlin and H. Liu for valuable critical comments.
We are also grateful to Drs. R. Harnik and R. Plestid for earlier
collaboration on the mCP project. M.P. would like to thank Dr. T. Bringmann
for earlier discussions of related ideas. M.P. is supported in part by U.S.
Department of Energy (Grant No. desc0011842).

##  Appendix A mCP-atom scatterting using Hulthen potential

For a mCP of mass $M$, and charge $\epsilon$ (with electron charge 1), the
perturbative expression for an atomic scattering can be written as

$\frac{d\sigma}{dq^{2}}=4\pi
Z^{2}\frac{\alpha^{2}\epsilon^{2}}{v^{2}q^{4}}F^{2}(q),$ (33)

where $F(q)$ is atomic form factor normalized to 1 at high-momentum transfer
(which corresponds to elastic scattering on un-screened nucleus). The same
cross section can also be written as,

$\frac{d\sigma}{d\Omega}=4Z^{2}\frac{\alpha^{2}\epsilon^{2}}{q^{4}}F^{2}(q)$
(34)

We will take the simplest ansatz for $F(q)$ that nevertheless captures the
main physics regimes:

$F(q)=\frac{a_{z}^{2}q^{2}}{1+a_{z}^{2}q^{2}}$ (35)

where

$a_{z}=\frac{a_{0}}{4}\left(\frac{9\pi^{2}}{2Z}\right)^{\frac{1}{3}}\sim
0.89\frac{a_{0}}{Z^{\frac{1}{3}}}\sim\frac{1}{4.2Z^{\frac{1}{3}}\textrm{keV}}.$
(36)

The momentum transfer cross-section is relevant to calculate the overburden
required for thermalization as well as the eventual terminal velocity. It is
given by

$\sigma_{T}=\frac{\pi}{2(\mu v)^{4}}\int_{0}^{2\mu
v}\frac{d\sigma}{d\Omega}q^{2}dq^{2}.$ (37)

Setting $a_{z}\mu v=R$, we get,

$\displaystyle\sigma_{T}$
$\displaystyle=2\pi\alpha^{2}Z^{2}\epsilon^{2}\frac{\log\left(4R^{2}+1\right)-\frac{4R^{2}}{4R^{2}+1}}{\mu^{2}v^{4}}$
$\displaystyle=2\pi a_{z}^{2}{\cal
B}^{2}\left(\log\left(4R^{2}+1\right)-\frac{4R^{2}}{4R^{2}+1}\right)$ (38)

Here ${\cal B}=\frac{Z\alpha\epsilon}{a_{z}\mu v^{2}}$. It is easy to see that
at $R\gg 1$ the cross section scales as $v^{-4}$, while in the opposite
regime, $R\ll 1$, it does not depend on velocity.

The perturbative answer is valid only if either $Z\epsilon\alpha a_{z}\mu\ll
1$, which for $Z_{\rm rock}\sim 20$ translates to,

$\epsilon\ll 7.5\times 10^{-5}\frac{\rm GeV}{\mu}$ (39)

or if $\frac{Z\epsilon\alpha}{v}\ll 1$, for $Z_{\rm rock}\sim 20$,

$\epsilon\ll 6.3v$ (40)

For larger coupling, expressions are available in the classical limit, $\mu
a_{z}v\gg 1$ which is not valid here, $0.6\sqrt{\frac{\mu}{\rm GeV}}\sim 1$.

We will instead use the expressions derived from Hulthen potential in Ref.
Tulin:2013teo .

$\sigma_{T}=\frac{4\pi}{\mu^{2}v^{2}}\sin^{2}\delta_{0}$ (41)

Here,

$\delta_{0}=\textrm{arg}\left(\frac{i\Gamma\left(\frac{i2\mu
a_{z}v}{\kappa}\right)}{\Gamma(\lambda_{+})\Gamma(\lambda_{-})}\right)$ (42)
$\displaystyle\lambda_{\pm}$ $\displaystyle=1+\frac{i\mu
a_{z}v}{\kappa}\pm\sqrt{\frac{Z\epsilon\alpha 2\mu}{\kappa
m_{\phi}}-\frac{\mu^{2}a_{z}^{2}v^{2}}{\kappa^{2}}}\quad\textrm{attractive}$
$\displaystyle=1+\frac{i\mu a_{z}v}{\kappa}\pm i\sqrt{\frac{Z\epsilon\alpha
2\mu}{\kappa
m_{\phi}}+\frac{\mu^{2}a_{z}^{2}v^{2}}{\kappa^{2}}}\quad\textrm{repulsive}$
(43)

We find that a very good approximation for both attractive and repulsive
interactions,

$\displaystyle\langle\sigma_{T}\rangle_{\rm th}$
$\displaystyle\sim\textrm{Min}\left(\frac{16\pi
Z^{2}\alpha^{2}\epsilon^{2}}{\mu^{2}v_{\rm th}^{4}},\frac{4\pi}{\mu_{\rm
rock,Q}^{2}v_{\rm th}^{2}}\right)$
$\displaystyle\langle\sigma_{T}v\rangle_{\rm th}$
$\displaystyle\sim\textrm{Min}\left(\frac{8\pi
Z^{2}\alpha^{2}\epsilon^{2}}{\mu^{2}v_{\rm th}^{4}}\times v_{\rm
th},\frac{2.2\pi}{\mu_{\rm rock,Q}^{2}v_{\rm th}^{2}}\times v_{\rm th}\right)$
$\displaystyle\langle\sigma_{T}v^{3}\rangle_{\rm th}$
$\displaystyle\sim\textrm{Min}\left(\frac{5\pi
Z^{2}\alpha^{2}\epsilon^{2}}{\mu^{2}v_{\rm th}^{4}}\times v_{\rm
th}^{3},\frac{2.2\pi}{\mu_{\rm rock,Q}^{2}v_{\rm th}^{2}}\times v_{\rm
th}^{3}\right)$ (44)

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

$v_{\rm th}=\sqrt{\frac{2\textrm{T}}{m_{Q}}}\sim 7\times
10^{-6}\sqrt{\frac{\rm GeV}{m_{Q}}}$ (45)

The differential angular flux coming in per infinitesimal pipe length is

$\frac{d\Phi}{dld\cos{\theta}d\varphi}=4r\frac{f_{Q}\eta\rho_{Q}}{m_{Q}}v_{\rm
th}$ (46)

Here $\theta\in\\{0,\frac{\pi}{2}\\}$, is the angle between the incoming
particle velocity and the beam axis, while $\varphi\in\\{0,\frac{\pi}{2}\\}$
subtended by the velocity on the radial direction.

The time spent by the particles inside the pipe is,

$\tau(\theta,\varphi)=\frac{2r}{v_{\rm th}\sin{\theta}\cos{\varphi}}$ (47)

The maximum time the particles can spend (because of acceleration along the
beam axis) is,

$\tau_{\rm max}=\sqrt{\frac{2l}{a}}=\sqrt{2lL}\sqrt{\frac{m_{Q}}{\epsilon
e\Delta V}}$ (48)

If a particle spends time $\tau=\textrm{Min}[\tau_{\rm
max},\tau(\theta,\varphi)]$, inside, it is accelerated to,

$E_{Q}=\frac{1}{2}m_{Q}v_{f}^{2}=\frac{1}{2}m_{Q}a^{2}\tau^{2}$ (49)

now,

$E_{Q}[{\rm max}]=\frac{1}{2}m_{Q}a^{2}\tau_{\rm max}^{2}=\frac{l}{L}\epsilon
e\Delta V$ (50)

Given a threshold for subsequent detection $E_{\rm thr}$, there is an $l_{\rm
min}$,

$l_{\rm min}=\frac{E_{\rm thr}}{\epsilon e\Delta V}L$ (51)

Furthermore, we require that the particles enter at the correct angle, this is
satisfied if,

$\sin{\theta}\cos{\varphi}<\frac{r\epsilon e\Delta V}{L\sqrt{TE_{\rm thr}}}$
(52)

Putting this together we get,

$\displaystyle\Phi[E>E_{\rm thr}]=$ $\displaystyle\int_{l_{\rm min}}^{L}dl\int
d\cos{\theta}d\varphi\frac{d\Phi}{dld\cos{\theta}d\varphi}$ (53)
$\displaystyle\times\Theta\left(\frac{r\epsilon e\Delta V}{L\sqrt{TE_{\rm
thr}}}-sin{\theta}\cos{\varphi}\right)$ (54)

This can be simplified further to

$\displaystyle\Phi[E>E_{\rm thr}]=$ $\displaystyle 2\pi rL\left(1-\frac{E_{\rm
thr}}{\epsilon e\Delta V}\right)\frac{f_{Q}\eta\rho_{Q}}{m_{Q}}$ (55)
$\displaystyle\times v_{\rm th}\textrm{Min}\left[1,\frac{r\epsilon e\Delta
V}{L\sqrt{TE_{\rm thr}}}\right].$ (56)

##  Appendix C Other probes of terrestrial milli-charged dark matter

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 Neufeld:2018slx ;
Neufeld:2019xes ; Rajendran:2020tmw . 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

$\lambda_{\textrm{MFP}}=\frac{1}{n_{\rm shield}\sigma_{T}}\approx
10^{-6}\textrm{cm}\frac{\mu}{\rm GeV}$ (57)

Here we approximate $n_{\rm shield}=\frac{10^{22}}{{\rm cm}^{3}}$ and take the
cross-section in Eqn. 5 to be $\sigma_{T}\approx\frac{4\pi}{\mu^{2}v_{\rm
th}^{2}}$. This MFP is much smaller than any realistic shielding. As a result
mCDM will not be constrained by cryogens.

LHC beam lifetime

In Neufeld:2018slx , 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 $\epsilon^{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}\sim n_{\rm atoms}\epsilon^{-2}\sim
10^{10}\textrm{cm}^{-3}\times(10^{-2}/\epsilon)^{2},$ (58)

where we took a realistic residual gas density at $10^{6}$ particles per cm3.
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: $\mathbf{180^{\textrm{m}}}$ Ta.

In Lehnert:2019tuw ; Pospelov:2019vuf  the non-observation of the decay of
isomeric Tantalum was used to set limit on terrestrial DM. The form-factor to
scatter with $180^{\textrm{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.

$P_{g}(T)=e^{-\frac{E_{g}}{T}}\approx e^{-11\frac{\epsilon^{2}}{v_{\rm
th}^{2}}}$ (59)

where $E_{g}=2\mu\left(\pi\alpha Z_{\rm Ta}\epsilon\right)^{2}$ and $T$ is the
ambient temperature. This factor evaluates to tiny values rendering the cross-
section too small for Lehnert:2019tuw  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 $180$m $\to 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 $\epsilon$ is very small (and the mass is in $\sim 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 Earth-bound exotic species Neufeld:2018slx . 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.

## References

  * [1] Andrew Haas, Christopher S. Hill, Eder Izaguirre, and Itay Yavin.  Looking for milli-charged particles with a new experiment at the LHC.  Phys. Lett. B, 746:117–120, 2015.
  * [2] Eder Izaguirre and Itay Yavin.  New window to millicharged particles at the LHC.  Phys. Rev. D, 92(3):035014, 2015.
  * [3] Gabriel Magill, Ryan Plestid, Maxim Pospelov, and Yu-Dai Tsai.  Millicharged particles in neutrino experiments.  Phys. Rev. Lett., 122(7):071801, 2019.
  * [4] Kevin J. Kelly and Yu-Dai Tsai.  Proton fixed-target scintillation experiment to search for millicharged dark matter.  Phys. Rev. D, 100(1):015043, 2019.
  * [5] Roni Harnik, Zhen Liu, and Ornella Palamara.  Millicharged Particles in Liquid Argon Neutrino Experiments.  JHEP, 07:170, 2019.
  * [6] R. Acciarri et al.  Improved Limits on Millicharged Particles Using the ArgoNeuT Experiment at Fermilab.  Phys. Rev. Lett., 124(13):131801, 2020.
  * [7] Ryan Plestid, Volodymyr Takhistov, Yu-Dai Tsai, Torsten Bringmann, Alexander Kusenko, and Maxim Pospelov.  New Constraints on Millicharged Particles from Cosmic-ray Production.  2 2020.
  * [8] Itay M. Bloch, Andrea Caputo, Rouven Essig, Diego Redigolo, Mukul Sholapurkar, and Tomer Volansky.  Exploring New Physics with O(keV) Electron Recoils in Direct Detection Experiments.  6 2020.
  * [9] Saeid Foroughi-Abari, Felix Kling, and Yu-Dai Tsai.  FORMOSA: Looking Forward to Millicharged Dark Sectors.  10 2020.
  * [10] Roni Harnik, Ryan Plestid, Maxim Pospelov, and Harikrishnan Ramani.  Millicharged Cosmic Rays and Low Recoil Detectors.  10 2020.
  * [11] Giacomo Marocco and Subir Sarkar.  Blast from the past: Constraints on the dark sector from the BEBC WA66 beam dump experiment.  SciPost Phys., 10:043, 2021.
  * [12] Bob Holdom.  Two U(1)’s and Epsilon Charge Shifts.  Phys. Lett. B, 166:196–198, 1986.
  * [13] Rennan Barkana, Nadav Joseph Outmezguine, Diego Redigolo, and Tomer Volansky.  Strong constraints on light dark matter interpretation of the EDGES signal.  Phys. Rev. D, 98(10):103005, 2018.
  * [14] Kimberly K. Boddy, Vera Gluscevic, Vivian Poulin, Ely D. Kovetz, Marc Kamionkowski, and Rennan Barkana.  Critical assessment of CMB limits on dark matter-baryon scattering: New treatment of the relative bulk velocity.  Phys. Rev. D, 98(12):123506, 2018.
  * [15] Hongwan Liu, Nadav Joseph Outmezguine, Diego Redigolo, and Tomer Volansky.  Reviving Millicharged Dark Matter for 21-cm Cosmology.  Phys. Rev. D, 100(12):123011, 2019.
  * [16] E. Aprile et al.  Excess electronic recoil events in XENON1T.  Phys. Rev. D, 102(7):072004, 2020.
  * [17] Yasaman Farzan and Meshkat Rajaee.  Pico-charged particles explaining 511 keV line and XENON1T signal.  7 2020.
  * [18] Gaia Lanfranchi, Maxim Pospelov, and Philip Schuster.  The Search for Feebly-Interacting Particles.  11 2020.
  * [19] Ely D. Kovetz, Vivian Poulin, Vera Gluscevic, Kimberly K. Boddy, Rennan Barkana, and Marc Kamionkowski.  Tighter limits on dark matter explanations of the anomalous EDGES 21 cm signal.  Phys. Rev. D, 98(10):103529, 2018.
  * [20] Cyril Creque-Sarbinowski, Lingyuan Ji, Ely D. Kovetz, and Marc Kamionkowski.  Direct millicharged dark matter cannot explain the EDGES signal.  Phys. Rev. D, 100(2):023528, 2019.
  * [21] Tracy R. Slatyer and Chih-Liang Wu.  Early-Universe constraints on dark matter-baryon scattering and their implications for a global 21 cm signal.  Phys. Rev. D, 98(2):023013, 2018.
  * [22] Albert Stebbins and Gordan Krnjaic.  New Limits on Charged Dark Matter from Large-Scale Coherent Magnetic Fields.  JCAP, 12:003, 2019.
  * [23] M. Kachelriess and J. Tjemsland.  Reacceleration of charged dark matter.  JCAP, 10:001, 2020.
  * [24] Sacha Davidson, Steen Hannestad, and Georg Raffelt.  Updated bounds on millicharged particles.  JHEP, 05:003, 2000.
  * [25] Jae Hyeok Chang, Rouven Essig, and Samuel D. McDermott.  Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle.  JHEP, 09:051, 2018.
  * [26] Hendrik Vogel and Javier Redondo.  Dark Radiation constraints on minicharged particles in models with a hidden photon.  JCAP, 02:029, 2014.
  * [27] Maxim Pospelov, Surjeet Rajendran, and Harikrishnan Ramani.  Metastable Nuclear Isomers as Dark Matter Accelerators.  Phys. Rev. D, 101(5):055001, 2020.
  * [28] Björn Lehnert, Harikrishnan Ramani, Mikael Hult, Guillaume Lutter, Maxim Pospelov, Surjeet Rajendran, and Kai Zuber.  Search for dark matter induced de-excitation of 180Tam.  Phys. Rev. Lett., 124(18):181802, 2020.
  * [29] Surjeet Rajendran and Harikrishnan Ramani.  A Composite Solution to the Neutron Bottle Anomaly.  8 2020.
  * [30] David A. Neufeld, Glennys R. Farrar, and Christopher F. McKee.  Dark Matter that Interacts with Baryons: Density Distribution within the Earth and New Constraints on the Interaction Cross-section.  Astrophys. J., 866(2):111, 2018.
  * [31] Giorgio Busoni, Andrea De Simone, Pat Scott, and Aaron C. Vincent.  Evaporation and scattering of momentum- and velocity-dependent dark matter in the Sun.  JCAP, 10:037, 2017.
  * [32] Valerio De Luca, Andrea Mitridate, Michele Redi, Juri Smirnov, and Alessandro Strumia.  Colored Dark Matter.  Phys. Rev. D, 97(11):115024, 2018.
  * [33] Rebecca K. Leane and Juri Smirnov.  Exoplanets as New Sub-GeV Dark Matter Detectors.  9 2020.
  * [34] Maxim Pospelov and Adam Ritz.  Resonant scattering and recombination of pseudo-degenerate WIMPs.  Phys. Rev. D, 78:055003, 2008.
  * [35] Haipeng An, Maxim Pospelov, and Josef Pradler.  Direct constraints on charged excitations of dark matter.  Phys. Rev. Lett., 109:251302, 2012.
  * [36] Bartosz Fornal, Benjamin Grinstein, and Yue Zhao.  Dark Matter Capture by Atomic Nuclei.  Phys. Lett. B, 811:135869, 2020.
  * [37] C. Broggini, D. Bemmerer, A. Caciolli, and D. Trezzi.  LUNA: Status and Prospects.  Prog. Part. Nucl. Phys., 98:55–84, 2018.
  * [38] Asher Berlin, Raffaele Tito D’Agnolo, Sebastian A.R. Ellis, Philip Schuster, and Natalia Toro.  Directly Deflecting Particle Dark Matter.  Phys. Rev. Lett., 124(1):011801, 2020.
  * [39] LD Landau, EM Lifshitz, and LP Pitaevskij.  Course of theoretical physics. vol. 10: Physical kinetics.  Oxford, 1981.
  * [40] Timon Emken, Rouven Essig, Chris Kouvaris, and Mukul Sholapurkar.  Direct Detection of Strongly Interacting Sub-GeV Dark Matter via Electron Recoils.  JCAP, 09:070, 2019.
  * [41] Rouven Essig, Marivi Fernandez-Serra, Jeremy Mardon, Adrian Soto, Tomer Volansky, and Tien-Tien Yu.  Direct Detection of sub-GeV Dark Matter with Semiconductor Targets.  JHEP, 05:046, 2016.
  * [42] Nima Arkani-Hamed, Douglas P. Finkbeiner, Tracy R. Slatyer, and Neal Weiner.  A Theory of Dark Matter.  Phys. Rev. D, 79:015014, 2009.
  * [43] Maxim Pospelov and Adam Ritz.  Astrophysical Signatures of Secluded Dark Matter.  Phys. Lett. B, 671:391–397, 2009.
  * [44] Haipeng An, Mark B. Wise, and Yue Zhang.  Effects of Bound States on Dark Matter Annihilation.  Phys. Rev. D, 93(11):115020, 2016.
  * [45] Marco Cirelli, Paolo Panci, Kalliopi Petraki, Filippo Sala, and Marco Taoso.  Dark Matter’s secret liaisons: phenomenology of a dark U(1) sector with bound states.  JCAP, 05:036, 2017.
  * [46] P.F. Smith, J.R.J. Bennett, G.J. Homer, J.D. Lewin, H.E. Walford, and W.A. Smith.  A SEARCH FOR ANOMALOUS HYDROGEN IN ENRICHED D-2 O, USING A TIME-OF-FLIGHT SPECTROMETER.  Nucl. Phys. B, 206:333–348, 1982.
  * [47] T.K. Hemmick et al.  A Search for Anomalously Heavy Isotopes of Low $Z$ Nuclei.  Phys. Rev. D, 41:2074–2080, 1990.
  * [48] E. Barbosa de Souza et al.  First search for a dark matter annual modulation signal with NaI(Tl) in the Southern Hemisphere by DM-Ice17.  Phys. Rev. D, 95(3):032006, 2017.
  * [49] A Sen, G Domínguez-Cañizares, NC Podaru, DJW Mous, M Junker, G Imbriani, and V Rigato.  A high intensity, high stability 3.5 mv singletron™ accelerator.  Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 450:390–395, 2019.
  * [50] WeiPing Liu, ZhiHong Li, JiangJun He, XiaoDong Tang, Gang Lian, Zhu An, JianJun Chang, Han Chen, QingHao Chen, XiongJun Chen, et al.  Progress of jinping underground laboratory for nuclear astrophysics (juna).  SCIENCE CHINA Physics, Mechanics & Astronomy, 59(4):642001, 2016\.
  * [51] Daniel Robertson, Manoel Couder, Uwe Greife, Frank Strieder, and Michael Wiescher.  Underground nuclear astrophysics studies with CASPAR.  EPJ Web Conf., 109:09002, 2016.
  * [52] A.H. Abdelhameed et al.  First results from the CRESST-III low-mass dark matter program.  Phys. Rev. D, 100(10):102002, 2019.
  * [53] L. D. Landau and L. M. Lifshitz.  Quantum Mechanics Non-Relativistic Theory, Third Edition: Volume 3.  Butterworth-Heinemann, 3 edition, January 1981.
  * [54] Simon Knapen, Tongyan Lin, and Kathryn M. Zurek.  Light Dark Matter: Models and Constraints.  Phys. Rev. D, 96(11):115021, 2017.
  * [55] Gonzalo Alonso-Álvarez, Fatih Ertas, Joerg Jaeckel, Felix Kahlhoefer, and Lennert J. Thormaehlen.  Hidden Photon Dark Matter in the Light of XENON1T and Stellar Cooling.  JCAP, 11:029, 2020.
  * [56] Haipeng An, Maxim Pospelov, Josef Pradler, and Adam Ritz.  New limits on dark photons from solar emission and keV scale dark matter.  6 2020.
  * [57] Masha Baryakhtar, Asher Berlin, Hongwan Liu, and Neal Weiner.  Electromagnetic Signals of Inelastic Dark Matter Scattering.  6 2020.
  * [58] Haipeng An and Daneng Yang.  Direct detection of freeze-in inelastic dark matter.  6 2020.
  * [59] Sean Tulin, Hai-Bo Yu, and Kathryn M. Zurek.  Beyond Collisionless Dark Matter: Particle Physics Dynamics for Dark Matter Halo Structure.  Phys. Rev. D, 87(11):115007, 2013.
  * [60] David A. Neufeld and Daniel J. Brach-Neufeld.  Dark Matter that Interacts with Baryons: Experimental Limits on the Interaction Cross-section for 27 Atomic Nuclei, and Resultant Constraints on the Particle Properties.  Astrophys. J., 877(1):8, 2019.