# Low-Energy Signals from the Formation of Dark Matter-Nuclear Bound States

Asher Berlin<EMAIL_ADDRESS>Center for Cosmology and Particle Physics,
Department of Physics, New York University, New York, NY 10003, USA  Hongwan
Liu<EMAIL_ADDRESS>Center for Cosmology and Particle Physics,
Department of Physics, New York University, New York, NY 10003, USA
Department of Physics, Princeton University, Princeton, New Jersey, 08544, USA
Maxim Pospelov<EMAIL_ADDRESS>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 matter particles may bind with nuclei if there exists an attractive force
of sufficient strength. We show that a dark photon mediator of mass
$\sim(10-100)\ \text{MeV}$ that kinetically mixes with Standard Model
electromagnetism at the level of $\sim 10^{-3}$ generates keV-scale binding
energies between dark matter and heavy elements, while forbidding the ability
to bind with light elements. In underground direct detection experiments, the
formation of such bound states liberates keV-scale energy in the form of
electrons and photons, giving rise to mono-energetic electronic signals with a
time-structure that may contain daily and seasonal modulations. We show that
data from liquid-xenon detectors provides exquisite sensitivity to this
scenario, constraining the galactic abundance of such dark particles to be at
most $\sim 10^{-18}-10^{-12}$ of the galactic dark matter density for masses
spanning $\sim(1-10^{5})\ \text{GeV}$. However, an exponentially small
fractional abundance of these dark particles is enough to explain the observed
electron recoil excess at XENON1T.

_Introduction._ Over the years, dark matter (DM) direct detection experiments
have become extraordinarily sensitive to sub-MeV energy deposition by exotic
sources, with thresholds recently extending down to sub-keV energies. While
the primary motivation for these experiments is to search for the elastic
scattering of weakly interacting massive particles (WIMPs) off nuclei, the
scope of these searches now includes, e.g., electron scattering, DM
absorption, and exo- and endo-thermic inelasticity in DM-nucleus scattering
(see, e.g., Ref. Lin (2019) for a review). Each of these searches seeks to
measure the energy deposited from either the kinetic (scattering) or mass
energy (absorption, inelasticity) of the incoming DM particle. In this paper,
we point out a third, distinct alternative: the energy released due to the
formation of a bound state between DM and a Standard Model (SM) nucleus.

The phenomenology of DM-nuclear bound state formation has been studied
previously in the literature. For instance,
$\text{\,}\mathrm{Me\kern-0.50003ptV}$-scale DM-nuclear binding was
investigated in Refs. Pospelov and Ritz (2007); An  _et al._ (2012); Fornal
_et al._ (2020), while Refs. Wallemacq (2013, 2014); Wallemacq and Cudell
(2015); Laletin and Cudell (2019) focused on
$\text{\,}\mathrm{ke\kern-0.50003ptV}$-scale bound states involving dark atoms
to explain the DAMA anomaly Bernabei  _et al._ (2008, 2010, 2018). In contrast
to these previously studied models, all of which require some degree of
intricate model-building, we instead consider bound states arising from one of
the simplest and most studied models in the literature within the last fifteen
years Boehm and Fayet (2004); Arkani-Hamed  _et al._ (2009); Pospelov  _et
al._ (2008): a DM particle $\chi$ charged under a massive kinetically-mixed
dark photon. We focus on a scenario where $\chi$ possesses a sizeable
interaction with normal matter and constitutes a small fraction $f_{\chi}\ll
1$ of the total galactic DM density. This model naturally leads to _i)_
$\text{\,}\mathrm{ke\kern-0.50003ptV}$-scale DM-nuclear binding energies
$E_{B}$ and _ii)_ preferential binding to heavy nuclei, such that upon
penetrating the terrestrial overburden, $\chi$ only binds with the much
heavier nuclei (such as xenon and thallium) commonly found in underground DM
detectors.

In a direct detection experiment whose target material consists of atoms $A$
of sufficiently large atomic number and mass, the process of DM-atom
“recombination” $\chi+A\to(\chi A)+E_{B}$ releases electromagnetic energy
$E_{B}$ equal to the binding energy of $(\chi A)$.111We adopt the notation
where a bound state is denoted by parentheses surrounding the names of the
constituent particles. For this reaction to occur, the minimum required
coupling between $\chi$ and $A$ is sufficiently large such that the galactic
$\chi$ population quickly thermalizes upon encountering Earth’s environment,
cooling down to terrestrial temperatures. The implications of thermalizing
with the terrestrial environment are two-fold. First, the terrestrial $\chi$
density is drastically enhanced compared to the galactic population, due to
conservation of flux (the “traffic-jam” scenario discussed in Refs. Pospelov
_et al._ (2020); Pospelov and Ramani (2020); Neufeld  _et al._ (2018)).
Second, there are no observable elastic scattering signals of $\chi$ despite
its large couplings and enhanced terrestrial density, since the thermal energy
of underground laboratory environments
($\lesssim$300\text{\,}\mathrm{K}$\sim$25\text{\,}\mathrm{me\kern-0.50003ptV}$$)
is well below existing kinematic thresholds. However, since $E_{B}$ naturally
lies near the keV-scale, the formation of DM-nuclear bound states is readily
detectable at large scale experiments designed to search for WIMP-nuclear
scattering.

Among such direct detection experiments, the suite of large scale dual-phase
xenon detectors plays an especially important role. For instance, ionization-
only data from XENON10 and XENON100 place some of the strongest constraints on
MeV-scale DM-electron scattering Essig  _et al._ (2017), and the large
exposure and low background counts (below $10^{-5}\text{ per kg-day-keV}$) of
the XENON1T experiment enable new benchmark sensitivity not only to WIMP-
nuclear scattering but also to sub-keV electronic recoils Aprile  _et al._
(2018). Intriguingly, the XENON1T collaboration recently reported an excess of
events consistent with electron recoils with an energy deposition of $(2-3)\
\text{keV}$ Aprile  _et al._ (2020). This may be consistent with a variety of
recently proposed new physics models, all of which invoke a substantial flux
of particles that feebly interact with normal matter (e.g., dark photons
Alonso-Álvarez  _et al._ (2020); An  _et al._ (2020), neutrinos and dark
radiation Bloch  _et al._ (2020); McKeen  _et al._ (2020); Farzan and Rajaee
(2020); Brdar  _et al._ (2021); Kuo  _et al._ (2021), and exothermic DM
Baryakhtar  _et al._ (2020); Bramante and Song (2020); Bloch  _et al._ (2020);
Harigaya  _et al._ (2020)). In this Letter, we find that the observed rate at
XENON1T may be explained as a result of a strongly-coupled particle that makes
up an extremely small fraction $f_{\chi}$ of the galactic DM density. In
particular, we find that a DM subcomponent that binds to xenon nuclei with
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ is a viable explanation to
this anomaly. More generally, these signals are significantly constrained by
XENON1T for fractional abundances greater than $10^{-18}\lesssim
f_{\chi}\lesssim 10^{-12}$ and particle masses spanning
$$1\text{\,}\mathrm{Ge\kern-0.50003ptV}$\lesssim
m_{\chi}\lesssim$30\text{\,}\mathrm{Te\kern-0.50003ptV}$$.

_Model and Bound State Parameter Space._ We consider a subcomponent of DM
$\chi$ that is directly charged under a new massive dark $U(1)$ gauge boson
$A^{\prime}_{\mu}$ that kinetically mixes with the SM photon,222For
concreteness, our calculations assume that $\chi$ is fermionic, although all
considerations in this paper apply equally well to scalar DM.

$\mathcal{L}\supset-\frac{\epsilon}{2}\,F_{\mu\nu}^{\prime}\,F^{\mu\nu}+\frac{m_{A^{\prime}}^{2}}{2}\,A_{\mu}^{\prime\,2}~{},$
(1)

where $m_{A^{\prime}}$ is the dark photon mass and $\epsilon\ll 1$ controls
the strength of kinetic mixing Holdom (1986). If $\epsilon$ is generated
radiatively from particles charged under both the SM and dark sector, the
natural expectation is $\epsilon\sim(\alpha_{D}\,\alpha)^{1/2}/4\pi$, where
$\alpha_{D}$ and $\alpha$ are the dark photon and SM fine-structure constant,
respectively. $\chi$ interacts with normal matter through a small effective
coupling $eq_{\text{eff}}$ where
$q_{\text{eff}}\equiv\epsilon\sqrt{\alpha_{D}/\alpha}$.

The dark photon also mediates attractive self-interactions, such that
resonances and capture to $(\chi\bar{\chi})$ bound states can significantly
reduce the cosmological $\chi$ density Arkani-Hamed  _et al._ (2009); Pospelov
and Ritz (2009); Feng  _et al._ (2010); Cirelli  _et al._ (2017); An  _et al._
(2016). It is therefore reasonable to consider a small fraction of the DM
energy density $f_{\chi}\equiv\rho_{\chi}/\rho_{{}_{\text{DM}}}$ that is
composed of such particles. In the local vicinity of the galaxy, we consider
$f_{\chi}$ to be a free parameter, noting that deviations from a standard
thermal cosmological history could result in $f_{\chi}\ll 1$.333As a concrete
example, arbitrarily small abundances of $\chi$ are cosmologically generated
provided that the reheat temperature of the universe $T_{\text{RH}}$ is
significantly smaller than $m_{\chi}$. In this case, electron annihilations
freeze-in a fractional abundance of $\chi$ corresponding to
$f_{\chi}\sim(\alpha
q_{\text{eff}})^{2}e^{-2m_{\chi}/T_{\text{RH}}}m_{\chi}\,m_{\text{pl}}/(T_{\text{RH}}\,T_{\text{mre}})$,
where $m_{\text{pl}}\sim 10^{19}\ \text{GeV}$ is the Planck mass and
$T_{\text{mre}}\sim 1\ \text{eV}$ is the temperature at matter-radiation
equality.

Figure 1: The minimum coupling $q_{\text{eff}}$ required for dark matter to
bind with various nuclei, as a function of dark matter mass $m_{\chi}$, fixing
$m_{A^{\prime}}=50\ \text{MeV}$. A dark matter-nuclear bound state exists with
an element for values of $q_{\text{eff}}$ above the corresponding line.

A massive dark photon only allows binding with heavy nuclei. Intuitively, this
selection arises because the characteristic size of the bound state
$(q_{\text{eff}}Z\alpha\mu)^{-1}$ should be smaller than the range of the
interaction $m_{A^{\prime}}^{-1}$, with both the atomic number $Z$ of the
nucleus and the DM-nuclear reduced mass $\mu$ increasing with larger nuclei.
To incorporate effects associated with the finite size of the nucleus, we use
Bargmann’s limit Bargmann (1952), which can be used to show that a DM-nuclear
bound state exists only if

$q_{\text{eff}}Z\alpha\,\mu\gtrsim
m_{A^{\prime}}\times\begin{cases}1/2&(m_{A^{\prime}}R_{\text{nuc}}\ll 1)\\\
m_{A^{\prime}}R_{\text{nuc}}/3&(m_{A^{\prime}}R_{\text{nuc}}\gg
1)~{},\end{cases}$ (2)

where the radius of a nucleus of atomic mass $A$ is
$R_{\text{nuc}}\sim$1.1\text{\,}\mathrm{fm}$\times A^{1/3}$. In Fig. 1, we
show the minimum value of $q_{\text{eff}}$ that is required for $\chi$ to bind
with various nuclei as a function of the DM mass $m_{\chi}$, fixing
$m_{A^{\prime}}=$50\text{\,}\mathrm{Me\kern-0.50003ptV}$$. We see a clear
preference to bind to heavier elements for all $\chi$ masses. This is most
pronounced for large DM masses
$m_{\chi}\gg$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$$, in which case
$\mu\simeq m_{N}$ and the minimal coupling to bind strongly depends on the
nuclear mass.

Figure 2: In the $\\{m_{A^{\prime}},q_{\text{eff}}\\}$ plane, contours of
fixed binding energy $E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ of
$(\chi\,\text{Xe}^{+})$ for different choices of $m_{\chi}$. Along the dashed
parts of the contours, $\chi$ also binds with Fe.

This leads to the intriguing possibility that binding to heavy elements, such
as xenon and thallium, is allowed, while for lighter elements, such as
nitrogen and iron, no bound state exists. In this case, $\chi$ does not bind
to the light elements it encounters when traversing Earth’s atmosphere and
crust, but does bind to heavy elements employed in direct detection
experiments. In the Bohr-like regime ($m_{A^{\prime}}\ll
q_{\text{eff}}Z\alpha\mu$), the binding energy is

$E_{B}^{(\text{Bohr})}\sim$10\text{\,}\mathrm{ke\kern-0.50003ptV}$\times\bigg{(}\frac{q_{\text{eff}}}{10^{-3}}\bigg{)}^{2}\bigg{(}\frac{Z}{54}\bigg{)}^{2}\bigg{(}\frac{\mu}{$122\text{\,}\mathrm{Ge\kern-0.50003ptV}$}\bigg{)}~{},$
(3)

thus opening up an opportunity to search for the
$\mathcal{O}($\text{\,}\mathrm{ke\kern-0.50003ptV}$)$ energy release in the
formation of such bound states in xenon-based targets. From Eq. (2), for
$m_{\chi}\gg$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$$ this occurs when the
dark photon mass is $m_{A^{\prime}}\sim q_{\text{eff}}\times\mathcal{O}(10)$
$\text{\,}\mathrm{Ge\kern-0.50003ptV}$. In Fig. 2, we show the value of
$q_{\text{eff}}$ that is required to achieve a binding energy of
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ with xenon (corresponding to
the recoil energy of the observed excess at XENON1T Aprile  _et al._ (2020))
as a function of $m_{A^{\prime}}$ for various choices of $m_{\chi}$. This is
determined by numerically solving the Schrödinger equation (see the
Supplemental Material for additional details). For $m_{\chi}\lesssim
m_{\text{Xe}}$, larger $\chi$ masses enhance $\mu$ and thus require smaller
values of $q_{\text{eff}}$ for fixed $E_{B}$, while this effect saturates for
$m_{\chi}\gg m_{\text{Xe}}$.444For
$m_{\chi}\lesssim$20\text{\,}\mathrm{Ge\kern-0.50003ptV}$$ and
$\epsilon\lesssim 10^{-3}$ (see the right panel of Fig. 3), a dark sector
coupling of $\alpha_{D}\gtrsim 1$ is necessary for our choice of model
parameters, which is reasonable if $\chi$ is a composite state. As expected,
the required value of $q_{\text{eff}}$ decreases for longer-ranged dark
photons (smaller $m_{A^{\prime}}$), saturating once $m_{A^{\prime}}$ is
smaller than the inverse Bohr radius of the DM-nuclear bound state. However,
for even smaller values of $m_{A^{\prime}}$, $\chi$ also binds with iron, the
heaviest abundant element in Earth’s crust (denoted by dashed lines in Fig.
2), thus preventing $\chi$ particles from reaching underground detectors. We
are thus motivated to consider dark photon masses of
$m_{A^{\prime}}\sim(10-100)\ \text{MeV}$.

_Terrestrial Thermalization._ Before ultimately binding to a heavy nucleus in
an underground direct detection experiment, $\chi$ thermalizes to terrestrial
temperatures after elastically scatters many times off the much lighter
elements in Earth’s atmosphere and crust. This is governed by the transfer
cross section for elastic DM-nuclear scattering $\chi N\to\chi N$, which in
the perturbative Born limit ($m_{A^{\prime}}\gg q_{\text{eff}}Z\alpha\mu$) is
approximately

$\displaystyle\sigma_{T}^{(\text{Born})}$
$\displaystyle\simeq\frac{64\pi\,Z^{2}\alpha^{2}\,q_{\text{eff}}^{2}\,\mu^{2}}{3\,m_{A^{\prime}}^{4}}$
$\displaystyle\simeq$3\text{\times}{10}^{-26}\text{\,}{\mathrm{cm}}^{2}$\times\bigg{(}\frac{q_{\text{eff}}}{10^{-3}}\bigg{)}^{2}\bigg{(}\frac{$50\text{\,}\mathrm{Me\kern-0.50003ptV}$}{m_{A^{\prime}}}\bigg{)}^{4}~{},$
(4)

where in the second line we have taken $m_{\chi}\gg m_{N}$ and set the nuclear
parameters equal to that of silicon, one of the most abundant elements in
Earth’s crust. Taking a terrestrial silicon density of
$n_{\text{Si}}\sim${10}^{22}\text{\,}{\mathrm{cm}}^{-3}$$, the typical
distance
$\ell_{\text{therm}}\sim(m_{\chi}\,m_{N}/\mu^{2})\,(n_{\text{Si}}\,\sigma_{T})^{-1}$
for $\chi$ to equilibrate to room temperature is much smaller than $1\
\text{km}$ for sub-TeV DM masses.

Earth’s gravitational field $g$ induces a radially inward bulk flow of the
thermalized $\chi$ particles. The drift velocity of this flow is
parametrically of size
$v_{\text{drift}}\sim(m_{\chi}/\mu)\,g\,\ell_{\text{mfp}}/v_{\text{rel}}$,
where $\ell_{\text{mfp}}\sim 1/(n_{N}\,\sigma_{T})$ is the mean free path for
scattering off nuclei $N$ and $v_{\text{rel}}$ is the relative thermal
velocity between $\chi$ and $N$. In the parameter space of interest, this
drift is very slow, e.g., $v_{\text{drift}}\ll 10^{-10}$ for sub-TeV DM
masses. Since the virialized galactic $\chi$ population continually bombards
Earth with a characteristic speed set by the much faster galactic wind
$v_{\text{vir}}\sim 10^{-3}$, conservation of DM flux implies that the
terrestrial energy density $\rho_{\chi}$ is greatly enhanced555The rate to
form bound $(\chi\bar{\chi})$ states on Earth is much too slow to affect the
signals discussed here. This is because although such dark interactions are
unsuppressed by $\epsilon\ll 1$, the largest terrestrial $\chi$ densities that
we consider in this work are smaller than the density of normal matter by at
least thirteen orders of magnitude. compared to the galactic density Pospelov
and Ramani (2020). Specifically,
$\rho_{\chi}\sim(v_{\text{vir}}/v_{\text{drift}})\,f_{\chi}\,\rho_{{}_{\text{DM}}}$,
where
$\rho_{{}_{\text{DM}}}\simeq$0.3\text{\,}\mathrm{Ge\kern-0.50003ptV}\text{\,}{\mathrm{cm}}^{-3}$$
is the local DM energy density.666For our choice of model parameters, if
$m_{\chi}\gg 10\ \text{TeV}$, then the thermalization distance
$\ell_{\text{therm}}$ is much larger than the detector depth $h$ (which
consists of $h\simeq 1.4\ \text{km}$ of rock at Gran Sasso). As a result,
$\chi$ particles bombarding the Earth from above do not thermalize before
reaching the detector. However, if $\ell_{\text{therm}}$ is much smaller than
the radius of the Earth, $\chi$ particles approaching the Earth from below
travel more distance through Earth’s crust and hence can thermalize above the
detector. In this case, the local value of $\rho_{\chi}$ is reduced by the
small region of solid angles corresponding to such trajectories, suppressing
the thermalized $\chi$ density at an underground detector by $\sim
h/(2\ell_{\text{therm}})$. Taking the limestone rock of the Gran Sasso
overburden to be composed of an equal mixture of calcium carbonate and
magnesium carbonate with density
$\sim$3\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$$ Miramonti (2005), we
find that the terrestrial number density of $\chi$ is approximately

$n_{\chi}\sim f_{\chi}\times 10^{8}\
\text{cm}^{-3}\times\begin{cases}\big{(}\frac{\text{GeV}}{m_{\chi}}\big{)}^{1.5}&(m_{\chi}\ll
m_{N})\\\ \big{(}\frac{\text{GeV}}{m_{\chi}}\big{)}^{2.1}&(m_{\chi}\gg
m_{N})~{}.\end{cases}$ (5)

In performing this calculation, instead of the Born-like estimate in Eq. (Low-
Energy Signals from the Formation of Dark Matter-Nuclear Bound States), we
have incorporated a full thermal average of the scattering rate utilizing the
semi-analytic results outlined in Ref. Colquhoun  _et al._ (2021). A more
detailed discussion will be provided in forthcoming work Berlin  _et al._
(2022).

_Rate of Bound State Formation._ The formation of DM-nuclear bound states
occurs by the Migdal/Auger-like ejection of an atomic electron,
$\chi+A\to(\chi A^{+})+e^{-}$, followed by subsequent relaxation of the non-
ejected atomic electrons to the ground state. Bound state formation by
bremsstrahlung emission of a final state photon, $\chi+A\to(\chi A)+\gamma$,
is subdominant, as is the case for the elastic scattering of light DM, since
the rate is suppressed by $\sim$\text{\,}\mathrm{ke\kern-0.50003ptV}$/m_{N}$
Bell  _et al._ (2020); Ibe  _et al._ (2018).

We estimate the probability of the transition $\chi+A\to(\chi A^{+})+e^{-}$
using quantum mechanical perturbation theory. In the Supplemental Material, we
derive the cross section for this process. We focus on the $s$-wave DM-nuclear
final state since it is guaranteed to exist if a bound state is allowed. For
an incoming DM particle that is also $s$-wave, the cross section for this
$s$-wave to $s$-wave transition is

$\sigma_{B}v_{\text{rel}}\simeq\frac{4\pi}{9}\,\frac{(Z\alpha\,m_{e})^{5}}{(2\mu\,E_{B})^{7/2}}\,\Big{(}\frac{\mu}{m_{N}}\Big{)}^{4}\,F_{\chi}^{2}\,F_{e}^{2}~{}.$
(6)

The numerical factors $F_{\chi}$ and $F_{e}$ depend on the initial and final
state wavefunctions of $\chi$ and the ejected electron, respectively, and must
be determined numerically. For a binding energy of
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ in xenon, we find
$F_{e}\simeq 0.7$. Additionally fixing $q_{\text{eff}}$ as in Fig. 2 and
$m_{A^{\prime}}=$50\text{\,}\mathrm{Me\kern-0.50003ptV}$$, we find that
$F_{\chi}\sim\mathcal{O}(10)$ with the precise value depending on the DM mass
(see the Supplemental Material for additional details). Upon evaluating
$F_{\chi}$ and $F_{e}$, we find that the cross section for forming DM-xenon
bound states with $E_{B}=2.5\ \text{keV}$ and
$m_{A^{\prime}}=$50\text{\,}\mathrm{Me\kern-0.50003ptV}$$ is well fit by the
functional form

$\sigma_{B}v_{\text{rel}}\simeq$6\text{\times}{10}^{-34}\text{\,}{\mathrm{cm}}^{2}$\times\left(\frac{\mu}{100\
\text{GeV}}\right)^{0.55}~{}.$ (7)

In the Supplemental Material, we show that although $p$-wave (and higher) to
$s$-wave transitions also occur, they are suppressed by the small temperature
of the thermalized DM population and are thus subdominant to the $s$-wave to
$s$-wave transition discussed above.

_Recombination Signal._ Direct detection experiments that employ heavy
elements, such as liquid noble targets, have remarkable sensitivity to the
formation of these DM-nuclear bound states. For a detector employing a target
of atomic mass $A$, the signal rate per unit target mass per unit time is

$R_{\text{sig}}=P_{\text{surv}}\,n_{\chi}\,\sigma_{B}v_{\text{rel}}\,\left(\frac{N_{A}}{A\
\text{grams}}\right)~{},$ (8)

where $N_{A}$ is Avogadro’s number and $P_{\text{surv}}$ is the survival
probability that an incoming $\chi$ is not captured by naturally occurring
elements in the terrestrial overburden. Since we are interested in dark photon
masses that forbid $\chi$ from binding with elements much lighter than xenon,
premature capture in the terrestrial crust can only occur from binding with
rare heavy elements.

The most abundant of such elements is barium, which possesses an atomic number
and mass slightly greater than xenon. The average density of barium in Earth’s
crust is $n_{\text{Ba}}\sim${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$$ with an
exact value that varies geographically. For our estimates below, we adopt a
range $n_{\text{Ba}}\in(0.7-6)\times${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$$,
whose log-central value of $2\times${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$$ is
representative of limestone Wedepohl (2014), the dominant rock in the
overburden at Gran Sasso National Laboratory. The survival probability is then
given by
$P_{\text{surv}}\simeq\text{exp}[-n_{\text{Ba}}\,\sigma_{\text{cap}}v_{\text{rel}}\,\Delta
t]$, where $\sigma_{\text{cap}}$ is the cross section for $\chi$ to bind with
barium and $\Delta t$ is the time it takes for $\chi$ to penetrate the XENON1T
overburden, which consists of $h\simeq 1.4\ \text{km}$ of rock. The transit
time is given by the minimum of either the diffusion or gravitational-drift
timescales, $\Delta t\sim\min(t_{\text{diff}},t_{\text{drift}})$, where the
diffusion time $t_{\text{diff}}$ and the gravitational-drift timescale
$t_{\text{drift}}$ are approximately

$t_{\text{diff}}\sim
h^{2}\,\frac{\Gamma_{\text{drag}}\,m_{\chi}}{$300\text{\,}\mathrm{K}$}~{}~{},~{}~{}t_{\text{drift}}\sim
h\,\frac{\Gamma_{\text{drag}}}{g}~{},$ (9)

and $\Gamma_{\text{drag}}\sim(\mu/m_{\chi})n_{N}\sigma_{T}v_{\text{rel}}$ is
the drag rate (the inverse timescale for a particle to change its momentum by
an $\mathcal{O}(1)$ amount) from elastic scattering off of light nuclei with
density $n_{N}$ in the crust. Fixing $E_{B}=2.5\ \text{keV}$ in xenon and
$m_{A^{\prime}}=50\ \text{MeV}$, we find that diffusion is efficient
($t_{\text{diff}}\ll t_{\text{drift}}$) for
$m_{\chi}\ll$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$$, while gravitational
drift is most efficient ($t_{\text{drift}}\ll t_{\text{diff}}$) for
$m_{\chi}\gg$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$$, such that
$P_{\text{surv}}\simeq 1$ in either case. Instead, for DM masses of
$m_{\chi}\sim$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$$, $P_{\text{surv}}\ll
1$ and premature capture by naturally occurring barium may exponentially
reduce the DM flux that is able to reach underground detectors.777Even if
capture by barium is significant, a flux of DM particles can be inadvertently
delivered to a lab through other means, such as air ventilation, bypassing the
need to diffuse through the rock. We conservatively neglect such processes in
our analysis.

Figure 3: (Left) Limits derived from XENON1T (green) on the fractional
abundance $f_{\chi}$, in the case that the binding energy between $\chi$ and
xenon is $E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ and
$m_{A^{\prime}}=$50\text{\,}\mathrm{Me\kern-0.50003ptV}$$. The solid green
line assumes an abundance of barium in the Gran Sasso overburden of
$2\times${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$$, whereas the upper and lower
dashed green lines assume a barium density that is larger or smaller by a
factor of three, respectively. Also shown are existing constraints (shaded
gray) from searches for $\chi$-nuclear elastic scattering at CRESST Angloher
_et al._ (2017), CDMS Abrams  _et al._ (2002), and XENON1T Aprile  _et al._
(2017). (Right) Contours of $(\chi\,\text{Xe}^{+})$ binding energy
$$2\text{\,}\mathrm{ke\kern-0.50003ptV}$\leq
E_{B}\leq$3\text{\,}\mathrm{ke\kern-0.50003ptV}$$ for various choices of
$\alpha_{D}$, fixing $m_{\chi}=1\ \text{TeV}$, in the
$\\{m_{A^{\prime}},\epsilon\\}$ plane. Along the dashed parts of the contours,
$\chi$ also binds with Fe. The gray regions are currently excluded by
accelerator and beam dump searches for visibly decaying dark photons
Lanfranchi  _et al._ (2020). Future acclerator searches, for instance APEX
Abrahamyan  _et al._ (2011), FASER Ariga  _et al._ (2019), HPS Adrian  _et
al._ (2018), LHCb Ilten  _et al._ (2015, 2016), NA64 Gninenko (2014); Andreas
_et al._ (2013), and SeaQuest Berlin  _et al._ (2018), will explore nearly all
of the currently allowed parameter space shown.

The time dependence of these signals depends on the terrestrial $\chi$
density, which scales linearly with the velocity of the virialized galactic
population as $n_{\chi}\propto v_{\text{vir}}$. Since in Earth’s frame
$v_{\text{vir}}$ varies due to the relative motion of Earth around the Sun, we
expect an annual modulation of $n_{\chi}$. Furthermore, the short mean free
path in Earth’s crust implies a sizeable diurnal modulation. For the signals
discussed in this work, modulation over a timescale $t_{\text{mod}}$ is not
present if the diffusion time for $\chi$ to travel a length
$v_{\text{drift}}\,t_{\text{mod}}$ is shorter than $t_{\text{mod}}$. We
estimate that for masses $m_{\chi}\ll 10\ \text{GeV}$ or $m_{\chi}\ll 100\
\text{GeV}$, diffusion is strong enough to wash out an annual or daily
modulation in XENON1T data, respectively. Since, unlike a WIMP, $\chi$ does
not free stream throughout Earth, the phase of the annual modulation depends
on the relative direction between the vector normal to Earth’s surface at Gran
Sasso and Earth’s galactic motion. As a result, we find that the annual
modulation for such strongly interacting relics peaks near the spring and fall
equinox, which is consistent with the findings of Ref. Laletin and Cudell
(2019).

In Fig. 3 (left), we illustrate the existing sensitivity of XENON1T to
fractional abundances $f_{\chi}$ as a function of $m_{\chi}$, fixing
$m_{A^{\prime}}=50\ \text{MeV}$ and
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ in xenon. Following Ref.
Bloch  _et al._ (2020), an S1-S2 signal rate of $R_{\text{sig}}\simeq 60\
\text{tonne}^{-1}\ \text{yr}^{-1}$ accounts for the excess of electron recoil
events as measured by XENON1T. This model space is shown as the green lines in
Fig. 3 (left) for various assumptions regarding the possible barium abundance
in the Gran Sasso overburden (from bottom to top, these contours assume a
barium density of $0.7\times${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$$,
$2\times${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$$, and
$6\times${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$$, respectively). Thus, values
of $f_{\chi}$ that lie significantly above these lines are presently excluded
by XENON1T. We see that if a component of DM can bind to heavy nuclei, XENON1T
probes fractional abundances as small as $f_{\chi}\sim 10^{-18}$. Previously
existing constraints on this parameter space are also shown as shaded gray
regions. These are derived from searches at a surface-level CRESST run
Angloher  _et al._ (2017), as well as underground CDMS Abrams  _et al._ (2002)
and XENON1T Aprile  _et al._ (2017) runs that are sensitive to the nuclear
elastic scattering of $\chi$ particles that do not thermalize with the
terrestrial environment. For small $\chi$ masses, these limits are
significantly weakened due to the increased likelihood for $\chi$ to
thermalize in the overburden, and thus are typically many orders of magnitude
weaker than those derived from bound state formation in XENON1T.

We also note that a search for anomalously heavy gold atoms (in this case, due
to heavy $\chi$ particles bound to the nucleus) can be recast as constraints
in this model space Javorsek  _et al._ (2002). Conservatively assuming that
the gold nuggets used in these tests originated below the point of
thermalization for $\chi$ throughout the entire age of the Earth, and that the
capture rate in gold is similar to that of xenon, the corresponding limits are
subdominant to the direct detection constraints described above. In
particular, we estimate that these searches are sensitive to fractional
abundances of $f_{\chi}\gtrsim 10^{-9}$ for
$m_{\chi}\sim$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$$ and $f_{\chi}\gtrsim
10^{-6}$ for $m_{\chi}\sim$\text{\,}\mathrm{Te\kern-0.50003ptV}$$, with no
sensitivity outside this mass range.

_Discussion._ We have explored the intriguing possibility that the excess of
events at XENON1T is explained by an exponentially small subcomponent of DM
that binds with heavy nuclei through a dark photon mediator of mass
$m_{A^{\prime}}\sim\mathcal{O}(10-100)\ \text{MeV}$ that generates a sizeable
DM-SM coupling $q_{\text{eff}}=\epsilon\sqrt{\alpha_{D}/\alpha}\gtrsim
10^{-3}$. For a perturbative model in which
$\alpha_{D}\lesssim\mathcal{O}(1)$, we are thus motivated to consider
$\epsilon\gtrsim 10^{-4}$. This parameter space can be efficiently probed by
searching for visible decays of dark photons that are produced in the
collisions of dedicated accelerator experiments.888For $m_{A^{\prime}}\ll 1\
\text{MeV}$, accelerator searches for GeV-scale $\chi$ particles exclude
$q_{\text{eff}}\gtrsim 10^{-2}$ Acciarri  _et al._ (2020), but lack
sensitivity to models with much larger $A^{\prime}$ masses, as considered in
this work. This is illustrated in Fig. 3 (right) in the
$\\{m_{A^{\prime}},\varepsilon\\}$ plane. For two representative choices of
$\alpha_{D}$ shown as red and blue bands, the $\chi-\text{Xe}$ binding energy
is varied from $(2-3)\ \text{keV}$, while fixing
$m_{\chi}=$1\text{\,}\mathrm{Te\kern-0.50003ptV}$$. Also shown in dark gray
are existing constraints from searches for the production and visible decay of
dark photons in accelerator and beam dump experiments Lanfranchi  _et al._
(2020). Near future accelerator experiments, such as APEX Abrahamyan  _et al._
(2011), FASER Ariga  _et al._ (2019), HPS Adrian  _et al._ (2018), LHCb Ilten
_et al._ (2015, 2016), NA64 Gninenko (2014); Andreas  _et al._ (2013), and
SeaQuest Berlin  _et al._ (2018), will be able to decisively probe this region
of parameter space Agrawal  _et al._ (2021).

We have focused on exhibiting the salient features of how DM-SM bound state
formation can give rise to signals in the XENON1T S1-S2 data of Ref. Aprile
_et al._ (2020) since this search has the largest exposure among experiments
sensitive to $\mathcal{O}(\text{keV})$ electron recoils. However, other
experiments may have sensitivity to new parameter space for two reasons.
First, when threshold energies are reduced to below $\sim 1\ \text{keV}$ at
the expense of reduced exposure or larger backgrounds, smaller couplings
$q_{\text{eff}}$ can be probed for larger $f_{\chi}$. Examples in this vein
are the XENON1T S2-only Aprile  _et al._ (2019) and XENON10 Angle  _et al._
(2011) searches, which are sensitive to energy depositions as small as
$\mathcal{O}(100)\ \text{eV}$ and $\mathcal{O}(10)\ \text{eV}$, respectively.
Second, nuclear targets which are heavier than xenon are sensitive to smaller
$q_{\text{eff}}$. Examples among these are tungsten in the CRESST experiment
Angloher  _et al._ (2017) and thallium dopant present in the DAMA experiment
Bernabei  _et al._ (2008, 2010, 2018). However, since CRESST vetoes on
electron recoil events and reports limits only on nuclear recoil events, its
public results are at present insensitive to this model.

The DAMA experiment has reported a long-standing annual modulation signal in
the $(1-6)\ \text{keV}$ energy range Bernabei  _et al._ (2008, 2010, 2018).
While the primary target materials (sodium and iodine) are lighter than xenon,
thallium is present at the $\sim 10^{-3}$ level. Since thallium nuclei are
much heavier than xenon, this raises the tantalizing possibility that
$q_{\text{eff}}$ is large enough such that $\chi$ binds deeply in thallium but
is too small to deposit a significant energy above threshold in xenon-based
experiments, thus explaining the DAMA observation. For $m_{\chi}\gtrsim 100\
\text{GeV}$, explaining the modulated rate at DAMA requires fractional
abundances of $f_{\chi}\gtrsim 10^{-9}$, which is in conflict with limits
derived from CDMS Abrams  _et al._ (2002) (see Fig. 3 (left)). Instead, for
$m_{\chi}\lesssim 100\ \text{GeV}$, $f_{\chi}\sim 10^{-12}-10^{-9}$ predicts a
rate consistent with DAMA, provided that $q_{\text{eff}}$ lies in a narrow
range to facilitate a $\chi-\text{Tl}$ binding energy of $E_{B}\sim 2\
\text{keV}$ without being in tension with an S2-only search at XENON1T Aprile
_et al._ (2019). However, as discussed above, the phase of the annual
modulation of these signals is inconsistent with DAMA’s observation of peak
rates near the beginning of June and December. Although strongly interacting
relics cannot accommodate this excess of events, it is interesting to note
that DAMA can place some of the most stringent constraints on models in which
$\chi$ preferentially binds to thallium (and heavier elements).

_Acknowledgements._ The authors would like to acknowledge Carlos Blanco, Oren
Slone, and Neal Weiner for fruitful conversations. AB is supported by the
James Arthur Fellowship. HL is supported by the DOE under Award Number DE-
SC0007968, NSF grant PHY-1915409, and the Simons Foundation. MP is supported
in part by U.S. Department of Energy (Grant No. DE-SC0011842). HR is supported
in part by NSF Grant PHY-1720397 and the Gordon and Betty Moore Foundation
Grant GBMF7946.

## References

  * Lin (2019) T. Lin, PoS 333, 009 (2019), arXiv:1904.07915 [hep-ph] .
  * Pospelov and Ritz (2007) M. Pospelov and A. Ritz, Phys. Lett. B 651, 208 (2007), arXiv:hep-ph/0703128 .
  * An  _et al._ (2012) H. An, M. Pospelov, and J. Pradler, Phys. Rev. Lett. 109, 251302 (2012), arXiv:1209.6358 [hep-ph] .
  * Fornal  _et al._ (2020) B. Fornal, B. Grinstein, and Y. Zhao, Phys. Lett. B 811, 135869 (2020), arXiv:2005.04240 [hep-ph] .
  * Wallemacq (2013) Q. Wallemacq, Phys. Rev. D 88, 063516 (2013), arXiv:1307.7623 [astro-ph.CO] .
  * Wallemacq (2014) Q. Wallemacq, Adv. High Energy Phys. 2014, 525208 (2014), arXiv:1401.5243 [hep-ph] .
  * Wallemacq and Cudell (2015) Q. Wallemacq and J.-R. Cudell, JCAP 02, 011 (2015), arXiv:1411.3178 [hep-ph] .
  * Laletin and Cudell (2019) M. Laletin and J.-R. Cudell, JCAP 07, 010 (2019), arXiv:1903.04637 [hep-ph] .
  * Bernabei  _et al._ (2008) R. Bernabei _et al._ (DAMA), Eur. Phys. J. C 56, 333 (2008), arXiv:0804.2741 [astro-ph] .
  * Bernabei  _et al._ (2010) R. Bernabei _et al._ (DAMA, LIBRA), Eur. Phys. J. C 67, 39 (2010), arXiv:1002.1028 [astro-ph.GA] .
  * Bernabei  _et al._ (2018) R. Bernabei _et al._ , Nucl. Phys. Atom. Energy 19, 307 (2018), arXiv:1805.10486 [hep-ex] .
  * Boehm and Fayet (2004) C. Boehm and P. Fayet, Nucl. Phys. B 683, 219 (2004), arXiv:hep-ph/0305261 .
  * Arkani-Hamed  _et al._ (2009) N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, Phys. Rev. D 79, 015014 (2009), arXiv:0810.0713 [hep-ph] .
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Lett. B 662, 53 (2008), arXiv:0711.4866 [hep-ph] .
  * Pospelov  _et al._ (2020) M. Pospelov, S. Rajendran, and H. Ramani, Phys. Rev. D 101, 055001 (2020), arXiv:1907.00011 [hep-ph] .
  * Pospelov and Ramani (2020) M. Pospelov and H. Ramani, (2020), arXiv:2012.03957 [hep-ph] .
  * Neufeld  _et al._ (2018) D. A. Neufeld, G. R. Farrar, and C. F. McKee, Astrophys. J. 866, 111 (2018), arXiv:1805.08794 [astro-ph.CO] .
  * Essig  _et al._ (2017) R. Essig, T. Volansky, and T.-T. Yu, Phys. Rev. D 96, 043017 (2017), arXiv:1703.00910 [hep-ph] .
  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Aprile  _et al._ (2020) E. Aprile _et al._ (XENON), Phys. Rev. D 102, 072004 (2020), arXiv:2006.09721 [hep-ex] .
  * Alonso-Álvarez  _et al._ (2020) G. Alonso-Álvarez, F. Ertas, J. Jaeckel, F. Kahlhoefer, and L. J. Thormaehlen, JCAP 11, 029 (2020), arXiv:2006.11243 [hep-ph] .
  * An  _et al._ (2020) H. An, M. Pospelov, J. Pradler, and A. Ritz, (2020), arXiv:2006.13929 [hep-ph] .
  * Bloch  _et al._ (2020) I. M. Bloch, A. Caputo, R. Essig, D. Redigolo, M. Sholapurkar, and T. Volansky, (2020), arXiv:2006.14521 [hep-ph] .
  * McKeen  _et al._ (2020) D. McKeen, M. Pospelov, and N. Raj, Phys. Rev. Lett. 125, 231803 (2020), arXiv:2006.15140 [hep-ph] .
  * Farzan and Rajaee (2020) Y. Farzan and M. Rajaee, (2020), arXiv:2007.14421 [hep-ph] .
  * Brdar  _et al._ (2021) V. Brdar, A. Greljo, J. Kopp, and T. Opferkuch, JCAP 01, 039 (2021), arXiv:2007.15563 [hep-ph] .
  * Kuo  _et al._ (2021) J.-L. Kuo, M. Pospelov, and J. Pradler, (2021), arXiv:2102.08409 [hep-ph] .
  * Baryakhtar  _et al._ (2020) M. Baryakhtar, A. Berlin, H. Liu, and N. Weiner, (2020), arXiv:2006.13918 [hep-ph] .
  * Bramante and Song (2020) J. Bramante and N. Song, Phys. Rev. Lett. 125, 161805 (2020), arXiv:2006.14089 [hep-ph] .
  * Harigaya  _et al._ (2020) K. Harigaya, Y. Nakai, and M. Suzuki, Phys. Lett. B 809, 135729 (2020), arXiv:2006.11938 [hep-ph] .
  * Holdom (1986) B. Holdom, Phys. Lett. B 166, 196 (1986).
  * Pospelov and Ritz (2009) M. Pospelov and A. Ritz, Phys. Lett. B 671, 391 (2009), arXiv:0810.1502 [hep-ph] .
  * Feng  _et al._ (2010) J. L. Feng, M. Kaplinghat, and H.-B. Yu, Phys. Rev. D 82, 083525 (2010), arXiv:1005.4678 [hep-ph] .
  * Cirelli  _et al._ (2017) M. Cirelli, P. Panci, K. Petraki, F. Sala, and M. Taoso, JCAP 05, 036 (2017), arXiv:1612.07295 [hep-ph] .
  * An  _et al._ (2016) H. An, M. B. Wise, and Y. Zhang, Phys. Rev. D 93, 115020 (2016), arXiv:1604.01776 [hep-ph] .
  * Bargmann (1952) V. Bargmann, Proceedings of the National Academy of Sciences 38, 961 (1952), https://www.pnas.org/content/38/11/961.full.pdf .
  * Miramonti (2005) L. Miramonti, AIP Conf. Proc. 785, 3 (2005), arXiv:hep-ex/0503054 .
  * Colquhoun  _et al._ (2021) B. Colquhoun, S. Heeba, F. Kahlhoefer, L. Sagunski, and S. Tulin, Phys. Rev. D 103, 035006 (2021), arXiv:2011.04679 [hep-ph] .
  * Berlin  _et al._ (2022) A. Berlin, H. Liu, M. Pospelov, and H. Ramani, (2022), To appear.
  * Bell  _et al._ (2020) N. F. Bell, J. B. Dent, J. L. Newstead, S. Sabharwal, and T. J. Weiler, Phys. Rev. D 101, 015012 (2020), arXiv:1905.00046 [hep-ph] .
  * Ibe  _et al._ (2018) M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, JHEP 03, 194 (2018), arXiv:1707.07258 [hep-ph] .
  * Wedepohl (2014) K. H. Wedepohl,  _Handbook of Geochemistry, Volume II/4: Elements Kr(36) to Ba(56)_ (Springer Berlin, Berlin, 2014).
  * Angloher  _et al._ (2017) G. Angloher _et al._ (CRESST), Eur. Phys. J. C 77, 637 (2017), arXiv:1707.06749 [astro-ph.CO] .
  * Abrams  _et al._ (2002) D. Abrams _et al._ (CDMS), Phys. Rev. D 66, 122003 (2002), arXiv:astro-ph/0203500 .
  * Aprile  _et al._ (2017) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 119, 181301 (2017), arXiv:1705.06655 [astro-ph.CO] .
  * Lanfranchi  _et al._ (2020) G. Lanfranchi, M. Pospelov, and P. Schuster,  (2020), 10.1146/annurev-nucl-102419-055056, arXiv:2011.02157 [hep-ph] .
  * Abrahamyan  _et al._ (2011) S. Abrahamyan _et al._ (APEX), Phys. Rev. Lett. 107, 191804 (2011), arXiv:1108.2750 [hep-ex] .
  * Ariga  _et al._ (2019) A. Ariga _et al._ (FASER), Phys. Rev. D 99, 095011 (2019), arXiv:1811.12522 [hep-ph] .
  * Adrian  _et al._ (2018) P. H. Adrian _et al._ (HPS), Phys. Rev. D 98, 091101 (2018), arXiv:1807.11530 [hep-ex] .
  * Ilten  _et al._ (2015) P. Ilten, J. Thaler, M. Williams, and W. Xue, Phys. Rev. D 92, 115017 (2015), arXiv:1509.06765 [hep-ph] .
  * Ilten  _et al._ (2016) P. Ilten, Y. Soreq, J. Thaler, M. Williams, and W. Xue, Phys. Rev. Lett. 116, 251803 (2016), arXiv:1603.08926 [hep-ph] .
  * Gninenko (2014) S. N. Gninenko, Phys. Rev. D 89, 075008 (2014), arXiv:1308.6521 [hep-ph] .
  * Andreas  _et al._ (2013) S. Andreas _et al._ , (2013), arXiv:1312.3309 [hep-ex] .
  * Berlin  _et al._ (2018) A. Berlin, S. Gori, P. Schuster, and N. Toro, Phys. Rev. D 98, 035011 (2018), arXiv:1804.00661 [hep-ph] .
  * Javorsek  _et al._ (2002) D. Javorsek, D. Elmore, E. Fischbach, D. Granger, T. Miller, D. Oliver, and V. Teplitz, Phys. Rev. D 65, 072003 (2002).
  * Acciarri  _et al._ (2020) R. Acciarri _et al._ (ArgoNeuT), Phys. Rev. Lett. 124, 131801 (2020), arXiv:1911.07996 [hep-ex] .
  * Agrawal  _et al._ (2021) P. Agrawal _et al._ , (2021), arXiv:2102.12143 [hep-ph] .
  * Aprile  _et al._ (2019) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 251801 (2019), arXiv:1907.11485 [hep-ex] .
  * Angle  _et al._ (2011) J. Angle _et al._ (XENON10), Phys. Rev. Lett. 107, 051301 (2011), [Erratum: Phys.Rev.Lett. 110, 249901 (2013)], arXiv:1104.3088 [astro-ph.CO] .
  * Landau and Lifshitz (1981) L. D. Landau and L. M. Lifshitz,  _Quantum Mechanics Non-Relativistic Theory, Third Edition: Volume 3_, 3rd ed. (Butterworth-Heinemann, 1981).
  * Bunge  _et al._ (1993) C. F. Bunge, J. A. Barrientos, and A. V. Bunge, Atom. Data Nucl. Data Tabl. 53, 113 (1993).
  * McLean and McLean (1981) A. D. McLean and R. S. McLean, Atom. Data Nucl. Data Tabl. 26, 197 (1981).
  * Xeni (2018) C. Xeni, “Sensitivity of dark matter direct detection experiments to a cosmological relativistic background through electron recoils,”  (2018).
  * Cardona and Ley (1978) M. Cardona and L. Ley,  _Photoemission in solids I_ , Vol. 26 (1978).
  * Sakurai and Napolitano (2020) J. J. Sakurai and J. Napolitano,  _Modern Quantum Mechanics_, Quantum physics, quantum information and quantum computation (Cambridge University Press, 2020).
  * West and Morton (1978) J. B. West and J. Morton, Atom. Data Nucl. Data Tabl. 22, 103 (1978).
  * Veigele (1973) W. J. Veigele, Atom. Data Nucl. Data Tabl. 5, 51 (1973).

oneΔ

Low-Energy Signals from the Formation of Dark Matter-Nuclear Bound States

Supplemental Material
Asher Berlin, Hongwan Liu, Maxim Pospelov, and Harikrishnan Ramani

In this Supplemental Material, we provide a detailed derivation of the cross
section for bound state formation $\chi+A\to(\chi A^{+})+e^{-}$.

## Perturbation Theory Set-Up

The Hamiltonian for the system can be written as
$\hat{H}=\hat{H}_{0}+\hat{V}_{\text{pert}}$, with $\hat{V}_{\text{pert}}$
being a small perturbation to the unperturbed Hamiltonian $\hat{H}_{0}$. The
cross section for transitioning from an initial state $\ket{i}$ to a continuum
final state $\ket{f}$, both eigenstates of $\hat{H}_{0}$, is given by Fermi’s
golden rule. In the process we are considering, an electron is ejected into a
continuum final state, with a change in energy given by $E_{B}$, the binding
energy of $(\chi A)$, so that Fermi’s golden rule reads Landau and Lifshitz
(1981)

$\displaystyle\textrm{d}\sigma=\frac{\textrm{d}^{3}\mathbf{p}_{e}}{(2\pi)^{3}}\,2\pi\,\left|\braket{f}{\hat{V}_{\text{pert}}}{i}\right|^{2}\delta\left(E_{e}-E_{B}+\omega_{e,i}\right)=\frac{\textrm{d}\Omega_{e}}{(2\pi)^{2}}\,m_{e}p_{e}\,\left|\braket{f}{\hat{V}_{\text{pert}}}{i}\right|^{2}\,,$
(S1)

where $\mathbf{p}_{e}$ and $E_{e}$ are the momentum and energy of the final
state electron respectively, and $\omega_{e,i}$ is the binding energy of the
initial state electron. In the second equality, we have integrated over the
delta function, which leaves just the solid angle of the electron momentum
$\Omega_{e}$ and fixes $p_{e}=\sqrt{2m_{e}(E_{B}-\omega_{e,i})}$. We neglect
the kinetic energy of the incoming thermalized $\chi$ particle, which is much
smaller than $E_{B}$. The appropriate normalization of the states $\ket{i}$
and $\ket{f}$ will be discussed below.

The total Hamiltonian involving the nucleus $N$, the DM particle $\chi$, and
the ejected electron $e$ contains the following terms:

$\displaystyle\hat{H}=\hat{T}+\hat{V}_{\chi N}+\hat{V}_{eA}\,,$ (S2)

where $\hat{T}$ is the kinetic energy terms of all the particles involved, and
the various potential energy contributions are $\hat{V}_{\chi N}$ between
$\chi$ and the nucleus $N$ and $\hat{V}_{eA}$ between the ejected electron and
the atom (including both $N$ and the other unejected electrons).999The other
potential energy terms involving solely the unejected electrons are
unimportant, since either the ejected electron states or the $\chi$ states are
eigenstates of these operators, and the initial and final wavefunctions are
orthogonal. Furthermore, we do not include the potential arising from
interactions between $\chi$ and electrons because $\chi$ becomes tightly bound
to the atom and thus dominantly experiences just the bare nucleus: the typical
size of the bound state, given in Eq. (S6), is much smaller than the size of
the xenon atom. Let us now discuss both of these contributions.

We model $V_{eA}$ by adopting the Thomas-Fermi model for the neutral atom,
which gives an effective screened potential energy of the nucleus and the
atomic electrons as a function of the distance to the nucleus. In this model,
the effective screened potential energy of the entire atom $V_{A}$ is Landau
and Lifshitz (1981)

$\displaystyle
V_{A}(|\mathbf{r}_{e}-\mathbf{r}_{N}|)=-\,\frac{Z\alpha}{|\mathbf{r}_{e}-\mathbf{r}_{N}|}~{}\phi\left(\frac{|\mathbf{r}_{e}-\mathbf{r}_{N}|}{b}\right)\,,$
(S3)

where $b\equiv(9\pi^{2}/2Z)^{1/3}/(4\alpha m_{e})$, and $\phi(\xi)$ satisfies
the Thomas-Fermi equation,

$\displaystyle\frac{d^{2}\phi}{d\xi^{2}}=\frac{\phi^{3/2}}{\sqrt{\xi}}\,,$
(S4)

with the boundary conditions $\phi(0)=1$ and $\lim_{\xi\to\infty}\phi(\xi)=0$,
which can be solved numerically. $V_{eA}$, the potential energy between the
ejected electron and the rest of the atom, can be obtained by subtracting the
contribution from the initial bound state electron, i.e.,

$\displaystyle
V_{eA}(|\mathbf{r}_{e}-\mathbf{r}_{N}|)=V_{A}(|\mathbf{r}_{e}-\mathbf{r}_{N}|)-\alpha\int\textrm{d}^{3}\mathbf{r}^{\prime}\frac{|\psi_{e,i}(\mathbf{r}^{\prime})|^{2}}{|(\mathbf{r}_{e}-\mathbf{r}_{N})-\mathbf{r}^{\prime}|}\,,$
(S5)

where $\psi_{e,i}(\mathbf{r})$ is the wavefunction of the initial bound state
electron before it is ejected. The explicit form of $\psi_{e,i}(r)$ will be
discussed below near Eq. (S24).

The characteristic size of the $(\chi N)$ bound state is

$\displaystyle\frac{1}{\sqrt{2\mu
E_{B}}}\simeq$9\text{\,}\mathrm{fm}$\times\left(\frac{$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$}{\mu}\right)^{1/2}\left(\frac{$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$}{E_{B}}\right)^{1/2}\,,$
(S6)

which is comparable to typical nuclear radii
$R_{\text{nuc}}\sim$5\text{\,}\mathrm{fm}$\times(A/100)^{1/3}$. Here, $\mu$ is
the reduced mass of $\chi$ and $N$. Finite size effects of the nucleus are
therefore important for describing the potential energy $V_{\chi N}$ between
$\chi$ and $N$. Modeling the nucleus as a uniformly charged sphere of radius
$R_{\text{nuc}}$ centered at $\mathbf{r}_{N}$, the potential at
$\mathbf{r}_{\chi}$ can be written as

$\displaystyle V_{\chi
N}(|\mathbf{r}_{\chi}-\mathbf{r}_{N}|)=-\int_{N}\textrm{d}^{3}\mathbf{r}^{\prime}~{}\frac{\rho}{|(\mathbf{r}_{\chi}-\mathbf{r}_{N})-\mathbf{r}^{\prime}|}~{}e^{-m_{A^{\prime}}|(\mathbf{r}_{\chi}-\mathbf{r}_{N})-\mathbf{r}^{\prime}|}~{},$
(S7)

where $\rho\equiv q_{\text{eff}}\,Z\,\alpha\,/(4\pi R_{\text{nuc}}^{3}/3)$,
and the integral is performed over the volume of the nucleus. This integral
can be evaluated analytically to give

$\displaystyle V_{\chi
N}(r)=-\,\frac{4\pi\,\rho}{m_{A^{\prime}}^{3}\,r}\begin{cases}m_{A^{\prime}}\,r-e^{-m_{A^{\prime}}R_{\text{nuc}}}(1+m_{A^{\prime}}\,R_{\text{nuc}})\,\sinh\left(m_{A^{\prime}}\,r\right)&(r<R_{\text{nuc}})\\\
e^{-m_{A^{\prime}}r}\Big{[}m_{A^{\prime}}\,R_{\text{nuc}}\cosh(m_{A^{\prime}}\,R_{\text{nuc}})-\sinh(m_{A^{\prime}}\,R_{\text{nuc}})\Big{]}&(r\geq
R_{\text{nuc}})\,,\end{cases}$ (S8)

where we have defined $\mathbf{r}\equiv\mathbf{r}_{N}-\mathbf{r}_{\chi}$.
Defining the origin to be located at the $(\chi N)$ center-of-mass,
$\mathbf{r}_{N}=(\mu/m_{N})\,\mathbf{r}$ and
$\mathbf{r}_{\chi}=-(\mu/m_{\chi})\,\mathbf{r}$. Hence, after bound state
formation, the position of $N$ is parametrically $\mu/m_{N}$ times the spatial
extent of the bound state, i.e., $r_{N}\sim(\mu/m_{N})/\sqrt{2\mu E_{B}}$.
Since this is always much smaller than the typical size of the atom, we may
expand $V_{eA}(|\mathbf{r}_{e}-\mathbf{r}_{N}|)$ in powers of the small
quantity $\mathbf{r}_{N}$ to obtain

$\displaystyle V_{eA}(|\mathbf{r}_{e}-\mathbf{r}_{N}|)$
$\displaystyle=V_{eA}(r_{e})+\mathbf{r}_{N}^{i}\,\partial^{i}\,V_{eA}(r_{e})+\frac{1}{2}\,\mathbf{r}_{N}^{i}\mathbf{r}_{N}^{j}\,\partial^{i}\partial^{j}\,V_{eA}(r_{e})+\cdots$
$\displaystyle\equiv V_{eA}^{(0)}+V_{eA}^{(1)}+V_{eA}^{(2)}+\cdots~{},$ (S9)

where $i$, $j$ denote spatial components, with repeated indices summed over.
We now take the unperturbed Hamiltonian to be

$\displaystyle\hat{H}_{0}\equiv\hat{T}+\hat{V}_{\chi N}+\hat{V}_{eA}^{(0)}\,.$
(S10)

The initial and final states of both $\chi$ and the ejected electron $e$ are
bound states or asymptotically free states that are eigenstates of
$\hat{H}_{0}$; the actual wavefunctions will be worked out below. The full
perturbative portion of the Hamiltonian is therefore
$\hat{V}_{eA}^{(1)}+\hat{V}_{eA}^{(2)}+\cdots$. We therefore define the
perturbative correction to the Hamiltonian as

$\displaystyle\hat{V}_{\text{pert}}\equiv\hat{V}_{eA}^{(1)}+\hat{V}_{eA}^{(2)}\,,$
(S11)

keeping just the first two terms in the expansion of
$V_{eA}(|\mathbf{r}_{e}-\mathbf{r}_{N}|)$. We assume that $\chi$ is captured
into the $s$-wave ground state of the potential $V_{\chi N}$, since it is the
most deeply bound; as we show below, under this assumption,
$\hat{V}_{eA}^{(1)}$ leads to a $p$-to-$s$ transition with a unit change in
angular momentum between the initial and final $\chi$ states, while
$\hat{V}_{eA}^{(2)}$ facilitates an $s$-to-$s$ transition instead, with no
change in $\chi$ angular momentum. We compute the rate for each angular state
by incorporating each contribution to $\hat{V}_{\text{pert}}$ separately,
estimating the total bound state capture cross section as
$\sigma_{B}\approx\sigma_{s}+\sigma_{p}$, where $\sigma_{s,p}$ is the cross
section for an $s$-to-$s$ or $p$-to-$s$ transition, respectively. This
treatment is appropriate as long as $\sigma_{s}\gg\sigma_{p}$ or vice-versa,
which is true across almost all relevant parameter space.

## Wavefunctions in a Central Potential

Before we investigate each transition separately, we note that $\chi$ and $e$
are both subject to a two-body central potential with the nucleus. In each
case, the Schrödinger equation is separable into an angular and a radial
component. For a central potential $V$, and a well-defined energy $E$ and
angular momentum quantum number $l$, the radial wavefunction
$R(\text{\LARGE\textcalligra{r}})$ satisfies the following equation:

$\displaystyle\frac{1}{\text{\LARGE\textcalligra{r}}^{\,2}}\,\frac{\textrm{d}}{\textrm{d}\text{\LARGE\textcalligra{r}}}\left(\text{\LARGE\textcalligra{r}}^{\,2}\,\frac{\textrm{d}R}{\textrm{d}\text{\LARGE\textcalligra{r}}}\right)-\left[\frac{l(l+1)}{\text{\LARGE\textcalligra{r}}^{\,2}}+2\mu\,V\right]R=-2\mu
E\,R~{},$ (S12)

where r is the relative coordinate between the two bodies, and $\mu$ is the
reduced mass of two bodies. For a bound state wavefunction $E<0$, whereas for
a continuum wavefunction $E\simeq k^{2}/2\mu>0$ such that $\mathbf{k}$ is the
momentum of the particle far away from the nuclear potential.

For bound wavefunctions, the full wavefunction is $R$ times the appropriate
spherical harmonic for angular momentum $l$ and azimuthal quantum number $m$.
Instead, continuum wavefunctions with definite total momentum $\mathbf{k}$ can
be expanded as Landau and Lifshitz (1981)

$\displaystyle\psi(\text{\LARGE\textcalligra{r}}\,,\theta)=\sum_{l=0}^{\infty}\,\frac{2l+1}{2k}\,i^{l}\,P_{l}(\cos\theta)\,R_{kl}(\text{\LARGE\textcalligra{r}})\,,$
(S13)

where $R_{kl}(\text{\LARGE\textcalligra{r}})$ is the solution to Eq. (S12)
with $E=k^{2}/2\mu$, $P_{l}$ is the $l^{\text{th}}$ Legendre polynomial, and
$\cos{\theta}\equiv\hat{\text{\LARGE\textcalligra{r}}\,}\cdot\hat{\mathbf{k}}$
(we have neglected the phase shift factor $\delta_{\ell}$, consistent with the
perturbative treatment used here). The constant prefactors correctly normalize
the asymptotic wavefunction to a single particle plane wave.

## s-to-s Capture Cross Section

We first consider the case of $s$-to-$s$ capture through $\hat{V}_{eA}^{(2)}$.
Although this term is higher order in $r_{N}$ than $\hat{V}_{eA}^{(1)}$, there
is no velocity suppression to this cross section, since there is no change in
$l$ between the initial and final states. From Eq. (S1), the matrix element
for this process in the position basis is

$\displaystyle\textrm{d}\sigma=\frac{d\Omega_{e}}{(2\pi)^{2}}\,m_{e}p_{e}\,\left|\int\textrm{d}^{3}\mathbf{r}~{}\psi_{\chi,f}^{*}(\mathbf{r})\,\psi_{\chi,i}(\mathbf{r})\int\textrm{d}^{3}\mathbf{r}_{e}~{}\psi_{e,f}^{*}(\mathbf{r}_{e})\,\psi_{e,i}(\mathbf{r}_{e})~{}V_{eA}^{(2)}(r_{e})\right|^{2}\,,$
(S14)

where the subscripts $i$ and $f$ denote initial and final states for the
position-space wavefunctions $\psi$ for $\chi$ and $e$. We can simplify this
integral with the following identity that we will utilize several times in
this Supplemental Material:

$\displaystyle\int\textrm{d}^{3}\mathbf{x}~{}f(x)\,\mathbf{x}^{i}\mathbf{x}^{j}=\frac{1}{3}\int\textrm{d}^{3}\mathbf{x}~{}f(x)\,x^{2}\,\delta^{ij}\,,$
(S15)

for some function $f(x)$. Since the $s$-wave initial and final $\chi$ states
are spherically symmetric, in Eq. (S14) we can therefore replace

$\displaystyle
V_{eA}^{(2)}=\frac{1}{2}\,\mathbf{r}_{N}^{i}\mathbf{r}_{N}^{j}\,\partial^{i}\partial^{j}\,V_{eA}(r_{e})\to\frac{r_{N}^{2}}{6}\,\nabla^{2}V_{eA}(r_{e})\,.$
(S16)

By Gauss’s law, we can relate $\nabla^{2}V_{eA}$ to the charge density of the
atom, such that

$\displaystyle
V_{eA}^{(2)}(\mathbf{r}_{e})\to\frac{2\pi}{3}\,\alpha\,r_{N}^{2}\big{(}-Z\delta^{3}(\mathbf{r}_{e})+n_{\text{TF}}(\mathbf{r}_{e})-|\psi_{e,i}(\mathbf{r}_{e})|^{2}\big{)}\,,$
(S17)

where $n_{\text{TF}}(\mathbf{r}_{e})-|\psi_{e,i}(\mathbf{r}_{e})|^{2}$ is the
number density of electrons in the Thomas-Fermi model minus the number density
of the initial bound state electron, and the Dirac delta function arises from
the nucleus, which we take to be a point charge from the point of view of the
ejected electron. The expression for the electronic number density
$n_{\text{TF}}$ is Landau and Lifshitz (1981)

$\displaystyle
n_{\text{TF}}(\mathbf{r}_{e})=\frac{1}{3\pi^{2}}\left[\frac{2\,Z\alpha\,m_{e}}{r_{e}}~{}\phi\left(\frac{r_{e}}{b}\right)\right]^{3/2}\,,$
(S18)

where $\phi$ is determined numerically from Eq. (S4). With this, Eq. (S14)
becomes

$\displaystyle\textrm{d}\sigma$
$\displaystyle=\frac{\textrm{d}\Omega_{e}}{(2\pi)^{2}}\,m_{e}p_{e}\left(\frac{2\pi\,\alpha}{3}\right)^{2}\,\left|Z\,\psi_{e,f}^{*}(0)\,\psi_{e,i}(0)-\int\textrm{d}^{3}\mathbf{r}_{e}~{}\psi_{e,f}^{*}(\mathbf{r}_{e})\,\psi_{e,i}(\mathbf{r}_{e})\,\big{(}n_{\text{TF}}(\mathbf{r}_{e})-|\psi_{e,i}(\mathbf{r}_{e})|^{2}\big{)}\right|^{2}$
$\displaystyle\qquad\qquad\qquad\qquad\quad~{}~{}\times\left|\int\textrm{d}^{3}\mathbf{r}~{}\Big{(}\frac{\mu}{m_{N}}\Big{)}^{2}r^{2}\,\psi_{\chi,f}^{*}(\mathbf{r})\,\psi_{\chi,i}(\mathbf{r})\right|^{2}\,,$
(S19)

where we have used $r_{N}=(\mu/m_{N})\,r$.

We now turn our attention to obtaining the wavefunctions in Eq. (s-to-s
Capture Cross Section), beginning with the initial and final state $\chi$
wavefunctions. These wavefunctions are eigenstates of $\hat{T}+\hat{V}_{\chi
N}$, which are naturally obtained in the $\chi$-$N$ center-of-mass frame using
the variables defined above:
$\mathbf{r}\equiv\mathbf{r}_{N}-\mathbf{r}_{\chi}$, such that
$\mathbf{r}_{N}=(\mu/m_{N})\,\mathbf{r}$ and
$\mathbf{r}_{\chi}=-(\mu/m_{\chi})\,\mathbf{r}$. The initial state $\chi$
wavefunction $\psi_{\chi,i}$ is an $l=0$ state with energy $E=k^{2}/2\mu$ and
momentum $k$ (in the $\chi$-$N$ frame); based on the decomposition of
continuum wavefunctions in a central potential given in Eq. (S13), we have

$\displaystyle\psi_{\chi,i}(r)=\frac{R_{k0}^{\chi}(r)}{2k\,\sqrt{v_{\text{rel}}}}\,,$
(S20)

where $v_{\text{rel}}$ is the velocity of the incoming $\chi$ particle in the
$\chi$-$N$ center-of-mass frame, and $1/\sqrt{v_{\text{rel}}}$ normalizes
$\psi_{\chi,i}$ to the “one particle unit current density” prescription, which
is necessary to ensure that the right-hand-side of Eq. (S1) is correctly
scaled to give the cross section on the left-hand-side Landau and Lifshitz
(1981). $R_{k0}^{\chi}(r)$ is the solution to Eq. (S12) with $V=V_{\chi N}$,
$\mu$ the $\chi-N$ reduced mass, and $E=k^{2}/2\mu\simeq 0$ (we neglect the
kinetic energy of the incoming $\chi$, which is always small relative to the
binding energy $E_{B}$). Asymptotically far away from the nucleus,
$R_{k0}^{\chi}$ should tend to the free continuum solution of Eq. (S12) (i.e.,
the solution with $V=0$ and $l=0$), which is $2kj_{0}(kr)$ Landau and Lifshitz
(1981). In the $k\ll 1/r$ limit, we therefore expect
$\lim_{r\to\infty}\psi_{\chi,i}(r)=1/\sqrt{v_{\text{rel}}}$. To incorporate
the effect of the potential $V_{\chi N}$ on $\psi_{\chi,i}(r)$ at smaller
radii, we parametrize the general solution as

$\psi_{\chi,i}(r)=G(r)/\sqrt{v_{\text{rel}}}~{}~{},~{}~{}\lim_{r\to\infty}G(r)=1~{},$
(S21)

for some dimensionless function $G(r)$ that we compute numerically by solving
Eq. (S12) for the $\chi-N$ system.

The wavefunction of the final state $\chi$, $\psi_{\chi,f}$, is a bound
eigenstate of the potential $V_{\chi N}$ with $l=0$; we obtain the radial
component $R_{\text{b.s.}}^{\chi}(r)$ of $\psi_{\chi,f}$ by again numerically
solving Eq. (S12) for the $\chi-N$ system but with $E=-E_{B}$ and adjusting
the effective coupling $q_{\text{eff}}$ until we obtain a solution that goes
to zero as $r\to\infty$ with no nodes. Including the angular piece of the
wavefunction (i.e., the spherical harmonic $Y_{00}=1/\sqrt{4\pi}$), we have

$\displaystyle\psi_{\chi,f}(r)=R_{\text{b.s.}}^{\chi}(r)/\sqrt{4\pi}\,.$ (S22)

Next, we consider the initial and final state wavefunctions for the ejected
electron. From Eq. (s-to-s Capture Cross Section), we see that the capture
cross section is enhanced by $Z^{2}$ for electronic wavefunctions that have
non-vanishing weight at the origin. Since wavefunctions of angular momentum
$\ell$ scale as $\sim r^{\ell}$ at small radii, we focus on $s$-wave
electronic wavefunctions. The electron is treated as initially occupying a
bound eigenstate of $\hat{T}+\hat{V}_{eA}^{(0)}$. For simplicity, we adopt the
Roothaan-Hartree-Fock electronic wavefunctions for the atomic orbitals, where
the radial part of the wavefunctions is decomposed into a linear combination
of Slater orbitals, as computed in Refs. Bunge  _et al._ (1993); McLean and
McLean (1981) (see also Ref. Xeni (2018)). This decomposition can be written
as

$\displaystyle
R_{nl}^{e}(r_{e})=\sum_{j}C_{jnl}\,S_{jl}(r_{e})~{}~{},~{}~{}S_{jl}(r_{e})=\frac{(2Z_{jl})^{n_{jl}+1/2}}{\sqrt{(2n_{jl})!}}\,a_{0}^{-3/2}\left(\frac{r_{e}}{a_{0}}\right)^{n_{jl}-1}e^{-Z_{jl}\,r_{e}/a_{0}}~{},$
(S23)

where $C_{jnl}$ is the weight given to each Slater orbital (indexed by $j$),
$Z_{jl}$ is an effective charge, $n_{jl}$ is the principal quantum number of
that Slater orbital (see, e.g., Table 4.1 of Ref. Xeni (2018)), and
$a_{0}=1/\alpha m_{e}$ is the Bohr radius. For a $\chi-\text{Xe}$ binding
energy of $E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$, only the $3s$,
$4s$, and $5s$ electrons in xenon are shallowly bound enough to be ejected in
an $s$-to-$s$ transition. The full initial state wavefunction of the electron
for each of these states is then simply

$\displaystyle\psi_{e,i}=\frac{R_{n0}^{e}(r_{e})}{\sqrt{4\pi}}\,.$ (S24)

after including the appropriate spherical harmonic.

For the outgoing electron, we obtain the radial wavefunction
$R_{p_{e}0}^{e}(r_{e})$ by numerically solving Eq. (S12) for the $e-A$ system,
i.e., with $\mu=m_{e}$, $V=V_{eA}^{(0)}$, $l=0$, and
$E=p_{e}^{2}/2m_{e}=E_{B}-\omega_{e,i}$. Asymptotically far away from the
nucleus, we again expect $R_{p_{e}0}^{e}(r_{e})$ to approach the free
continuum solution $2p_{e}j_{0}(p_{e}r)$ up to a phase, which sets the
normalization of $R^{e}_{p_{e}0}(r_{e})$. The full final electron wavefunction
is then (see Eq. (S13) for the normalization factor)

$\displaystyle\psi_{e,f}(r_{e})=\frac{R^{e}_{p_{e}0}(r_{e})}{2p_{e}}\,.$ (S25)

Having determined the wavefunctions in Eq. (s-to-s Capture Cross Section), we
are now ready to compute $\sigma_{s}$, the $s$-to-$s$ capture cross section
after summing over all possible initial electron states. Returning to Eq.
(s-to-s Capture Cross Section), we find that the contribution to
$V_{eA}^{(2)}$ from the screening electrons (the second term in the first set
of vertical brackets) is subdominant (at the level of $\sim 0.5\%$) to the
$Z^{2}$ contribution from the nucleus itself (the first term in the first set
of vertical brackets). Keeping only this $Z^{2}$ contribution yields

$\displaystyle\textrm{d}\sigma_{s}$
$\displaystyle=2\,\sum_{n}\frac{\textrm{d}\Omega_{e}}{(2\pi)^{2}}\,m_{e}p_{e}\,\left(\frac{2\pi\,Z\alpha}{3}\right)^{2}\left(\frac{\mu}{m_{N}}\right)^{4}\,\left|\psi_{e,f}^{*}(0)\,\frac{R_{n0}^{e}(0)}{\sqrt{4\pi}}\,\int\textrm{d}r~{}4\pi
r^{4}\,\frac{R_{\text{b.s.}}^{\chi}(r)}{\sqrt{4\pi}}\,\frac{G(r)}{\sqrt{v_{\text{rel}}}}\right|^{2}\,,$
(S26)

where for $E_{B}=2.5\ \text{keV}$ in xenon the sum is over the $n=$ 3, 4, and
5 electrons, and the factor of two accounts for the pair of electrons in each
of these $s$ orbitals. Also note that
$p_{e}=\sqrt{2m_{e}(E_{B}-\omega_{e,i})}$ depends on $n$ through the orbital-
dependent electron binding energy $\omega_{e,i}$. Integrating over the
outgoing electron solid angle, the above expression simplifies to

$\displaystyle\sigma_{s}v_{\text{rel}}=\frac{4\pi\,(Z\alpha)^{2}}{9}\,\left(\frac{\mu}{m_{N}}\right)^{4}\left|\int\textrm{d}r~{}r^{4}\,R_{\text{b.s.}}^{\chi}(r)\,G(r)\right|^{2}~{}\sum_{n}\,2\,m_{e}p_{e}\,\left|\psi_{e,f}^{*}(0)\,R_{n0}(0)\right|^{2}\,.$
(S27)

Noting that the continuum wavefunctions $G$ and $\psi_{e,f}$ are
dimensionless, while the bound state wavefunctions $R_{\text{b.s.}}^{\chi}$
and $R_{n0}$ have dimension $[\text{length}]^{-3/2}$, we can use the
characteristic length scale of the $\chi$-$N$ bound state $1/\sqrt{2\mu
E_{B}}$ as well as the size of the atom $1/(Z\alpha m_{e})$ to construct the
following dimensionless form factors:

$\displaystyle F_{\chi,s}^{2}=(2\mu
E_{B})^{7/2}\,\left|\int\textrm{d}r~{}r^{4}\,R_{\text{b.s.}}^{\chi}(r)\,G(r)\right|^{2}~{}~{},~{}~{}F_{e}^{2}=\frac{2}{(Z\alpha
m_{e})^{3}}\sum_{n}\frac{p_{e}}{m_{e}}\,\left|\psi^{*}_{e,f}(0)\,R_{n0}(0)\right|^{2}~{}.$
(S28)

Our final result for the cross section is then

$\displaystyle\sigma_{s}v_{\text{rel}}$
$\displaystyle=\frac{4\pi}{9}\,\frac{(Z\alpha m_{e})^{5}}{(2\mu
E_{B})^{7/2}}\,\Big{(}\frac{\mu}{m_{N}}\Big{)}^{4}\,F_{\chi,s}^{2}\,F_{e}^{2}\,,$
(S29)

or numerically

$\displaystyle\sigma_{s}v_{\text{rel}}$
$\displaystyle\simeq$7\text{\times}{10}^{-34}\text{\,}{\mathrm{cm}}^{2}$\times\left(\frac{Z}{54}\right)^{5}\bigg{(}\frac{$122\text{\,}\mathrm{Ge\kern-0.50003ptV}$}{m_{N}}\bigg{)}^{4}\bigg{(}\frac{$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$}{E_{B}}\bigg{)}^{7/2}\bigg{(}\frac{\mu}{$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$}\bigg{)}^{1/2}\bigg{(}\frac{F_{\chi,s}^{2}}{49}\bigg{)}\bigg{(}\frac{F_{e}^{2}}{0.5}\bigg{)}\,.$
(S30)

Taking the ionization potentials for the $3s$, $4s$, and $5s$ xenon states to
be $1148.7\text{\,}\mathrm{e\kern-0.50003ptV}$,
$213.2\text{\,}\mathrm{e\kern-0.50003ptV}$, and
$23.3\text{\,}\mathrm{e\kern-0.50003ptV}$ Cardona and Ley (1978),
respectively, for $E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ we find
$F_{e}^{2}=0.5$ with a relative contribution of 81%, 17% and 2% from each
state. For $F_{\chi,s}^{2}$, we find only a weak dependence on $m_{\chi}$,
ranging from $F_{\chi,s}^{2}\approx 35$ for
$m_{\chi}=$1\text{\,}\mathrm{Ge\kern-0.50003ptV}$$ to $F_{\chi,s}^{2}\approx
62$ for $m_{\chi}\gg m_{\text{Xe}}$. A numerical fit to the capture cross
section in Xe with $E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ gives

$\displaystyle\sigma_{\text{Xe},s}v_{\text{rel}}\simeq$6\text{\times}{10}^{-34}\text{\,}{\mathrm{cm}}^{2}$\times\left(\frac{\mu}{$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$}\right)^{0.55}\,.$
(S31)

Fig. S1 shows $\sigma_{\text{Xe},s}v_{\text{rel}}$ as a function of
$m_{\chi}$. In addition, we also show _i)_ the capture cross section in
barium, assuming $\chi$ binds to Xe with
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$, which determines the
survival probability of an incoming $\chi$ travelling through the terrestrial
overburden of underground detectors, and _(ii)_ the capture cross section in
thallium, assuming $\chi$ binds to Tl with
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$, which is relevant for
signals in DAMA (see the main body for further discussion).

Figure S1: The capture cross sections for the $s$-to-$s$ process (solid) and
$p$-to-$s$ process (dashed), for capture of $\chi$ onto xenon (red), barium
(green), and thallium (blue), as a function of dark matter mass $m_{\chi}$.
Cross sections for Xe and Ba assume the existence of a $\chi-\text{Xe}$ bound
state with $E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$, while the
calculation for Tl assumes that $\chi$ binds to Tl with
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$.

## p-to-s Capture Cross Section

We now consider the case of $p$-to-$s$ capture through $\hat{V}_{eA}^{(1)}$.
This process results in a dipole transition of the final state electron, and
bears many similarities to the photoelectric cross section of the atom.

To begin, let us consider the matrix element of Eq. (S1), which in position
space involves $V_{eN}^{(1)}=\mathbf{r}_{N}^{i}\,\partial^{i}\,V_{eA}^{(0)}$.
In the language of operators, this corresponds to
$\hat{V}_{eN}^{(1)}=i\,\hat{\mathbf{r}}_{N}^{i}\,\,[\,\hat{\mathbf{p}}_{e}^{i},\hat{V}_{eA}^{(0)}\,]$.
The matrix element of $\hat{V}_{eN}^{(1)}$ between the initial state
$|i\rangle=|\chi_{i}\rangle\,|e_{i}\rangle$ and final state
$|f\rangle=|\chi_{f}\rangle\,|e_{f}\rangle$ is then

$\langle
f|\hat{V}_{eN}^{(1)}|i\rangle=i\,\langle\chi_{f}|\hat{\mathbf{r}}_{N}^{i}|\chi_{i}\rangle~{}\langle
e_{f}|\,[\,\hat{\mathbf{p}}_{e}^{i},\hat{V}_{eA}^{(0)}\,]\,|e_{i}\rangle=i\,\langle\chi_{f}|\hat{\mathbf{r}}_{N}^{i}|\chi_{i}\rangle~{}\langle
e_{f}|\,[\,\hat{\mathbf{p}}_{e}^{i},\hat{H}^{(0)}\,]\,|e_{i}\rangle=-iE_{B}\,\langle\chi_{f}|\hat{\mathbf{r}}_{N}^{i}|\chi_{i}\rangle~{}\langle
e_{f}|\hat{\mathbf{p}}_{e}^{i}|e_{i}\rangle~{},$ (S32)

where in the second and third equalities we have made use of the fact that the
only contribution to $\hat{H}_{0}$ that does not commute with
$\hat{\mathbf{p}}_{e}$ is $\hat{V}_{eA}$, and that the initial and final
electron states are both eigenstates of $\hat{H}_{0}$, with the difference in
energy between the two states being given by the $\chi$-$N$ binding energy,
$E_{B}$. The square of this matrix element can be simplified by averaging over
the direction of the electron momentum and using Eq. (S15), which allows us to
replace

$\big{|}\langle
f|\hat{V}_{eN}^{(1)}|i\rangle\big{|}^{2}\to\frac{1}{3}\,E_{B}^{2}\,\big{|}\langle\chi_{f}|\hat{\mathbf{r}}_{N}|\chi_{i}\rangle\big{|}^{2}~{}\big{|}\langle
e_{f}|\hat{\mathbf{p}}_{e}|e_{i}\rangle\big{|}^{2}~{}.$ (S33)

Notice that we have successfully factored the contributions into a DM-only
piece and an electron-only piece. The electron part of the matrix element is
the same one that appears in electromagnetic dipole transitions. In
particular, the photoelectric cross section in the long wavelength limit
(averaging over photon polarization) is $\sigma_{\text{pe}}=2\alpha
p_{e}\,\big{|}\langle
e_{f}|\hat{\mathbf{p}}_{e}|e_{i}\rangle\big{|}^{2}/(3m_{e}E_{B})$ Sakurai and
Napolitano (2020) for an incoming photon with energy $E_{B}$ and outgoing
electron with momentum $p_{e}$. This allows Eq. (S33) to be rewritten as

$\big{|}\langle
f|\hat{V}_{eN}^{(1)}|i\rangle\big{|}^{2}\to~{}\frac{m_{e}E_{B}^{3}}{2\alpha
p_{e}}\,\sigma_{\text{pe}}~{}\big{|}\langle\chi_{f}|\hat{\mathbf{r}}_{N}|\chi_{i}\rangle\big{|}^{2}~{}.$
(S34)

The capture cross section is determined from this matrix element as in Eq.
(S1). Rewriting the DM matrix element in position basis and using
$r_{N}=(\mu/m_{N})\,r$, we have

$\displaystyle\textrm{d}\sigma_{p}=\frac{\textrm{d}\Omega_{e}}{(2\pi)^{2}}\,\frac{m_{e}^{2}E_{B}^{3}}{2\alpha}\,\left(\frac{\mu}{m_{N}}\right)^{2}\,\sigma_{\text{pe}}\,~{}\left|\int\textrm{d}^{3}\mathbf{r}~{}\psi_{\chi,f}^{*}(\mathbf{r})\,\psi_{\chi,i}(\mathbf{r})~{}\mathbf{r}\right|^{2}\,.$
(S35)

The DM wavefunctions can be derived in a similar manner to the the previous
section. The final bound state wavefunction $\psi_{\chi,f}$ is identically
given by the $s$-wave state in Eq. (S22). Since $\psi_{\chi,f}$ is spherically
symmetric, we see that the integral in Eq. (S35) vanishes if $\psi_{\chi,i}$
is also $s$-wave. Hence, we take the incoming wavefunction $\psi_{\chi,i}$ as
the $\ell=1$ term in the continuum wavefunction expansion in Eq. (S13), which
gives

$\displaystyle\psi_{\chi,i}=\frac{3i}{2k\sqrt{v_{\text{rel}}}}\,(\hat{\mathbf{k}}\cdot\hat{\mathbf{r}})\,R_{k1}(r)\,,$
(S36)

where $R_{k1}(r)$ is the solution to Eq. (S12) with $l=1$, $V=V_{\chi N}$, and
$E\simeq 0$. Far away from the nucleus, $R_{k1}$ should tend toward the free
$l=1$ solution $R_{k1}(r)\to 2kj_{1}(kr)\simeq 2k^{2}r/3$ for $k\ll 1/r$. As
in the last section, to incorporate the effect of the potential $V_{\chi N}$
on $\psi_{\chi,i}(r)$ at smaller radii, we parametrize the general solution as

$\displaystyle\psi_{\chi,i}=i(\mathbf{k}\cdot\mathbf{r})\,F(r)/\sqrt{v_{\text{rel}}}~{}~{},~{}~{}\lim_{r\to\infty}F(r)=1~{},$
(S37)

for some dimensionless function $F(r)$ that we compute numerically by solving
Eq. (S12) for the $\chi-N$ system.

Having determined the DM wavefunctions, Eq. (S35) reduces to

$\displaystyle\sigma_{p}v_{\text{rel}}=\frac{2m_{e}^{2}\,E_{B}^{3}\,\mu^{3}\,T}{3\alpha\,m_{N}^{2}}\,\sigma_{\text{pe}}\,\left|\int\textrm{d}r~{}R_{\text{b.s.}}^{\chi}(r)\,F(r)\,r^{4}\right|^{2}\,,$
(S38)

where we integrated over the outgoing electron solid angle and used Eq. (S15)
to replace $(\mathbf{k}\cdot\mathbf{r})\,\mathbf{r}\to(r^{2}/3)\,\mathbf{k}$
in the integral over $r$. Finally, we replaced $k^{2}$ with its thermally-
averaged value $k^{2}\to 3\mu T$, where $T\simeq 300\ \text{K}$ is the
temperature of the thermalized $\chi$ particle; as a result, the cross section
is suppressed by the small thermal velocity. As in the last section, we can
construct the following dimensionless quantity

$\displaystyle F_{\chi,p}^{2}=(2\mu
E_{B})^{7/2}\left|\int\textrm{d}r~{}R_{\text{b.s.}}^{\chi}(r)\,F(r)\,r^{4}\right|^{2}\,,$
(S39)

such that the $p$-to-$s$ capture cross section is given by

$\displaystyle\sigma_{p}v_{\text{rel}}=\frac{m_{e}^{2}\,T\,F_{\chi,p}^{2}\,\sigma_{\text{pe}}}{12\alpha\,m_{N}^{2}\,(2\mu
E_{B})^{1/2}}~{}.$ (S40)

For parameters that are representative of $E_{B}=2.5\ \text{keV}$ in xenon, we
find

$\displaystyle\sigma_{p}v_{\text{rel}}\simeq${10}^{-35}\text{\,}{\mathrm{cm}}^{2}$\times\bigg{(}\frac{T}{$300\text{\,}\mathrm{K}$}\bigg{)}\bigg{(}\frac{$122\text{\,}\mathrm{Ge\kern-0.50003ptV}$}{m_{N}}\bigg{)}^{2}\bigg{(}\frac{$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$}{E_{B}}\bigg{)}^{1/2}\bigg{(}\frac{$100\text{\,}\mathrm{Ge\kern-0.50003ptV}$}{\mu}\bigg{)}^{1/2}\bigg{(}\frac{F_{\chi,p}^{2}}{176}\bigg{)}\bigg{(}\frac{\sigma_{pe}}{$2.5\text{\times}{10}^{-19}\text{\,}{\mathrm{cm}}^{2}$}\bigg{)}\,.$
(S41)

The photoelectric cross section $\sigma_{\text{pe}}$ for various elements,
including xenon, with a photon energy of
$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$ can be interpolated from data in
Refs. West and Morton (1978); Veigele (1973). Fixing
$E_{B}=$2.5\text{\,}\mathrm{ke\kern-0.50003ptV}$$ in xenon, we find only a
weak dependence of $F_{\chi,p}^{2}$ on $m_{\chi}$, ranging from
$F_{\chi,p}^{2}\approx 81$ for
$m_{\chi}=$1\text{\,}\mathrm{Ge\kern-0.50003ptV}$$ to $F_{\chi,p}^{2}\approx
287$ for $m_{\chi}\gg m_{N}$.

The dashed lines in Fig. S1 show $\sigma_{p}v_{\text{rel}}$ for various
elements as a function of $m_{\chi}$. For xenon, we find that the $p$-to-$s$
rate is smaller than the $s$-to-$s$ rate, with more than an order of magnitude
suppression once $m_{\chi}\gtrsim$10\text{\,}\mathrm{Ge\kern-0.50003ptV}$$;
for simplicity, we can therefore neglect the contribution of $\sigma_{p}$ to
the full capture cross section in xenon. For capture in barium and thallium,
however, the $p$-to-$s$ cross section dominates over the $s$-to-$s$ rate for
$m_{\chi}\lesssim$5\text{\,}\mathrm{Ge\kern-0.50003ptV}$$.