#  Singularities in the Gravitational Capture of Dark Matter through Long-
Range Interactions

Cristian Gaidau  and Jessie Shelton

###### Abstract

We re-examine the gravitational capture of dark matter (DM) through long-range
interactions. We demonstrate that neglecting the thermal motion of target
particles, which is often a good approximation for short-range capture,
results in parametrically inaccurate results for long-range capture. When the
particle mediating the scattering process has a mass that is small in
comparison to the momentum transfer in scattering events, correctly
incorporating the thermal motion of target particles results in a quadratic,
rather than logarithmic, sensitivity to the mediator mass, which substantially
enhances the capture rate. We quantitatively assess the impact of this finite
temperature effect on the captured DM population in the Sun as a function of
mediator mass. We find that capture of DM through light dark photons, as in
e.g. mirror DM, can be powerfully enhanced, with self-capture attaining a
geometric limit over much of parameter space. For visibly-decaying dark
photons, thermal corrections are not large in the Sun, but may be important in
understanding long-range DM capture in more massive bodies such as Population
III stars. We additionally provide the first calculation of the long-range DM
self-evaporation rate.

##  1 Introduction

The gravitational capture of dark matter (DM) in massive bodies such as stars
and planets is a powerful probe of DM’s non-gravitational interactions with
the SM. Dark matter that undergoes scattering with the material in a massive
body can lose enough energy in the scattering event to become gravitationally
bound. Through such processes the massive body can build up a population of
captured dark matter, which can in turn give rise to a variety of
observational signatures depending on the detailed properties of DM. If the
captured DM can self-annihilate, then its annihilation products may be
detected either directly, if they escape the body [1, 2, 3], or indirectly,
through the energy they deposit into the capturing body [4, 5, 6, 7, 8, 9].
Even in the absence of annihilations, the captured DM population may still be
detectable through its gravitational effect on the capturing body [10, 11, 12,
13, 14, 15, 16, 17].

In computing the DM capture rate for non-degenerate systems, e.g., main
sequence stars such as our Sun, planets such as our Earth, and nuclei in white
dwarfs, it is common to neglect the temperature of the target particles. For
short-range interactions, incorporating the finite temperature of the targets
gives corrections of less than 10% to the zero-temperature nuclear capture
rate in the Sun for DM heavier than the evaporation mass [18, 19], and thus in
this case the zero-temperature approximation for DM capture is usually well
justified. However, we show here that for long-range interactions, the thermal
motion of target nuclei has a far more dramatic effect on the capture rate.
The zero-temperature capture rate through long-range interactions is
logarithmically divergent [20], and is thus only logarithmically sensitive to
the (model-dependent) regulator, which could be a screening scale [20] or an
intrinsic mediator mass [21]. We find that retaining the finite temperature of
the target nuclei instead results in a quadratically divergent capture rate in
the long-range regime, resulting in both a substantially enhanced capture rate
and increased sensitivity to the physics responsible for regulating the
capture rate in the infrared. Similar results hold for the self-capture of DM
through long-range interactions.

Capture via long-range interactions is particularly relevant for DM that
interacts with itself and/or the SM via a new dark force. DM self-interactions
in such models can lead to novel cosmological and/or astrophysical signatures
[22, 23, 24, 25, 26, 27]. These models generally feature secluded DM
annihilations into mediators that may be long-lived, thus providing novel
cosmic ray signatures [2, 3, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39].
Thermal target motion is known to be important in cases where short-range
scattering on electrons dominates gravitational capture in the Sun, due to the
much larger thermal velocities of electrons compared to nuclei [40]. The
effect we point out in this paper is qualitatively distinct, and applies for
long-range interactions between DM and any non-degenerate population of
targets, including electrons. For concreteness, we illustrate the general
point with numerical results for the familiar case of nuclear capture in the
Sun, which offers the most immediate detection prospects, but discuss how our
results apply to other systems.

In Secs. 2 and 3 we discuss long-range nuclear capture and long-range DM self-
capture, respectively. Both sections start with a review of the zero-
temperature calculation, followed by the finite-temperature calculation
incorporating the thermal motion of target particles. We quantify the impact
of the finite target temperature and discuss its observational and
phenomenological implications in Sec. 4. We conclude in Sec. 5, and provide
additional technical details for the computation of nuclear and self-capture
rates in Appendices A and B, respectively. Finally, we discuss DM self-
evaporation via long-range forces in Appendix C.

##  2 Nuclear capture with thermal targets

We begin with the calculation of the rate of DM gravitational capture in the
Sun by scattering against nuclei, given by the nuclear capture coefficient
$C_{\mathrm{c}}$. This rate depends on three ingredients: a model of the DM
halo in the vicinity of the Sun, a model of the solar interior, and a model
for the DM-nuclear interaction cross-section.

For DM particles in the halo we assume a Maxwell distribution of velocities,

$\displaystyle
f(u)d^{3}u=4\pi\left(\frac{3}{2\pi\bar{v}^{2}}\right)^{3/2}\exp{\left(-\frac{3u^{2}}{2\bar{v}^{2}}\right)}u^{2}du,$
(2.1)

where $u$ is the DM speed far away from the Sun and $\bar{v}^{2}$ is the DM
rms speed. To find the distribution of initial DM speeds as seen by the Sun,
we perform a Galilean transformation to the Sun’s rest frame and average over
the solid angle:

$\displaystyle
f_{\eta}(u)=\left(\frac{3}{2\pi\bar{v}^{2}}\right)^{3/2}\exp{\left(-\frac{3(u^{2}+\tilde{v}^{2})}{2\bar{v}^{2}}\right)}\frac{\sinh{(\frac{3}{2\bar{v}^{2}}2u\tilde{v})}}{\frac{3}{2\bar{v}^{2}}2u\tilde{v}},$
(2.2)

where $\tilde{v}$ is the Galilean boost from the Galactic rest frame to the
Sun rest frame. Assuming a virialized DM halo,
$\bar{v}^{2}=3/2\,v^{2}_{\text{LSR}}$, where $v_{\text{LSR}}$ is the Local
Standard of Rest at the position of the Sun. Refs. [41, 42, 43] estimate
$v_{\text{LSR}}=235\text{ km}/\text{s}$, so that $\bar{v}=288\text{
km}/\text{s}$. Given the Sun’s peculiar velocity
${\bf{v}}_{\text{pec}}=(11,12,7)\text{ km}/\text{s}$ [44], it moves at
$\tilde{v}=247\text{ km}/\text{s}$ with respect to the DM halo. Meanwhile, the
number density of dark matter in the solar neighborhood is
$n_{\text{DM}}=\rho_{\text{DM}}/M_{\text{DM}}$, where we take
$\rho_{\text{DM}}=0.4\,\mathrm{GeV}/\mathrm{cm}^{3}$ [45].

We model the density and composition of the Sun using the BS05(AGS, OP) solar
model [46]. For DM-nuclear interactions, we consider for concreteness a
minimal reference model that consists of a fermionic DM particle $\chi$ with
mass $M_{\text{DM}}$, which couples to a massive $U(1)$ gauge boson $A_{D}$
with mass $M_{\gamma_{D}}$. We take $A_{D}$ to couple to the SM via kinetic
mixing with the hypercharge gauge boson [47]. The diagonalization of the gauge
field kinetic terms results in the mostly-$A_{D}$ mass eigenstate picking up
interactions with the SM via small mixings with hypercharge. In this paper we
will focus on a light $A_{D}$ boson, such that its couplings to the SM are
photon-like. The Lagrangian in the mass basis then reads

$\displaystyle\mathcal{L}=-\frac{1}{4}F_{D\mu\nu}F^{\mu\nu}_{D}+\frac{1}{2}M^{2}_{\gamma_{D}}A^{2}_{D}+\bar{\chi}(i\not{D}-M_{\text{DM}})\chi-\sum_{f}q_{f}e(A_{\mu}+\epsilon
A_{D\mu})\bar{f}\gamma^{\mu}f,$ (2.3)

where $\not{D}=\not{\partial}+ig_{D}\not{A}_{D}$, $\epsilon$ controls the
kinetic mixing and $f$ is a SM fermion with electric charge $q_{f}$. For
previous studies of gravitational capture in this model, see Refs. [32, 33,
39].

###  2.1 Zero-temperature calculation

We begin with a review of the zero-temperature calculation. We construct the
capture coefficient following the procedure outlined in Gould’s work [48, 18,
49], where the differential capture rate is integrated over spherical shells
of the capturing body,

$\displaystyle C_{\text{c}}=\sum_{i}C_{\text{c},i}=\sum_{i}\int_{0}^{R}dr4\pi
r^{2}\frac{dC_{\text{c},i}}{dV},$ (2.4)

where $R$ is the radius of the Sun, and the sum runs over the contributions of
different nuclear species. The capture rate per unit volume is

$\displaystyle\frac{dC_{\text{c},i}}{dV}=n_{i}(r)n_{\text{DM}}\int
d^{3}\vec{u}\frac{w}{u}f_{\eta}(u)|\vec{w}|\sigma_{c,i},$ (2.5)

where $n_{i}$ is the number density of nuclear species $i$, $n_{\text{DM}}$ is
the local DM halo number density, $\sigma_{c,i}$ is the DM-nucleus cross
section for scatterings that lead to capture, and
$w^{2}=u^{2}+v^{2}_{\text{esc}}(r)$ is the DM speed before the collision,
reflecting the increase in kinetic energy as the incoming DM falls into the
Sun’s potential well.

We will consider the scattering process in the center of mass (CM) frame. Let
$\vec{v}^{\prime}_{1}$ and $\vec{v}_{n}^{CM}$ be the CM velocities of the DM
particle and nucleus before the collision. The Galilean transformation to the
CM frame reads:

$\displaystyle\vec{w}=\vec{v}^{\prime}_{1}+\vec{v}_{CM},\phantom{space}\vec{v}_{n}=0=\vec{v}^{CM}_{n}+\vec{v}_{CM},$
(2.6)

such that $m_{i}\vec{v}^{CM}_{n}+M_{\text{DM}}\vec{v}^{\prime}_{1}=0$. We note
that

$\displaystyle\vec{v}_{CM}=\frac{M_{\text{DM}}}{m_{i}}\vec{v}^{\prime}_{1}\quad\text{and}\quad\vec{w}=\left(1+\frac{M_{\text{DM}}}{m_{i}}\right)\vec{v}^{\prime}_{1}.$
(2.7)

The vectors $\vec{w},\vec{v}_{CM}$ and $\vec{v}^{\prime}_{1}$ are all
collinear, which greatly simplifies the angular structure of the integral in
Eq. 2.5. We now switch the integration variable to $d^{3}\vec{w}$, giving

$\displaystyle\frac{dC_{\text{c},i}}{dV}=n_{i}(r)n_{\text{DM}}\int
d^{3}\vec{w}\Theta(|\vec{w}|-v_{\text{esc}}(r))f_{\eta}(\sqrt{w^{2}-v^{2}_{\text{esc}}(r)})|\vec{w}|\sigma_{c,i}.$
(2.8)

We turn to the construction of the capture cross-section $\sigma_{c,i}$. DM
scattering off nuclei proceeds via the $t$-channel exchange of $A_{D}$. In the
non-relativistic limit, the matrix element-squared for this process reads

$\displaystyle|\mathcal{M}|^{2}=\frac{4g_{D}^{2}\epsilon^{2}e^{2}}{v_{1}^{{}^{\prime}4}}\frac{m_{n}^{2}}{M_{\text{DM}}^{2}}\frac{1}{(1-\cos\theta+\mu)^{2}},$
(2.9)

where $\theta$ is the scattering angle in the CM frame, $m_{n}$ is the mass of
the nucleon and $\mu$ regulates the forward singularity. For a massive dark
gauge boson of mass $M_{\gamma_{D}}$,

$\displaystyle\mu\equiv\frac{M^{2}_{\gamma_{D}}}{2M^{2}_{\text{DM}}v_{1}^{{}^{\prime}2}}.$
(2.10)

In the case of coherent scattering, the nuclear cross-section can be simply
written in terms of the nucleon cross-section

$\displaystyle\sigma_{c,i}=\left(\frac{\bar{m}_{i}}{\bar{m}}\right)^{2}Z_{i}^{2}\sigma_{c},$
(2.11)

where $Z_{i}$ is the number of protons in nucleus $i$, $\bar{m}\equiv
M_{\text{DM}}m_{n}/(M_{\text{DM}}+m_{n})$ is the reduced DM-nucleon mass and
$\bar{m}_{i}\equiv M_{\text{DM}}m_{i}/(M_{\text{DM}}+m_{i})$ is the reduced
DM-nucleus mass. Using the matrix element in Eq. 2.9 above, we now calculate
$\sigma_{c}$:

$\displaystyle\sigma_{c}=\frac{1}{16m_{n}^{2}M_{\text{DM}}^{2}}\frac{\bar{m}^{2}}{(2\pi)^{2}}\int
d\Omega_{f,\chi}|\mathcal{M}|^{2}\Theta\left(v_{\text{esc}}(r)-|\vec{v}^{lab}_{\text{DM},f}|\right),$
(2.12)

where $d\Omega_{f,\chi}=\sin\theta d\theta d\phi_{\chi}$ is the solid angle
specifying the direction of the CM-frame velocity of the DM particle after the
collision. Here the step function enforces the kinematic requirement that the
outgoing DM particle is gravitationally bound. Note that
$\vec{v}_{\text{DM},f}^{lab}=\vec{v}^{\prime}_{1,f}+\vec{v}_{CM}$ and
$|\vec{v}^{\prime}_{1,f}|=|\vec{v}_{1}^{\prime}|$. The step function thus
imposes an upper bound on the $\cos\theta$ integral:

$\displaystyle\cos\theta\leq\cos\theta_{+}\equiv\frac{m_{i}}{2M_{\text{DM}}}\left(\frac{v^{2}_{\text{esc}}(r)}{v^{{}^{\prime}2}_{1}}-1-\left(\frac{M_{\text{DM}}}{m_{i}}\right)^{2}\right)\leq
1,$ (2.13)

where the second inequality is a consequence of the fact that $w\geq
v_{\text{esc}}(r)$. For sufficiently energetic incoming DM particles,
$\cos\theta_{+}<-1$ and gravitational capture cannot occur for any scattering
angle $\theta$. Therefore, requiring $\cos\theta_{+}\geq-1$ imposes an upper
bound on $w$ for events leading to capture:

$\displaystyle w\leq
v_{\text{esc}}(r)\frac{m_{i}+M_{\text{DM}}}{|m_{i}-M_{\text{DM}}|}.$ (2.14)

The capture cross-section can be evaluated in closed form,

$\displaystyle\sigma_{c,i}=2\pi
Z_{i}^{2}\frac{\alpha\alpha_{D}\epsilon^{2}}{\bar{m}_{i}^{2}w^{4}}\frac{1+\cos\theta_{+}}{(2+\mu)(1-\cos\theta_{+}+\mu)},$
(2.15)

where $\alpha_{D}=g_{D}^{2}/(4\pi)$ and $\alpha=e^{2}/(4\pi)$. The kinematic
integral now reads:

$\displaystyle\frac{dC_{\text{c},i}}{dV}=n_{i}(r)n_{\text{DM}}4\pi\int_{v_{\text{esc}}(r)}^{v_{\text{esc}}(r)\frac{m_{i}+M_{\text{DM}}}{|m_{i}-M_{\text{DM}}|}}dwf_{\eta}(\sqrt{w^{2}-v^{2}_{\text{esc}}(r)})w^{3}\sigma_{c,i}.$
(2.16)

We note that $\sigma_{c,i}$ is finite in the limit $\mu\rightarrow 0$ except
at $w=v_{\text{esc}}(r)$, where $\cos\theta_{+}=1$. Upon integrating over the
incoming DM speed $w$ in Eq. 2.16, this singularity in $\sigma_{c,i}$ results
in $C_{\text{c}}$ having a logarithmic sensitivity on the regulator $\mu$ [20,
50].

Numerical calculations of $C_{\text{c}}$ with $T=0$ are shown as dashed lines
in Fig. 1, where we examine the dependence of the capture coefficient on
$M_{\gamma_{D}}$ at fixed $M_{\text{DM}}$ and vice-versa. Here we consider
fixed couplings $\alpha_{D}=10^{-3}$ and $\epsilon=10^{-3}$.

The left panel of Fig. 1 illustrates two distinct regimes for the capture rate
dependence on the mediator mass $M_{\gamma_{D}}$. When the dark photon mass is
sufficiently large compared to the typical momentum transfer in a capture
event, i.e., $\mu\gg 1$, the interaction is effectively pointlike, which in
this model yields a constant, velocity-independent cross section proportional
to $M_{\gamma_{D}}^{-4}$. In the opposite limit, when $M_{\gamma_{D}}$ is much
smaller than the typical momentum transfer and $\mu\ll 1$, the zero-
temperature capture rate depends logarithmically on $M_{\gamma_{D}}$.

###  2.2 Nuclear Capture at Finite Temperature

We extend our previous treatment of nuclear capture by “restoring” the thermal
motion of nuclear targets, which can be described by a Maxwell distribution

$\displaystyle f_{i}(\vec{v}_{n})=\left(\frac{m_{i}}{2\pi
T_{\text{Sun}}}\right)^{3/2}\exp{\left(-\frac{m_{i}v_{n}^{2}}{2T_{\text{Sun}}}\right)}$
(2.17)

in the Sun’s rest frame, where $m_{i}$ is the mass of the nucleus, $v_{n}$ is
the velocity of the nucleus and $T_{\text{Sun}}$ is the local temperature in
the Sun. We take $T_{\text{Sun}}$ to be the (constant) core temperature for
simplicity.

In addition to the integral over the distribution of incoming DM speeds, the
capture rate per unit shell now includes a convolution over the thermal
distribution of the target particle velocity,

$\displaystyle\frac{dC_{\text{c},i}}{dV}=n_{i}(r)n_{\text{DM}}\int
d^{3}\vec{u}\frac{w}{u}f_{\eta}(u)\int
d^{3}\vec{v}_{n}f_{i}(\vec{v}_{n})|\vec{w}-\vec{v}_{n}|\sigma_{c,i}.$ (2.18)

The scattering process itself remains simplest to describe in the CM frame;
however, the kinematic conditions for capture depend on the outgoing
velocities in the lab frame. This complicates the evaluation of the final
state phase space integrals, but as we will see being careful about the
capture conditions is important for obtaining the correct parametric
dependence of the capture rate. Transforming the velocity integrals $\int
d^{3}\vec{u}\int d^{3}\vec{v}_{n}$ to the CM frame, we have for the incoming
DM

$\displaystyle\int d^{3}\vec{u}\frac{w}{u}f_{\eta}(u)=\int
d^{3}\vec{w}\Theta(|\vec{w}|-v_{\text{esc}}(r))f_{\eta}(\sqrt{w^{2}-v^{2}_{\text{esc}}(r)}).$
(2.19)

Let $\vec{v}^{\prime}_{1}$ and $\vec{v}_{n}^{CM}$ be the CM velocities of the
DM particle and nucleus before the collision. The Galilean transformation
reads

$\displaystyle\vec{w}=\vec{v}^{\prime}_{1}+\vec{v}_{CM},\phantom{space}\vec{v}_{n}=\vec{v}_{n}^{CM}+\vec{v}_{CM},$
(2.20)

where

$\displaystyle\vec{v}_{CM}=\frac{m_{i}}{m_{i}+M_{\text{DM}}}\vec{v}_{n}+\frac{M_{\text{DM}}}{m_{i}+M_{\text{DM}}}\vec{w}.$
(2.21)

Next we switch variables to $(\vec{v}_{CM},\vec{v}^{\prime}_{1})$, where the
Jacobian for this transformation is
$J=\left(1+M_{\text{DM}}/m_{i}\right)^{3}$, giving

$\displaystyle\frac{dC_{\text{c},i}}{dV}$
$\displaystyle=n_{i}(r)n_{\text{DM}}\int d^{3}\vec{v}_{CM}\int
d^{3}\vec{v}^{\prime}_{1}J\Theta\left(|\vec{v}_{CM}+\vec{v}^{\prime}_{1}|-v_{\text{esc}}(r)\right)f_{\eta}\left(\sqrt{(\vec{v}_{CM}+\vec{v}_{1}^{\prime})^{2}-v^{2}_{\text{esc}}(r)}\right)\times$
$\displaystyle\times
f_{i}\left(\vec{v}_{CM}+\vec{v}_{n}^{CM}\right)\left(1+\frac{M_{\text{DM}}}{m_{i}}\right)|\vec{v^{\prime}}_{1}|\sigma_{c,i}.$
(2.22)

Without loss of generality, we orient the $z$-axis of our frame along
$\vec{v}^{\prime}_{1}$, so that $d^{3}\vec{v}^{\prime}_{1}=4\pi
v^{{}^{\prime}2}_{1}dv^{\prime}_{1}$. We let $\alpha$ and $\phi_{CM}$ specify
the polar and azimuthal angles of $\vec{v}_{CM}$ with respect to the $z$-axis.
Next we introduce the following definition to ease notation:

$\displaystyle\beta(r)\equiv\frac{v^{2}_{\text{esc}}(r)-v_{CM}^{2}-v_{1}^{{}^{\prime}2}}{2v_{CM}v^{\prime}_{1}}.$
(2.23)

Using the definitions above we write the capture rate per unit volume as an
integral over four kinematic variables,
$(v_{CM},v^{\prime}_{1},c\alpha,\phi_{CM})$, where for brevity we write
$c\alpha\equiv\cos\alpha$ and $s\alpha\equiv\sin\alpha$,

$\displaystyle\frac{dC_{\text{c},i}}{dV}$
$\displaystyle=n_{i}(r)n_{\text{DM}}\int v^{2}_{CM}dv_{CM}dc\alpha
d\phi_{CM}\int
v^{{}^{\prime}2}_{1}dv^{\prime}_{1}4\pi\Theta(c\alpha-\beta(r))\times$
$\displaystyle\times
f_{\eta}\left(\sqrt{(\vec{v}_{CM}+\vec{v}_{1}^{\prime})^{2}-v^{2}_{\text{esc}}(r)}\right)f_{i}\left(\vec{v}_{CM}+\vec{v}_{n}^{CM}\right)\left(1+\frac{M_{\text{DM}}}{m_{i}}\right)^{4}v^{\prime}_{1}\sigma_{c,i}.$
(2.24)

The two Maxwell distributions in the CM frame read

$\displaystyle
f_{i}\left(\vec{v}_{CM}+\vec{v}_{n}^{CM}\right)=\left(\frac{m_{i}}{2\pi
T_{\text{Sun}}}\right)^{3/2}$
$\displaystyle\exp{\left(-\frac{m_{i}}{2T_{\text{Sun}}}\left(v^{2}_{CM}+\frac{M_{\text{DM}}^{2}}{m_{i}^{2}}v_{1}^{{}^{\prime}2}-2\frac{M_{\text{DM}}}{m_{i}}v_{CM}v_{1}^{\prime}c\alpha\right)\right)}$
(2.25)

for the nuclei, and

$\displaystyle
f_{\eta}\left(\sqrt{(\vec{v}_{CM}+\vec{v}_{1}^{\prime})^{2}-v^{2}_{\text{esc}}(r)}\right)$
$\displaystyle=\frac{1}{\pi^{3/2}}\frac{e^{-\eta^{2}}}{2\eta}\frac{3}{2\bar{v}^{2}}\frac{\exp{\left(-\frac{3}{2\bar{v}^{2}}2v_{CM}v^{\prime}_{1}(c\alpha-\beta(r))\right)}}{\sqrt{2v_{CM}v^{\prime}_{1}(c\alpha-\beta(r))}}\times$
$\displaystyle\times\sinh{\left(\sqrt{\frac{3}{2\bar{v}^{2}}}2\eta\sqrt{2v_{CM}v^{\prime}_{1}(c\alpha-\beta(r))}\right)}$
(2.26)

for the DM particles, where $\eta^{2}\equiv 3\tilde{v}^{2}/(2\bar{v}^{2})$.

Next we consider the calculation of the capture cross-section $\sigma_{c,i}$.
Similar to the zero-temperature calculation we will consider coherent nuclear
scattering, such that the nuclear cross-section is given in terms of the
nucleon cross-section, Eq. 2.11, with the matrix element squared again given
by Eq. 2.9. The major difference in the finite temperature calculation is that
the velocities $(\vec{v}_{CM},\vec{v}_{1}^{\prime},\vec{v}_{1,f}^{\prime})$
are no longer coplanar in general. The argument of the kinematic step function
that enforces the requirement that the final state DM particle is
gravitationally bound is therefore significantly more involved. After some
algebraic manipulations of Eq. 2.12, the capture cross-section can be written

$\displaystyle\sigma_{c}=\frac{\bar{m}^{2}}{M_{\text{DM}}^{4}}\frac{\alpha\alpha_{D}\epsilon^{2}}{v_{1}^{{}^{\prime}4}}\int
d\cos\theta
d\phi_{\chi}\frac{1}{(1-\cos\theta+\mu)^{2}}\Theta(\beta(r)-\cos\gamma),$
(2.27)

where $\gamma$ is the angle between $\vec{v}_{CM}$ and the final state DM
velocity in the CM frame, $\vec{v}_{1,f}^{\prime}$:

$\displaystyle\cos\gamma=\sin\alpha\sin\theta\cos(\phi_{\chi}-\phi_{CM})+\cos\alpha\cos\theta.$
(2.28)

Here we can see that the singularity structure is (as expected for cross-
sections with forward singularities) dependent on the integral over the polar
scattering angle, while the criteria for capture make the integration in Eq.
2.27 over the azimuthal angle $\phi_{\chi}$ nontrivial. As forward scattering
cannot be neglected for long-range capture, the simpler average treatment of
[51] is insufficient to capture the parametric dependence of the capture
cross-section in this regime.

Before we put all the calculations together, let us introduce one more
definition that collects together three of the four angular integrals:

$\displaystyle
J_{c}\equiv\frac{1}{4\pi}\int_{0}^{2\pi}d\phi_{CM}\int_{-1}^{1}d\cos\theta\int_{0}^{2\pi}d\phi_{\chi}\frac{1}{(1-\cos\theta+\mu)^{2}}\Theta(\beta(r)-\cos\gamma).$
(2.29)

The step functions in Eqs. 2.24 and 2.27 require
$-1\leq\cos\gamma\leq\beta(r)\leq c\alpha\leq 1$. We also note the reduced
functional dependence
$J_{c}(M_{\text{DM}},M_{\gamma_{D}},r,v_{CM},v^{\prime}_{1},c\alpha)=J_{c}(\beta(r),c\alpha,\mu)$.
We exploit this property to precompute $J_{c}$ numerically, before computing
the remaining kinematic integrals over $v_{CM},v^{\prime}_{1}$ and $c\alpha$.
We delegate the details of this calculation to Appendix A. Putting everything
together, the capture rate per unit volume reads

$\displaystyle\frac{dC_{\text{c},i}}{dV}$
$\displaystyle=n_{i}(r)n_{\text{DM}}(4\pi)^{2}Z_{i}^{2}\frac{\alpha\alpha_{D}\epsilon^{2}}{\bar{m}^{2}_{i}}\int_{0}^{\infty}dv_{CM}\int_{0}^{\infty}dv^{\prime}_{1}\int
dc\alpha\Theta(c\alpha-\beta(r))\Theta(1-\beta(r))\times$
$\displaystyle\times\Theta(\beta(r)+1)\frac{v^{2}_{CM}}{v^{\prime}_{1}}f_{i}\left(v_{CM},v^{\prime}_{1},c\alpha\right)f_{\eta}\left(v_{CM},v^{\prime}_{1},c\alpha\right)J_{c}(\beta(r),c\alpha,\mu).$
(2.30)

The kinematic constraint $-1\leq\beta(r)\leq 1$ translates into bounds on the
speeds $v_{CM}$ and $v_{1}^{\prime}$:

$\displaystyle v_{\text{esc}}(r)\leq v_{CM}+v_{1}^{\prime}\quad\text{and}\quad
v_{\text{esc}}(r)\geq|v_{CM}-v_{1}^{\prime}|.$ (2.31)

Thus, we finally obtain the following capture rate per volume shell:

$\displaystyle\frac{dC_{\text{c},i}}{dV}$
$\displaystyle=n_{i}(r)n_{\text{DM}}(4\pi)^{2}Z_{i}^{2}\frac{\alpha\alpha_{D}\epsilon^{2}}{\bar{m}^{2}_{i}}\int_{0}^{\infty}dv_{CM}\int_{|v_{\text{esc}}(r)-v_{CM}|}^{v_{\text{esc}}(r)+v_{CM}}dv^{\prime}_{1}\times$
$\displaystyle\times\int
dc\alpha\Theta(c\alpha-\beta(r))\frac{v^{2}_{CM}}{v^{\prime}_{1}}f_{i}\left(v_{CM},v^{\prime}_{1},c\alpha\right)f_{\eta}\left(v_{CM},v^{\prime}_{1},c\alpha\right)J_{c}(\beta(r),c\alpha,\mu).$
(2.32)

Together with the precomputed $J_{c}$, we numerically calculate the capture
rate

$\displaystyle C_{\text{c}}=\sum_{i}C_{\text{c},i}=\sum_{i}\int_{0}^{R}dr4\pi
r^{2}\frac{dC_{\text{c},i}}{dV}.$ (2.33)

The resulting capture coefficient is shown in Fig. 1 as a function of
$M_{\text{DM}}$ and $M_{\gamma_{D}}$.

Figure 1: Left: Solar capture rate as a function of $M_{\gamma_{D}}$ with
(solid) and without (dashed) nuclear thermal motion. The different colors
correspond to different DM masses. Dotted lines are included to highlight the
parametric dependence on $M_{\gamma_{D}}$ in the long and short-range regimes.
Right: Solar capture rate as a function of $M_{\text{DM}}$ with (solid) and
without (dashed) nuclear thermal motion for fixed $M_{\gamma_{D}}=10^{-6}$
GeV. We fix $\alpha_{D}=10^{-3}$, $\epsilon=10^{-3}$ and a uniform temperature
$T_{\text{Sun}}=1.57\times 10^{7}$ K.

We again note the presence of two qualitatively distinct scattering regimes,
corresponding to different scalings of $C_{\text{c}}$ with $M_{\gamma_{D}}$.
For mediator masses large compared to typical momentum transfers, we again
obtain the $M_{\gamma_{D}}^{-4}$ scaling indicative of an effectively
pointlike interaction, similar to the $T=0$ calculation. Numerically, we find
that for a contact interaction the full thermal calculation gives a capture
coefficient which is $\mathcal{O}(10\%)$ smaller than the $T=0$ calculation,
in agreement with previous studies [19]. This small decrease in the capture
coefficient comes from the small increase in the typical CM energy of the
collision, which reduces the volume of phase space that meets the kinematic
criteria for capture.

Qualitatively new behavior occurs for small mediator masses, where we find
that $C_{\text{c}}$ is proportional to $M_{\gamma_{D}}^{-2}$. This parametric
dependence on the mediator mass is much stronger than the logarithmic
sensitivity obtained in the $T=0$ calculation and leads to a much larger
capture coefficient in the long-range regime. In the presence of target
thermal motion, the subspace of kinematic configurations that can lead to
capture at small momentum transfer is larger than the same subspace when
thermal motion is absent. In the $T=0$ calculation, the logarithmic scaling
with $M_{\gamma_{D}}$ is a consequence of the fact that capture with zero-
momentum transfer can occur only at one kinematic point,
$w=v_{\text{esc}}(r)$, for the incoming DM velocity. In the presence of
thermal motion of the targets, capture via zero-momentum transfer can occur in
the larger kinematic region
$w\in[v_{\text{esc}}(r),v_{\text{esc}}(r)+\mathcal{O}(\sqrt{3T_{\text{Sun}}/m_{i}})]$.
The forward singularity is now relevant in a finite (even if small) slice of
the kinematic integral rather than at a single point, thus promoting the
would-be logarithmic dependence on $\mu$ to a linear one. In turn, this leads
to the capture rate scaling as $\mu^{-1}\sim M^{-2}_{\gamma_{D}}$.

One important corollary of the discussion above is that the capture rate
calculation at finite $T$ remains insensitive to the high speed tail of the
halo DM distribution. Long-range gravitational capture is dominated by
scattering with small momentum transfer, which involves DM particles from the
small velocity tail of the halo distribution. We have checked numerically that
cutting the galactic frame DM velocity distribution off at the galactic escape
speed $v_{gal}\in[450\text{ km/s},650\text{ km/s}]$ [52] results in negligible
corrections to the capture rate. The dominance of small momentum transfer
events in the long-range capture rate also means that DM-nuclear scattering
will typically be well-approximated as coherent.

Finally, while we have used a constant core temperature for simplicity, the
Sun’s temperature varies with radius. Since most capture occurs near the core,
this is a reasonable first approximation. As the thermal velocity of the
target nuclei is important for controlling the thermal enhancement to the
long-range capture rate, we expect that accounting for the radial variation of
the solar temperature is more important for long-range capture than short-
range capture. We have checked numerically that incorporating the radial
temperature profile in the model of Ref. [46] reduces the finite-temperature
capture rate by up to tens of percent in the long-range regime, which can be
important for detailed studies of signals in any specific model realizing
long-range nuclear interactions.

##  3 Dark Matter self-capture and self-ejection with thermal corrections

In this section we discuss DM self-capture and self-ejection rates. Self-
capture occurs when an incoming DM particle scatters against a previously-
captured DM particle, such that both particles remain gravitationally bound
after the collision. Conversely, self-ejection occurs when both DM particles
are gravitationally unbound after scattering. The self-capture coefficient
$C_{\mathrm{sc}}$ gives the rate of gravitational self-capture per bound DM
particle, and similarly $C_{\mathrm{se}}$ for self-ejection. In our reference
DM model, the DM particle is Dirac; we assume the relic abundance of $\chi$
and $\bar{\chi}$ are equal. The calculation of $C_{\mathrm{sc}}$ proceeds
similarly to that for nuclear capture, except that the target particle is now
a previously-captured DM particle instead of a nucleus. The cross-section for
events leading to the net capture of DM is therefore subject to the additional
kinematic constraint that the target DM particle remain bound after the
collision; similarly, for self-ejection, an additional kinematic requirement
imposes that the target DM particle is unbound following the collision.

To compute the rate for either process, one must model the distribution of the
target DM population inside the capturing object. We consider the case where
the DM-SM interaction is strong enough to thermalize the captured DM particles
with the nuclei within the lifetime of the Sun.111See Ref. [53] for a
discussion of thermalization via short-range spin- and velocity-independent
nuclear interactions. The dark photon-mediated reference model used here
yields spin- and velocity-independent nuclear interactions, but depending on
the dark photon mass and the momentum transfer can produce either short- or
long-range interactions. The approach to thermalization involves multiple
scatterings with increasingly small momentum transfer; thus even if the
initial capture collision is in the long-range regime, the process of
thermalization may be dominated by short-range scatterings. The spatial
distribution of the DM is then given by the equilibrium quantity

$\displaystyle n_{c}(r)$
$\displaystyle=A\exp\left(-M_{\text{DM}}v^{2}_{\text{esc}}(r)/2T_{\text{Sun}}\right),$
(3.1)

which for simplicity we again take at the core temperature of the capturing
body. Here the normalization factor $A$ ensures $\int dVn_{c}(r)=1$. We will
again take the Sun’s temperature to be constant for simplicity. For self-
capture of heavy DM ($M_{\text{DM}}\gtrsim 10$ GeV), unlike nuclear capture,
this is a good approximation even in the long-range regime simply because the
physical volume occupied by the captured dark matter population is small,
occupying a sphere with radius at most $\mathcal{O}(2-3\%)$ of $R_{\odot}$,
where the variation of the solar temperature is small.

###  3.1 Zero-temperature calculation

We begin with a review of the self-capture calculation in the limit of no
thermal motion of the captured DM particles [54, 20]. The self-capture
coefficient is constructed in analogy to the nuclear capture coefficient, by
integrating contributions from all volume shells:

$\displaystyle\frac{dC_{\text{sc}}}{dV}=\frac{1}{2}n_{\text{DM}}n_{c}(r)\int_{0}^{\infty}d^{3}\vec{w}\Theta(|\vec{w}|-v_{\text{esc}}(r))f_{\eta}\left(\sqrt{w^{2}-v_{\text{esc}}^{2}(r)}\right)|\vec{w}|\sigma_{sc}.$
(3.2)

The differences with respect to Eq. 2.8 for the nuclear capture coefficient
are (i) the target particle distribution is now given by $n_{c}(r)$; (ii) an
overall factor of 1/2, stemming from the assumption that the relic abundance
of DM is symmetric; (iii) the relevant cross-section is now that for self-
capture, $\sigma_{sc}$.

At the microscopic level several processes contribute to $\sigma_{sc}$: self-
capture of incoming particles (antiparticles) by scattering against captured
particles (antiparticles), and of particles (antiparticles) by scattering on
target antiparticles (particles). In the non-relativistic limit the matrix
element squared, summing all contributions, can be written

$\displaystyle|\mathcal{M}|^{2}=\frac{4g^{4}_{D}}{v^{{}^{\prime}4}_{1}}f_{sc}(\cos\theta,\mu),$
(3.3)

with

$\displaystyle f_{sc}(\cos\theta,\mu)\equiv\frac{3}{(1-\cos\theta+\mu)^{2}}$
$\displaystyle+\frac{1}{(1+\cos\theta+\mu)^{2}}+\frac{1}{(1-\cos\theta+\mu)}\frac{1}{(1+\cos\theta+\mu)},$
(3.4)

where $\theta$ is the scattering angle in the CM frame and $\mu$ is defined in
Eq. 2.10. The self-capture condition demands

$\displaystyle|\vec{v}_{1,f}^{lab}|\leq v_{\text{esc}}(r)\quad\text{ and
}\quad|\vec{v}_{2,f}^{lab}|\leq v_{\text{esc}}(r).$ (3.5)

The Galilean transformation from the CM frame to the lab frame reads

$\displaystyle\vec{w}=\vec{v}^{\prime}_{1}+\vec{v}_{CM},\phantom{spsp}\vec{v}_{c}=0=-\vec{v}^{\prime}_{1}+\vec{v}_{CM}$
(3.6)

where $\vec{v}_{c}=0$ is a consequence of neglecting the thermal motion of the
target DM particles. Now we write the self-capture condition as a constraint
on the CM frame scattering angle $\theta$:

$\displaystyle\Theta_{sc}\equiv\Theta(\beta(r)-|\cos\theta|),$ (3.7)

noting that $\beta\geq 0$ is a necessary condition for self-capture to occur.
Furthermore, Eq. 2.23 simplifies to $\beta(r)=2v^{2}_{\text{esc}}(r)/w^{2}-1$,
ensuring that $\beta\leq 1$. The self-capture cross-section is found by
integrating the matrix element with $\Theta_{sc}$ specified in Eq. 3.7, over
the solid angle in the CM frame:

$\displaystyle\sigma_{sc}=\frac{1}{256\pi^{2}M^{2}_{\text{DM}}}\int
d\Omega_{f,\chi}|\mathcal{M}|^{2}\Theta_{sc}=\frac{4\pi\alpha^{2}_{D}}{M^{2}_{\text{DM}}w^{4}}\Theta(\beta(r))\int_{-\beta(r)}^{\beta(r)}d\cos\theta
f_{sc}(\cos\theta,\mu).$ (3.8)

The self-capture cross-section $\sigma_{sc}$ can be evaluated in closed form,

$\displaystyle\sigma_{sc}=\frac{4\pi\alpha^{2}_{D}}{M^{2}_{\text{DM}}w^{4}}\Theta(\beta(r))\left[\frac{8\beta(r)}{-\beta^{2}(r)+(1+\mu)^{2}}+\frac{2\text{ArcTanh}\left(\frac{\beta(r)}{1+\mu}\right)}{1+\mu}\right].$
(3.9)

Collecting terms, we obtain the following expression for self-capture per unit
volume:

$\displaystyle\frac{dC_{\text{sc}}}{dV}=\frac{1}{2}n_{\text{DM}}n_{c}(r)$
$\displaystyle\int_{v_{\text{esc}}(r)}^{\sqrt{2}v_{\text{esc}}(r)}\frac{dw}{w}4\pi
f_{\eta}\left(\sqrt{w^{2}-v_{\text{esc}}^{2}(r)}\right)\times$
$\displaystyle\times\frac{4\pi\alpha^{2}_{D}}{M^{2}_{\text{DM}}}\left[\frac{8\beta(r)}{-\beta^{2}(r)+(1+\mu)^{2}}+\frac{2\text{ArcTanh}\left(\frac{\beta(r)}{1+\mu}\right)}{1+\mu}\right].$
(3.10)

The term in the square brackets is singular at $w=v_{\text{esc}}(r)$ for
$\mu=0$, resulting in a logarithmic dependence of $C_{\text{sc}}$ on $\mu$,
and in turn on the dark photon mass $M_{\gamma_{D}}$, similar to zero-
temperature nuclear capture. The characteristic scattering regime is
determined by how large the mediator mass is compared to the typical momentum
transfer, $M_{\text{DM}}v_{\text{esc}}\sim(2\times 10^{-3})M_{\text{DM}}$:
heavy mediators correspond to pointlike interactions, giving rise to the
expected $M^{-4}_{\gamma_{D}}$ dependence; in the opposite limit, the
dependence is logarithmic.

###  3.2 Self-Capture at Finite Temperature

Now we turn to self-capture at finite temperature. We model the velocity
distribution of captured DM particles with a Maxwell distribution at the core
temperature, truncated at the local escape velocity:

$\displaystyle f_{c}(\vec{v}_{c})$
$\displaystyle=A_{\text{MB}}(r)\left(\frac{M_{\text{DM}}}{2\pi
T_{\text{Sun}}}\right)^{3/2}\exp{\left(-\frac{M_{\text{DM}}v_{c}^{2}}{2T_{\text{Sun}}}\right)},\qquad\qquad|\vec{v}_{c}|<v_{\text{esc}}(r),$
(3.11)

with $A_{\text{MB}}(r)$ chosen such that
$\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{c}f_{c}(\vec{v}_{c})=1$. The self-
capture coefficient now includes an integral over the initial velocity of the
captured DM particle, $\vec{v}_{c}$:

$\displaystyle\frac{dC_{\text{sc}}}{dV}$
$\displaystyle=\frac{1}{2}n_{\text{DM}}n_{c}(r)\int_{0}^{\infty}d^{3}\vec{u}\frac{w}{u}f_{\eta}(u)\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{c}f_{c}(v_{c})|\vec{w}-\vec{v}_{c}|\sigma_{sc}.$
(3.12)

The kinematic constraints leading to self-capture, given in Eq. 3.5, now
translate into the following step function on the scattering angle $\gamma$,
defined in Eq. 2.28:

$\displaystyle\Theta_{sc}\equiv\Theta_{sc}(\beta(r)-|\cos\gamma|)\times\Theta(\beta(r)).$
(3.13)

We make a Galilean transformation to bring the velocity integrals to the CM
frame, where $\vec{v}^{\prime}_{1}$ and $\vec{v}^{\prime}_{2}$ denote the CM
velocities of the incoming and captured DM particles before the collision,

$\displaystyle\int_{0}^{\infty}d^{3}\vec{u}\frac{w}{u}f_{\eta}(u)\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{c}f_{c}(v_{c})=\int_{v_{\text{esc}}(r)}^{\infty}d^{3}\vec{w}f_{\eta}\left(\sqrt{w^{2}-v_{\text{esc}}^{2}(r)}\right)\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{c}f_{c}(v_{c})=$
$\displaystyle=\int_{0}^{\infty}d^{3}\vec{w}f_{\eta}\left(\sqrt{w^{2}-v_{\text{esc}}^{2}(r)}\right)\Theta(|\vec{w}|-v_{\text{esc}}(r))\int_{0}^{\infty}d^{3}\vec{v}_{c}f_{c}(v_{c})\Theta(v_{\text{esc}}(r)-|\vec{v}_{c}|).$
(3.14)

Now we switch to CM frame velocities,

$\displaystyle\vec{w}=\vec{v}^{\prime}_{1}+\vec{v}_{CM},\phantom{spsp}\vec{v}_{c}=\vec{v}^{\prime}_{2}+\vec{v}_{CM}$
(3.15)

with $\vec{v}^{\prime}_{2}=-\vec{v}^{\prime}_{1}$. The Jacobian of the
transformation
$(\vec{w},\vec{v}_{c})\rightarrow(\vec{v}_{CM},\vec{v}^{\prime}_{1})$ is now
$J=8$. In the CM frame, the velocity integrals read

$\displaystyle\int_{0}^{\infty}d^{3}\vec{v}_{CM}\int_{0}^{\infty}d^{3}\vec{v}^{\prime}_{1}Jf_{\eta}\left(\sqrt{(\vec{v}^{\prime}_{1}+\vec{v}_{CM})^{2}-v_{\text{esc}}^{2}(r)}\right)f_{c}(|-\vec{v}^{\prime}_{1}+\vec{v}_{CM}|)\times$
$\displaystyle\times\Theta\left(|\vec{v}^{\prime}_{1}+\vec{v}_{CM}|-v_{\text{esc}}(r)\right)\Theta\left(v_{\text{esc}}(r)-|-\vec{v}^{\prime}_{1}+\vec{v}_{CM}|\right).$
(3.16)

We choose coordinates where $\vec{v}^{\prime}_{1}$ is along the $z$-axis, so
that the solid angle integration for $\vec{v}^{\prime}_{1}$ is trivial. We
denote the polar and azimuthal angles of $\vec{v}_{CM}$ as
$(\alpha,\phi_{CM})$. Using the definition of $\beta$ in Eq. 2.23 we write the
velocity integrals in a more compact form,

$\displaystyle\int_{0}^{\infty}v^{2}_{CM}dv_{CM}\int_{-1}^{1}dc\alpha\int_{0}^{2\pi}d\phi_{CM}\int_{0}^{\infty}4\pi
v^{{}^{\prime}2}_{1}dv^{\prime}_{1}Jf_{\eta}\left(\sqrt{2v_{CM}v^{\prime}_{1}(c\alpha-\beta(r))}\right)\times$
$\displaystyle\times
f_{c}\left(\sqrt{v^{2}_{CM}+v^{{}^{\prime}2}_{1}+2v_{CM}v^{\prime}_{1}c\alpha}\right)\Theta(c\alpha-\beta(r))\Theta(c\alpha+\beta(r)).$
(3.17)

We calculate the self-capture cross-section $\sigma_{sc}$ in a manner similar
to Eq. 3.8, where the matrix element given by Eq. 3.3 is convoluted with the
new self-capture condition in Eq. 3.13,

$\displaystyle\sigma_{sc}$
$\displaystyle=\frac{1}{256\pi^{2}M^{2}_{\text{DM}}}\int
d\cos\theta\,d\phi_{\chi}|\mathcal{M}|^{2}\Theta_{sc}=$
$\displaystyle=\frac{\alpha^{2}_{D}}{4M^{2}_{\text{DM}}(v_{1}^{\prime})^{4}}\Theta(\beta(r))\int
d\cos\theta d\phi_{\chi}f_{sc}(\cos\theta,\mu)\Theta(\beta(r)-|\cos\gamma|).$
(3.18)

Putting together this cross-section expression with the CM-frame velocity
integrals of Eq. 3.17 in Eq. 3.12 gives for the self-capture coefficient

$\displaystyle\frac{dC_{\text{sc}}}{dV}=\frac{1}{2}n_{\text{DM}}n_{c}(r)\int_{0}^{\infty}v^{2}_{CM}dv_{CM}\int_{0}^{\infty}v^{{}^{\prime}2}_{1}dv^{\prime}_{1}\int_{-1}^{1}dc\alpha(4\pi)^{2}f_{\eta}\left(\sqrt{2v_{CM}v^{\prime}_{1}(c\alpha-\beta(r))}\right)\times$
$\displaystyle\times
f_{c}\left(\sqrt{v^{2}_{CM}+v^{{}^{\prime}2}_{1}+2v_{CM}v^{\prime}_{1}c\alpha}\right)\Theta(c\alpha-\beta(r))\Theta(c\alpha+\beta(r))\Theta(\beta(r))\frac{2\alpha^{2}_{D}}{M^{2}_{\text{DM}}v_{1}^{{}^{\prime}3}}J_{sc}(\beta(r),c\alpha,\mu),$
(3.19)

where we introduced the angular integral $J_{sc}$,

$\displaystyle
J_{sc}\equiv\frac{1}{4\pi}\int_{0}^{2\pi}d\phi_{CM}\int_{0}^{2\pi}d\phi_{\chi}\int_{-1}^{1}d\cos\theta
f_{sc}(\cos\theta,\mu)\Theta(\beta(r)-|\cos\gamma|).$ (3.20)

We detail the calculation of $J_{sc}$ in Appendix B. The three step functions
in Eq. 3.19 restrict the integration domain as follows. The first and third
step functions make the second step function redundant. Additionally, the
first step function implicitly demands that $\beta(r)\leq 1$. Hence, we must
have $0\leq\beta(r)\leq 1$. Using the definition of $\beta$, this condition
turns into a finite integration region for the speeds:

$\displaystyle\frac{dC_{\text{sc}}}{dV}$
$\displaystyle=\frac{1}{2}n_{\text{DM}}n_{c}(r)\int_{0}^{v_{\text{esc}}(r)}v^{2}_{CM}dv_{CM}\int_{v_{\text{esc}}(r)-v_{CM}}^{\sqrt{v^{2}_{\text{esc}}(r)-v^{2}_{CM}}}v^{{}^{\prime}2}_{1}dv^{\prime}_{1}\int_{\beta(r)}^{1}dc\alpha\times$
$\displaystyle\times
f_{\eta}\left(\sqrt{2v_{CM}v^{\prime}_{1}(c\alpha-\beta(r))}\right)f_{c}\left(\sqrt{v^{2}_{CM}+v^{{}^{\prime}2}_{1}+2v_{CM}v^{\prime}_{1}c\alpha}\right)2\frac{(4\pi)^{2}\alpha^{2}_{D}}{M^{2}_{\text{DM}}v_{1}^{{}^{\prime}3}}J_{sc}(\beta(r),c\alpha,\mu).$
(3.21)

In Fig. 2 we present the results of numerically evaluating the remaining
integrals.

Figure 2: Left: Self-capture coefficient in the Sun as a function of mediator
mass $M_{\gamma_{D}}$. Solid lines correspond to the full calculation
including the target’s thermal motion, while dashed lines show the zero-
temperature calculation. The different colors correspond to different DM
masses. Dotted lines are included to highlight the parametric scaling with
$M_{\gamma_{D}}$ in the two regimes. Right: Self-capture coefficient in the
Sun, as a function of DM mass. We fix $M_{\gamma_{D}}=10^{-6}$ GeV,
corresponding to long-range scattering, and $\alpha_{D}=10^{-3}$.

The left panel shows the self-capture rate as a function of dark photon mass
for several different DM masses. The capture rate again exhibits two
qualitatively distinct regimes: the short-range regime with
$C_{\text{sc}}\propto M^{-4}_{\gamma_{D}}$ and the long-range regime where
$C_{\text{sc}}\propto M^{-2}_{\gamma_{D}}$, similar to long-range nuclear
capture. In the right panel, we plot the self-capture rate as a function of DM
mass, contrasting the full calculation to the zero temperature one. The
dependence of $C_{\text{sc}}$ on $M_{\text{DM}}$ comes from three sources: (i)
the local number density of DM in the halo; (ii) the scattering cross-section,
$\sigma_{sc}\propto M^{-2}_{\text{DM}}$; and the novel factor (iii)
$M^{2}_{\text{DM}}/M^{2}_{\gamma_{D}}$ from the $\mu^{-1}$ scaling of the
thermally-averaged self-capture cross-section. Combining these factors
produces a final $C^{T\neq 0}_{\text{sc}}\sim M^{-1}_{\text{DM}}$ scaling. By
contrast, in the zero-temperature calculation, the dependence on $\mu$ softens
to a logarithmic one, and to leading order, $C^{T=0}_{\text{sc}}\sim
M^{-3}_{\text{DM}}$.

Finally, we discuss the self-ejection coefficient $C_{\text{se}}$. Most of the
formalism and notation for this case carries over from self-capture and we
only highlight the qualitative differences between the two processes. First,
the self-ejection constraints on the final state speeds are reversed from
those for self-capture,

$\displaystyle|\vec{v}^{lab}_{1,f}|\geq
v_{\text{esc}}(r),\quad\quad|\vec{v}^{lab}_{2,f}|\geq v_{\text{esc}}(r).$
(3.22)

These conditions lead to the following step functions ensuring self-ejection:

$\displaystyle\Theta_{se}\equiv\Theta(-\beta(r)-|\cos\gamma|)\Theta(-\beta(r))\Theta(1+\beta(r)).$
(3.23)

Recall that self-capture demands $0\leq\beta(r)\leq 1$, while for self-
ejection $-1\leq\beta(r)\leq 0$. We minimize the computational cost of the
numerical integrals in this case by observing that the contribution to self-
ejection at $\bar{\beta}\equiv-\beta$ is equal to the contribution to self-
capture at $\beta$. This observation is a consequence of the fact that the
same microscopic scattering processes are responsible for both self-capture
and self-ejection, and the difference between the two processes is merely
kinematic. Then,

$\displaystyle\sigma_{se}(\beta(r),c\alpha,\mu)=\sigma_{sc}(\bar{\beta}(r),c\alpha,\mu).$
(3.24)

Therefore, we can construct the self-ejection coefficient per unit volume
using Eq. 3.19 by exploiting this symmetry:

$\displaystyle\frac{dC_{\text{se}}}{dV}=\frac{1}{2}n_{\text{DM}}n_{c}(r)\int_{0}^{\infty}v^{2}_{CM}dv_{CM}\int_{0}^{\infty}v^{{}^{\prime}2}_{1}dv^{\prime}_{1}\int_{-1}^{1}dc\alpha\Theta(1-\bar{\beta}(r))\Theta(\bar{\beta}(r))\Theta(c\alpha+\bar{\beta}(r))\times$
$\displaystyle\times
f_{\eta}\left(\sqrt{2v_{CM}v^{\prime}_{1}(c\alpha+\bar{\beta}(r))}\right)f_{c}\left(\sqrt{v^{2}_{CM}+v^{{}^{\prime}2}_{1}+2v_{CM}v^{\prime}_{1}c\alpha}\right)2\frac{(4\pi)^{2}\alpha^{2}_{D}}{M^{2}_{\text{DM}}v_{1}^{{}^{\prime}3}}J_{sc}(\bar{\beta}(r),c\alpha,\mu).$
(3.25)

At zero temperature, self-ejection is negligible compared to self-capture in
the Sun for both long- and short-range interactions [54, 20, 53]. We have
verified that this remains true for $T\neq 0$.

##  4 Implications for the Captured Dark Matter Population

As demonstrated in the two previous sections, retaining the finite thermal
motion of the target particles is critical in correctly evaluating capture
rates through long-range interactions ($\mu\ll 1)$. To help understand the
range of mediator masses where these thermal target effects are quantitatively
significant for a given astrophysical object, we begin by estimating the scale
of momentum transfer in DM-nucleus collisions. Using Eqs. 2.20 and 2.21, the
momentum exchange in a given scattering between halo DM particles and a
nuclear target is

$\displaystyle(\Delta\vec{p})^{2}$
$\displaystyle=M^{2}_{\text{DM}}v_{1}^{{}^{\prime}2}=M^{2}_{\text{DM}}|\vec{w}-\vec{v}_{CM}|^{2}=\frac{M^{2}_{\text{DM}}m_{i}^{2}}{(M_{\text{DM}}+m_{i})^{2}}|\vec{w}-\vec{v}_{n}|^{2}.$
(4.1)

Now, $v_{n}$ is drawn from a Maxwell distribution at the (core) temperature of
the massive body, hence $v^{2}_{n}\sim 3T/m_{i}$. Meanwhile,
$w^{2}=u^{2}+v^{2}_{\text{esc}}(r)$, where the natural scale for $u$ is
$\bar{v}$, the rms speed of the halo particles. Hence, the typical momentum
transfer in a scattering event can be estimated as

$\displaystyle\Delta
p_{typ}=\bar{m}_{i}\sqrt{v^{2}_{\text{esc}}(r)+\bar{v}^{2}+3T/m_{i}-2\sqrt{v^{2}_{\text{esc}}(r)+\bar{v}^{2}}\sqrt{3T/m_{i}}\cos(\measuredangle(\vec{w},\vec{v}_{n}))},$
(4.2)

up to factors of a few. In the Sun, the three velocity scales that control the
kinematics are

$\displaystyle\bar{v}=288\text{ km/s},$ $\displaystyle
v_{\text{esc}}(r)\in[600,1300]\text{ km/s},$
$\displaystyle\sqrt{3T_{\text{Sun}}/m_{i}}\sim 85\text{
km/s}\times\left(\frac{50\text{ GeV}}{m_{i}}\right)^{1/2}.$ (4.3)

Thus for scattering in the Sun, the dominant velocity scale is
$v_{\text{esc}}$, giving

$\displaystyle\Delta p_{typ}\approx 5\times 10^{-3}\bar{m}_{i}.$ (4.4)

For collisions with oxygen nuclei, $\Delta p_{typ}=16$ MeV for heavy dark
matter, $M_{\text{DM}}\gg 16$ GeV. For $M_{\text{DM}}\ll 16$ GeV, $\Delta
p_{typ}=1\,\text{MeV}(M_{\text{DM}}/\text{GeV})$. Thus for gravitational
capture in the Sun, typical DM-nuclear collisions are effectively short-range
for mediator masses heavier than the 10 MeV scale.

Electrons have the same temperature as nuclei, but owing to their smaller mass
will have a much larger thermal velocity. Thus for electrons the thermal
velocity can easily dominate over the other velocity scales in evaluating
$\Delta p_{typ}$. However, for DM heavier than the evaporation mass in the
Sun, the reduced DM-electron mass is $\bar{m}_{i}\approx m_{e}$, which
suppresses $\Delta p_{typ}$ in comparison to nuclear scattering. We thus
expect scattering on electrons to remain subleading in this DM mass range.

Meanwhile the typical momentum exchange in collisions between halo and bound
DM, relevant for self-capture, follows by replacing the DM-nucleus reduced
mass $\bar{m}_{i}$ with the reduced mass for the DM-DM system. Thus while the
typical momentum transfer in nuclear capture collisions saturates as the DM
mass increases, in the sense that $(\Delta
p_{typ})^{\text{c}}=\bar{m}_{i}v_{typ}\approx m_{i}v_{typ}$, the typical
momentum transfer in self-capture collisions $(\Delta
p_{typ})^{\text{sc}}=1/2M_{\text{DM}}v_{typ}$ increases with DM mass.
Consequently, for heavy DM, long-range self-capture rates can exhibit a
parametric enhancement for mediator masses where nuclear capture rates do not.

For mediator masses $M_{\gamma_{D}}\lesssim\Delta p_{typ}$, the short-range
approximation is no longer adequate to describe gravitational capture.
However, as we have seen, the biggest thermal enhancements to capture rates
arise from the forward scattering regime, where the momentum transfers can be
substantially smaller than $\Delta p_{typ}$. To understand precisely how large
the hierarchy $M_{\gamma_{D}}\ll\Delta p_{typ}$ needs to be for thermal
enhancements to be sizeable, we numerically evaluate the capture rates and
find the resulting captured dark matter population.

In the regime where DM reaches thermal equilibrium with the nuclei of the Sun,
the number of bound DM particles in the Sun, $N(t)$, can be obtained in terms
of the nuclear and self-capture coefficients, together with the annihilation
coefficient $C_{\text{ann}}$, by solving

$\frac{dN}{dt}=C_{\text{c}}+C_{\text{sc}}N-C_{\text{ann}}N^{2}.$ (4.5)

Here for brevity we have omitted subdominant processes such as self-ejection
and (self-)evaporation, which can be neglected for DM masses above $\sim 4$
GeV ($\sim 5$ GeV). In the generic situation where DM annihilations proceed
dominantly to light mediators, the annihilation coefficient $C_{\text{ann}}$
is given by

$\displaystyle C_{\text{ann}}$
$\displaystyle=\frac{1}{2}\langle\sigma_{\text{ann}}v\rangle\int
dVn^{2}_{c}(r)$ (4.6)
$\displaystyle=\frac{1}{2}\frac{\pi\alpha_{D}^{2}}{M^{2}_{\text{DM}}}\langle
S_{S}\rangle\int dVn^{2}_{c}(r),$ (4.7)

where $\langle S_{S}\rangle$ is the thermally-averaged Sommerfeld enhancement,
with $\langle S_{S}\rangle\propto\alpha_{D}M^{1/2}_{\text{DM}}$ [20], and
$n_{c}(r)$ is given by Eq. 3.1.

The solution to Eq. 4.5 is [54]

$\displaystyle
N(t)=\frac{C_{\text{c}}\text{tanh}\left(\frac{t}{\xi}\right)}{\frac{1}{\xi}-\frac{C_{\text{sc}}}{2}\text{tanh}\left(\frac{t}{\xi}\right)},$
(4.8)

in terms of the equilibration time scale

$\xi^{-1}=\sqrt{C_{\text{c}}C_{\text{ann}}+C^{2}_{\text{sc}}/4}.$ (4.9)

Parametrically, for heavy DM in the long-range regime, the nuclear capture
rate scales as

$\displaystyle
C_{\text{c}}\propto\frac{\epsilon^{2}\alpha_{D}}{M^{2}_{\text{DM}}M^{2}_{\gamma_{D}}},$
(4.10)

while for self-capture

$\displaystyle
C_{\text{sc}}\propto\frac{\alpha_{D}^{2}}{M_{\text{DM}}M^{2}_{\gamma_{D}}}.$
(4.11)

Meanwhile, since $\int dVn^{2}_{c}(r)\propto M^{3/2}_{\text{DM}}$,
$C_{\text{ann}}\propto\alpha_{D}^{3}$. Thus
$C_{\text{c}}C_{\text{ann}}\propto\epsilon^{2}\alpha_{D}^{4}M^{-2}_{\text{DM}}M^{-2}_{\gamma_{D}}$
and
$C_{\text{sc}}^{2}\propto\alpha_{D}^{4}M^{-2}_{\text{DM}}M^{-4}_{\gamma_{D}}$.
Dialing the dark gauge coupling does not affect the relative importance of
nuclear capture versus self-capture, since both $C_{\text{sc}}^{2}$ and
$C_{\text{c}}C_{\text{ann}}$ have the same dependence on $\alpha_{D}$ [33].
However, since at finite temperature $C_{\text{c}}C_{\text{ann}}$ and
$C_{\text{sc}}^{2}$ have different parametric dependences on $M_{\text{DM}}$
and $M_{\gamma_{D}}$, the importance of self-capture compared to nuclear
capture increases as the hierarchy between the dark photon and DM mass is
increased.

In the rest of this section we examine the finite-temperature DM population in
two different regimes for the mass of the dark photon.

###  4.1 Visibly decaying dark photons in the Sun and beyond

Figure 3: Ratio of the number of captured DM particles in the Sun at finite
temperature to the zero-temperature result, $N^{T\neq
0}_{\text{c}}/N^{T=0}_{\text{c}}$, as a function of dark gauge boson mass
$M_{\gamma_{D}}$. The solid (dashed) curves correspond to $M_{\text{DM}}=100$
$(1000)$ GeV. The different colors correspond to different values of kinetic
mixing parameter $\epsilon$. Shaded blue (red) regions are ruled out by
stellar cooling and supernova energy loss constraints for
$\epsilon=10^{-9}(10^{-10})$; for $\epsilon=10^{-7}$ and $\epsilon=10^{-8}$
these constraints plus those from the experiments LSND and E137 rule out the
entire mass range.

In Fig. 3 we show the ratio of the finite temperature calculation of the total
bound population today $N_{\text{c}}\equiv N(\tau_{\text{Sun}})$ in the Sun to
the zero-temperature calculation, as a function of $M_{\gamma_{D}}$.222Since
both $C_{\text{c}}C_{\text{ann}}$ and $C_{\text{sc}}^{2}$ are proportional to
the same power of $\alpha_{D}$, all dependence on $\alpha_{D}$ drops out of
this population ratio. Here $\tau_{\text{Sun}}=5\times 10^{9}$ years is the
age of the Sun. For $M_{\gamma_{D}}\gg 10$ MeV, capture collisions are short-
range and Fig. 3 accordingly displays the expected small negative thermal
corrections. Appreciable positive thermal corrections appear for
$M_{\gamma_{D}}\lesssim$ MeV, more than an order of magnitude below $\Delta
p_{typ}\sim 16$ MeV. For $\epsilon=10^{-7}$, self-capture makes negligible
contributions to the DM abundance, and the thermal enhancement exhibited for
$M_{\gamma_{D}}\lesssim$ MeV is due entirely to the enhanced nuclear capture
rate. As $\epsilon$ is decreased, self-capture plays a more important role in
determining the number of captured DM particles. Since long-range self-capture
also becomes more important at larger DM masses, we see the differences
between the $M_{\text{DM}}=100$ GeV (solid) and $M_{\text{DM}}=1$ TeV (dashed)
curves increase with decreasing $\epsilon$. Thus the differences between the
dashed and solid lines in Fig. 3 at given values of $\epsilon$ and
$M_{\gamma_{D}}$ help establish the relative importance of self-capture. We
also discuss the regimes where nuclear capture or self-capture dominates in a
broader range of parameter space in Fig. 4 below.

For the values of $\epsilon$ shown in Fig. 3, dark photons produced in the
annihilations of captured DM escape the Sun before decaying [3, 28]. Dark
photons with masses above the MeV scale can decay to charged particle-
antiparticle pairs outside the Sun, giving cosmic ray signals that can be
probed by AMS [33] and gamma-ray telescopes such as Fermi and HAWC [3, 31, 35,
36, 55, 38, 39]. In the parameter space relevant to AMS,
$10^{-10}\lesssim\epsilon\lesssim 10^{-6}$ and $M_{\gamma_{D}}\in[1\text{
MeV},20\text{ MeV}]$, thermal corrections modify the $T=0$ approximation to
the nuclear capture rate at the $\text{few }\%$ level. For these benchmark
choices of $\alpha_{D}$ and $M_{\mathrm{DM}}$ and setting $M_{\gamma_{D}}=100$
MeV, which corresponds to short-range nuclear capture, the values of the
kinetic mixing parameter $\epsilon$ shown in Fig. 3 are at least an order of
magnitude above the smallest coupling required for DM to thermalize with solar
nuclei following the estimate of Ref. [53]. For fixed $\epsilon$, going to
lower mediator masses leads to an overall larger cross-section, but for sub-
MeV mediators capture collisions are now in the transitional regime between
long-range and short-range. Thus in this mediator mass range the first few
collisions following capture may not be fully in the short-range regime.
However, the subsequent nuclear scatterings in the approach to thermalization
see successively smaller momentum transfers, and are thus well-described by
short-range interactions. The DM-nuclear cross-sections for the parameter
points shown in Fig. 3 are well above the lower limits for short-range
thermalization, and thus we expect the thermal DM population to be a self-
consistent description.

In general, we expect from Fig. 3 that visibly-decaying mediators are not
sufficiently light to realize large thermal enhancements to DM capture in the
Sun. Sub-MeV mass mediators, which can give rise to sizeable thermal
enhancements, are subject to stringent constraints. In Fig. 3 we show the
constraints from stellar cooling [56, 57, 58] and rapid energy loss in
Supernova 1987A [59] for $\epsilon=10^{-9}$ and $\epsilon=10^{-10}$ (for
$\epsilon=10^{-7}$ and $10^{-8}$ the entire mass range is ruled out after
imposing these constraints together with the exclusions from E137 [60, 61] and
LSND [62, 63]). For $\epsilon=10^{-9}$ and $\epsilon=10^{-10}$ there are
additional constraints on the minimal dark photon model from extragalactic
background light [64, 65] and supernovae [66] that can be circumvented in the
presence of additional invisible decay modes for the dark photon. In the
regions of dark photon parameter space shown in Fig. 3, stellar cooling and
supernova energy loss constraints arise from the free-streaming of dark
photons out of the star, and thus are not alleviated by invisible dark photon
decay modes.

While for concreteness we have focused on dark photon mediators, stellar
cooling constraints appear at similar mass scales for other possible mediators
[67, 68], and thus these general conclusions apply more broadly. However, as
Fig. 3 makes clear, sizeable thermal enhancements in solar capture begin to
appear a relatively small distance below the exclusions from stellar cooling
bounds. Thus in hotter, more massive stars, characterized by larger values of
$\Delta p_{typ}$, we can expect appreciable thermal enhancements to the
capture of DM via MeV-scale mediator exchange. Gravitational capture of DM in
Population III stars is especially interesting in this regard.

Population III stars are thought to form from the collapse of primordial gas
clouds at the centers of DM minihaloes. These first stars can have masses
ranging up to $1000M_{\odot}$, and are born in a DM-dense environment [69,
70]. The rate of (short-range) DM gravitational capture in these massive early
stars can potentially be large enough to give rise to detectable impacts on
the properties of the first-forming stars [71, 70]; see also [72]. Given the
large escape velocities within these massive stars, we expect thermal
enhancements to long-range capture to be appreciable in these systems even for
visibly-decaying mediators, and important for understanding their potential
evolution. Generally, the capture rate will be dominated by the central
regions of the star where the escape velocity and temperature are both
maximized; a model for the internal structure of Population III stars is
therefore important for establishing the quantitative impacts of long-range
capture in these objects.

###  4.2 Solar capture through ultralight dark photons

For dark photons, as distinct from other mediators, stellar cooling
constraints become weaker as the dark photon mass decreases [57]. Working in a
basis of mass eigenstates, the suppression of dark photon emission with
decreasing $M_{\gamma_{D}}$ from the plasma in the stellar core can be
understood as an in-medium suppression of the effective mixing angle with the
photon by the factor $M_{\gamma_{D}}^{2}/(M_{\gamma_{D}}^{2}-m_{A}^{2})$,
where $m_{A}$ is the effective mass for the photon (see, e.g., [68]). The
consequent suppression of the low-mass dark photon emission rate from stars
opens up a distinct regime where gravitational capture can be mediated by very
light (sub-eV) dark photons. In this regime gravitational capture in the Sun
does become long-range and can thus exhibit substantial thermal enhancements.

Figure 4: Parameter space for light dark photons and gravitational capture of
DM in the Sun. The gray region shows the constraints on minimal dark photons
from stellar cooling [56, 57, 58], Xenon10 [73], supernova cooling [59] and
$e^{+}e^{-}$ pair production [66], BBN [74], CMB spectral distortions [75],
the intergalactic diffuse photon background [64, 65], the E137 [60, 61] and
LSND [62, 63] experiments, and fifth force experiments [76, 77]. The cyan line
shows the solar Debye mass $m_{D}$, while blue lines indicate where the dark
photon has a lifetime comparable to the age of the universe (solid), the
radius of Earth’s orbit (dashed), and the Solar radius (dotted). Other
contours indicate aspects of nuclear capture and self-capture of DM for the
reference values $M_{\text{DM}}=100$ GeV and $\alpha_{D}=10^{-3}\,(10^{-4})$
with solid (dashed) lines, as follows: bright purple lines indicate where the
geometric limit for finite-temperature nuclear capture is saturated, while the
pale purple contour indicates the geometric limit in the zero temperature
calculation. Yellow lines indicate where the geometric limit for finite-
temperature DM self-capture is saturated; the corresponding geometric limits
in the zero-temperature calculation appear at values of $\epsilon$ above the
shown parameter space. The bright (pale) orange contours indicate where the
equilibration timescale $\xi$, Eq. 4.9, is equal to the age of the Sun at
finite (zero) temperature. Above the bright green line, nuclear capture
dominates the captured population dynamics and provides the dominant mechanism
of gravitational capture, while below the line self-capture dominates. The
pale green contour shows the analogous result in the zero-temperature
calculation.

For long-range scattering with $M_{\gamma_{D}}\gg m_{A}$, the parametric
dependence of the capture coefficient $C_{c}$ on the model parameters
$(M_{\text{DM}},M_{\gamma_{D}},\epsilon,\alpha_{D})$ is given in Eq. 4.10. We
determine the coefficient of proportionality numerically. For dark photons
with $M_{\gamma_{D}}\ll m_{A}$, we can incorporate the in-medium suppression
of the SM coupling with the replacement
$\epsilon^{2}/M_{\gamma_{D}}^{2}\rightarrow\epsilon^{2}/m_{A}^{2}$. This
scaling can be understood in the mass eigenstate basis by noting that within
the solar plasma, the mostly-SM photon eigenstate picks up a coupling to the
dark matter current given by $\epsilon
g_{D}m_{A}^{2}/(M_{\gamma_{D}}^{2}-m_{A}^{2})$, which becomes $-\epsilon
g_{D}$ in the limit $M_{\gamma_{D}}\ll m_{A}$. Scattering through the mostly-
SM eigenstate is then unsuppressed by the small ratio
$(M_{\gamma_{D}}/m_{A})$, and gives the leading contribution to the capture
rate in the light dark photon limit.

We estimate the size of the resulting IR cutoff by considering the leading
contribution to long-range scattering in the nonrelativistic limit. This
contribution comes from the longitudinal portion of the photon propagator and
in the nonrelativistic limit reduces to capture via an electric field. The
Debye mass, related to the plasma frequency by
$m_{D}^{2}=\omega_{p}^{2}(m/T)$, controls the long-range screening of electric
fields, and has also been discussed as a regulating mass scale for long-range
DM-SM scattering in [78].

We will assume, for simplicity, that all the nuclear species in the Sun are
ionized and impose local charge neutrality. Then the plasma frequency is

$\displaystyle\omega_{p}^{2}=\sum_{i}\frac{q_{i}^{2}n_{i}}{m_{i}}\approx\frac{e^{2}n_{e}}{m_{e}},$
(4.12)

where the electron density in the Sun is

$\displaystyle n_{e}(r)=\sum_{i}Z_{i}\frac{\rho_{i}(r)}{m_{i}}.$ (4.13)

Again, for simplicity we take the temperature and number densities at the core
of the Sun, where most capture occurs. Numerically this yields
$\omega_{p}=0.26$ keV and $m_{D}=5$ keV. For
$M_{\gamma_{D}}\gg\omega_{p},m_{D}$ the long-range capture rate scales as
$\epsilon^{2}/M_{\gamma_{D}}^{2}$; in the light dark photon parameter space,
where $M_{\gamma_{D}}\ll\omega_{p},m_{D}$, the capture rate becomes
independent of the dark photon mass, scaling instead as
$\epsilon^{2}/m_{D}^{2}$.

This scaling will hold up until reaching values of $\epsilon$ large enough to
saturate the geometric limit. The geometric limit corresponds to the maximum
physical capture rate where every DM particle incident on the Sun is captured.
We approximate this upper limit $C^{geom}_{\text{c}}$ as the incident flux of
halo DM on the cross-sectional area of the Sun:

$\displaystyle C^{geom}_{\text{c}}=\pi
R^{2}\frac{\rho_{\text{DM}}}{M_{\text{DM}}}\int
d^{3}\vec{u}\frac{w(R)}{u}f_{\eta}(u)|\vec{w}(R)|.$ (4.14)

For fixed values of $M_{\text{DM}}$ and $\alpha_{D}$, setting

$\displaystyle
C_{\text{c}}(M_{\text{DM}},M_{\gamma_{D}},\epsilon,\alpha_{D})=C^{geom}_{\text{c}}$
(4.15)

then determines the curve in the $(M_{\gamma_{D}},\epsilon)$ plane where the
geometric limit is saturated. We show the resulting estimate for the
saturation of the capture cross-section in Fig. 4 for a 100 GeV DM particle,
with bright solid (dashed) purple lines indicating
$\alpha_{D}=10^{-3}\,(10^{-4})$. Above the bright purple line, the Sun is
optically thick to DM and the capture cross-section has saturated; below the
purple line, the capture rate is proportional to $(\epsilon/m_{D})^{2}$. For
comparison, we show in pale purple the geometric saturation curve from the
zero-temperature calculation for $\alpha_{D}=10^{-3}$.

Models featuring very light dark photons are restricted by constraints on the
effective number of relativistic species during recombination,
$N_{\mathrm{eff}}$. These constraints can be ameliorated by simple models of
late-time asymmetric reheating [79, 80], but the condition that the dark
photon does not thermalize with the SM at late times requires the kinetic
mixing to satisfy $\epsilon\lesssim 10^{-9}$ [81, 82]. In this regime, the
captured DM population will certainly thermalize with itself, but may or may
not thermalize with the nuclei of the capturing body. Thus the temperature of
the captured DM may in general differ from that of the nuclei, but we
nonetheless expect a thermal distribution for both populations to remain an
excellent description.

For simplicity we will consider the interplay between capture and self-capture
assuming that the captured DM has thermalized with the nuclei at
$T_{\text{Sun}}$. We expect this to be a conservative assumption in the sense
that when a captured DM particle is unable to efficiently transfer momentum to
nuclei, the captured DM ball will have a higher temperature than the SM, and
thus self-capture will become even more important relative to nuclear capture.

Let us first consider the geometric limit of self-capture. Letting $r_{x}$
denote the radius of the captured DM ball, adapting Eq. 4.14 to self-capture
collisions gives

$\displaystyle C_{\text{sc}}^{geom}=\pi
r_{x}^{2}\frac{\rho_{\text{DM}}}{M_{\text{DM}}}\int
d^{3}\vec{u}\frac{w(r_{x})}{u}f_{\eta}(u)|\vec{w}(r_{x})|.$ (4.16)

We define $r_{x}$ as the radius of a sphere enclosing 95% of the captured DM
particles:

$\displaystyle 0.95=\int_{0}^{r_{x}}dVn_{c}(r).$ (4.17)

For $M_{\text{DM}}=100\text{ GeV}$, $r_{x}=0.02R$. The total self-capture rate
$C_{\text{sc}}N_{c}$ is bounded from above by the geometric self-capture rate:

$\displaystyle
C_{\text{sc}}(M_{\text{DM}},M_{\gamma_{D}},\alpha_{D})N_{c}(M_{\text{DM}},M_{\gamma_{D}},\alpha_{D},\epsilon)\leq
C_{\text{sc}}^{geom},$ (4.18)

where the right-hand side only depends on $M_{\text{DM}}$ and $N_{c}$ is given
by Eq. 4.8. In Fig. 4 yellow solid (dashed) lines show the curves in the
$(M_{\gamma_{D}},\epsilon)$ plane that saturate this finite-temperature
geometric limit for fixed $M_{\text{DM}}=100$ GeV and
$\alpha_{D}=10^{-3}\,(10^{-4})$. In the parameter space of interest, the
geometric limit is saturated at values of the dark photon masses that are much
larger than regulators from in-medium effects.

Also in Fig. 4, we show in orange where the equilibration time scale $\xi$ is
equal to $\tau_{\text{Sun}}$; above the orange curves, the DM population has
reached a steady-state value. Green lines indicate where
$C_{\text{c}}C_{\text{ann}}=C_{\text{sc}}^{2}/4$; above (below) the green
lines, $C_{\text{c}}C_{\text{ann}}\gg C_{\text{sc}}^{2}/4$
($C_{\text{c}}C_{\text{ann}}\ll C_{\text{sc}}^{2}/4$) and nuclear capture
(self-capture) is the dominant mechanism for gravitational DM capture.

The finite-temperature calculation shows three important differences from the
zero-temperature calculation. First, vast regions of parameter space are in
the geometric limit for self-capture, a more drastic change than exhibited for
the case of nuclear capture. An immediate consequence is that the much larger
self-capture rate is sufficient to equilibrate the captured population at much
smaller values of $\epsilon$ than the zero-temperature calculation would
suggest. Finally, owing to the different parametric scalings of capture and
self-capture at finite temperature, self-capture dominates over nuclear
capture in an expanded range of parameter space.

Ultra-light dark photons are cosmologically long-lived. In this light dark
photon regime, DM annihilation within massive bodies therefore produces no
visible signals from mediator decays; however, the potentially very large
captured DM population may leave an observable imprint in stellar structure or
energy transport. Concrete realizations of this scenario have been considered
in the context of mirror stars and give rise to astrophysically interesting
signatures and constraints already in the zero-temperature regime [27]; we
expect the thermal enhancements we find here will make such probes of mirror
dark matter even more powerful.

##  5 Conclusions

We have reexamined the gravitational capture of DM by non-degenerate
astrophysical bodies at finite temperature. When capture occurs via long-range
interactions, we find that retaining the thermal motion of the target
particles leads to a capture rate that is quadratically dependent on the mass
of the mediator. This result gives a potentially sizable enhancement to the DM
capture rate compared to expectations from the zero-temperature calculation,
which predicts a logarithmic sensitivity to the mediator mass. This parametric
difference between the finite-temperature and zero-temperature results can be
simply understood by observing that the small volume of phase space that
results in capture via arbitrarily soft scattering at finite target
temperature shrinks to a single point in the zero-temperature limit. As a
concrete example, we provide numerical calculations for the solar capture of a
Dirac DM particle interacting with the SM via a kinetically mixed dark photon,
and study the transition from short-range to long-range capture regimes.

Numerically we find that substantial thermal enhancements to the DM capture
rate begin to appear once the mediator mass is more than an order of magnitude
smaller than the typical momentum transfer in DM-nucleus collisions. Since we
estimate the typical momentum transfer scale to be $\Delta p_{typ}\sim$ 16 MeV
in the Sun, sizeable thermal enhancements to solar capture are not compatible
with visible decays of the mediator. However, for light dark photons,
substantial enhancements can be realized, in nuclear capture but particularly
in self-capture. We study the interplay between nuclear capture and self-
capture in the Sun and demonstrate that in almost all of the relevant
parameter space, the captured DM population becomes optically thick to self-
capture. Thus we expect stellar capture of mirror dark matter to be
substantially more efficient than previously understood. Even for MeV-scale
mediators, we expect the thermal enhancement to long-range capture to be
important for understanding DM capture in the first generation of stars, which
are substantially more massive than the Sun.

### Acknowledgements

We thank David Curtin for useful conversations. The work of CG and JS is
supported in part by DOE Early Career grant DE-SC0017840.

##  Appendix A The angular integral $J_{c}$

In this appendix we present the details of the analytical and numerical
calculation of the angular integral $J_{c}$, defined in Eq. 2.29:

$\displaystyle
J_{c}\equiv\frac{1}{4\pi}\int_{0}^{2\pi}d\phi_{CM}\int_{-1}^{1}d\cos\theta\int_{0}^{2\pi}d\phi_{\chi}\frac{1}{(1-\cos\theta+\mu)^{2}}\Theta(\beta(r)-\cos\gamma),$
(A.1)

with
$\cos\gamma=\sin\alpha\sin\theta\cos(\phi_{\chi}-\phi_{CM})+\cos\alpha\cos\theta$.
Note that $\cos\gamma$ depends on azimuthal angles only through the difference
$\phi_{\chi}-\phi_{CM}$. Switching variables to
$\phi_{\pm}=\phi_{CM}\pm\phi_{\chi}$, we write

$\displaystyle J_{c}$
$\displaystyle=\frac{1}{4\pi}\int_{-1}^{1}d\cos\theta\frac{1}{(1-\cos\theta+\mu)^{2}}\int
d\phi_{+}\int
d\phi_{-}\frac{1}{2}\Theta\left(2\pi-\frac{1}{2}(\phi_{+}+\phi_{-})\right)\times$
$\displaystyle\times\Theta\left(2\pi-\frac{1}{2}(\phi_{+}-\phi_{-})\right)\Theta\left(\frac{1}{2}(\phi_{+}+\phi_{-})\right)\Theta\left(\frac{1}{2}(\phi_{+}-\phi_{-})\right)\Theta(\beta(r)-\cos\gamma),$
(A.2)

where the new factor of 1/2 comes from the Jacobian of the transformation and
the first four step functions impose the original bounds of integration
$\phi_{CM},\phi_{\chi}\in[0,2\pi)$. The latter two of the aforementioned step
functions demand that $\phi_{+}>0$. The product of the first four step
functions in Eq. A.2 can be turned into upper and lower limits for the
$\phi_{\pm}$ integration,

$\displaystyle\int_{0}^{\infty}d\phi_{+}\int_{0}^{\infty}$ $\displaystyle
d\phi_{-}\Theta\left(4\pi-(\phi_{+}+\phi_{-})\right)\Theta\left(4\pi-(\phi_{+}-\phi_{-})\right)\Theta\left(\phi_{+}+\phi_{-}\right)\Theta\left(\phi_{+}-\phi_{-}\right)=$
$\displaystyle=\int_{0}^{\infty}d\phi_{-}\int_{\phi_{-}}^{4\pi-\phi_{-}}d\phi_{+}\Theta(4\pi-2\phi_{-})=\int_{0}^{2\pi}d\phi_{-}(4\pi-2\phi_{-}).$
(A.3)

We further exploit the symmetry of this integral by partitioning the
integration region into three pieces, $0\leq\phi_{-}<\pi/2$,
$\pi/2\leq\phi_{-}<3\pi/2$, and $3\pi/2\leq\phi_{-}<2\pi$. In the last region
we shift $\phi_{-}\rightarrow\phi_{-}-2\pi$, under which
$\cos\phi_{-}\rightarrow\cos\phi_{-}$ and $\cos\gamma\rightarrow\cos\gamma$.
In the second region, consider $\phi_{-}\rightarrow\phi_{-}-\pi$, which
implies $\cos\phi_{-}\rightarrow-\cos\phi_{-}$ and
$\cos\gamma\rightarrow\cos\bar{\gamma}\equiv-\sin\alpha\sin\theta\cos\phi_{-}+\cos\alpha\cos\theta$.
Then,

$\displaystyle\int_{0}^{2\pi}d\phi_{-}(4\pi-2\phi_{-})$
$\displaystyle\Theta(\beta(r)-\cos\gamma)=\int_{0}^{\pi/2}d\phi_{-}(4\pi)\Theta(\beta(r)-\cos\gamma)+$
$\displaystyle+\int_{0}^{\pi/2}d\phi_{-}(-2\phi_{-})\Theta(\beta(r)-\cos\gamma)+\int_{-\pi/2}^{\pi/2}d\phi_{-}(2\pi)\Theta(\beta(r)-\cos\bar{\gamma})+$
$\displaystyle+\int_{-\pi/2}^{\pi/2}d\phi_{-}(-2\phi_{-})\Theta(\beta(r)-\cos\bar{\gamma})+\int_{-\pi/2}^{0}d\phi_{-}(-2\phi_{-})\Theta(\beta(r)-\cos\gamma).$
(A.4)

The fourth integral vanishes by symmetry. The sum of the second and fifth
integrals vanishes for the same reason. The third integral is symmetric around
zero, so that we can simplify $J_{c}$ to the following expression:

$\displaystyle
J_{c}=\int_{0}^{\pi/2}d\phi_{-}\int_{-1}^{1}d\cos\theta\frac{1}{(1-\cos\theta+\mu)^{2}}\left[\Theta(\beta(r)-\cos\gamma)+\Theta(\beta(r)-\cos\bar{\gamma})\right].$
(A.5)

We note that $J_{c}=J_{c}(\beta(r),\cos\alpha,\mu)$, with $-1\leq\beta(r)\leq
1$ and $\beta(r)\leq\cos\alpha\leq 1$. We can take advantage of this property
to reduce the dimension of the numerical integration of $C_{\text{c}}$. We
precompute $J_{c}$ on the following discretized grid: $\beta\in[-1,1]$ in 60
steps, $\cos\alpha\in[-1,1]$ in 70 steps and $\mu\in[10^{-15.5},10^{12.5}]$ in
28 log steps. The additional step function $\Theta(\cos\alpha-\beta(r))$ in
Eq. 2.24 is included in the subsequent numerical computation over the particle
speeds. Loosely speaking, this procedure reduces a 6-dimensional integral to a
4$\oplus$2-dimensional integral, significantly reducing the computational
cost.

##  Appendix B The angular integral $J_{sc}$

In this appendix we discuss the details of the analytic and numerical
calculation of the angular integral introduced in the calculation of the self-
capture coefficient. We begin with the calculation of $J_{sc}$, which was
defined in Eq. 3.20,

$\displaystyle
J_{sc}=\frac{1}{4\pi}\int_{0}^{2\pi}d\phi_{CM}\int_{0}^{2\pi}d\phi_{\chi}\int_{-1}^{1}d\cos\theta
f_{sc}(\cos\theta,\mu)\Theta(\beta(r)-|\cos\gamma|),$ (B.1)

with Eq. 3.4 defining $f_{sc}$. We note again that only the difference of the
azimuthal angles appears in the integrand through $\cos\gamma$. Therefore, we
can apply the same argument used in Appendix A, by switching to
$\phi_{\pm}=\phi_{CM}\pm\phi_{\chi}$:

$\displaystyle
J_{sc}=\frac{1}{4\pi}\int_{0}^{2\pi}d\phi_{-}(4\pi-2\phi_{-})\int_{-1}^{1}d\cos\theta
f_{sc}(\cos\theta,\mu)\Theta(\beta(r)-|\cos\gamma|).$ (B.2)

We continue by applying the same logic as in Eq. A.4, where we partition and
shift the $\phi_{-}$ integration region to obtain

$\displaystyle J_{sc}=\int_{0}^{\pi/2}d\phi_{-}\int_{-1}^{1}d\cos\theta
f_{sc}(\cos\theta,\mu)\left(\Theta(\beta(r)-|\cos\gamma|)+\Theta(\beta(r)-|\cos\bar{\gamma}|)\right),$
(B.3)

where we recall that
$\cos\bar{\gamma}\equiv-\sin\alpha\sin\theta\cos\phi_{-}+\cos\alpha\cos\theta$.
Next we define
$f_{sc}^{sym}(\cos\theta,\mu)=f_{sc}(\cos\theta,\mu)+f_{sc}(-\cos\theta,\mu)$
and reduce the domain of integration of $\cos\theta$ by using the symmetry
properties of the integrand

$\displaystyle J_{sc}=\int_{0}^{\pi/2}d\phi_{-}\int_{0}^{1}d\cos\theta
f_{sc}^{sym}(\cos\theta,\mu)\left(\Theta(\beta(r)-|\cos\gamma|)+\Theta(\beta(r)-|\cos\bar{\gamma}|)\right).$
(B.4)

Note that $J_{sc}=J_{sc}(\beta(r),\cos\alpha,\mu)$, with
$0\leq\beta(r)\leq\cos\alpha\leq 1$ (see Eq. 3.19). Similar to $J_{c}$, we
calculate $J_{sc}$ by discretizing it on the following grid: $\beta\in[0,1]$
in 60 steps, $\cos\alpha\in[0,1]$ in 70 steps and
$\mu\in[10^{-15.5},10^{12.5}]$ in 28 log steps.

##  Appendix C Dark Matter self-evaporation

In this appendix we discuss the calculation of the self-evaporation rate.
Self-evaporation occurs when two captured DM particles scatter such that one
of the particles becomes gravitationally unbound after the collision. We
stress that self-evaporation does not occur in the zero-temperature
approximation, since the captured DM particles have no kinetic energy to
overcome the potential well. Therefore, one is confronted with incorporating
thermal motion of the target particles from the outset.

The calculation proceeds in analogy with Secs. 2 and 3. We will compute the
total self-evaporation rate, accounting for particle ($\chi$) and antiparticle
($\bar{\chi}$) self-evaporation by scattering against both particles and
antiparticles,

$\displaystyle C_{\text{sevap}}=\int_{0}^{R}dr4\pi
r^{2}\left(\frac{dC^{\chi\chi}_{\text{sevap}}}{dV}+\frac{dC^{\chi\bar{\chi}}_{\text{sevap}}}{dV}+\frac{dC^{\bar{\chi}\chi}_{\text{sevap}}}{dV}+\frac{dC^{\bar{\chi}\bar{\chi}}_{\text{sevap}}}{dV}\right),$
(C.1)

where

$\displaystyle\frac{dC^{\chi\bar{\chi}}_{\text{sevap}}}{dV}=n^{\chi}_{c}(r)n^{\bar{\chi}}_{c}(r)\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{1}f_{c}(\vec{v}_{1})\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{2}f_{c}(\vec{v}_{2})|\vec{v}_{1}-\vec{v}_{2}|\sigma^{\chi\bar{\chi}}_{\text{sevap}}$
(C.2)

is the self-evaporation per unit volume of particles by scattering against
antiparticles, and the other three $dC_{\text{sevap}}/dV$ terms are defined
analogously. The colliding DM particles, with initial velocities $\vec{v}_{1}$
and $\vec{v}_{2}$, are described by the same Maxwell distribution at the core
temperature of the capturing body, given by Eq. 3.11. In the case where the
particle and antiparticle populations are equal,
$n^{\chi}_{c}(r)=n^{\bar{\chi}}_{c}(r)=1/2\,n_{c}(r)$, where $n_{c}(r)$ is the
total number density of captured DM (see Eq. 3.1), and we write

$\displaystyle C_{\text{sevap}}=\frac{1}{4}\int_{0}^{R}dr4\pi
r^{2}\frac{dC_{\text{sevap}}}{dV},$
$\displaystyle\frac{dC_{\text{sevap}}}{dV}=n^{2}_{c}(r)$
$\displaystyle\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{1}f_{c}(\vec{v}_{1})\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{2}f_{c}(\vec{v}_{2})|\vec{v}_{1}-\vec{v}_{2}|\sigma_{\text{sevap}}.$
(C.3)

We begin by transforming the kinematic integrals to the CM frame, where we use
the definitions introduced in Eq. 3.15,

$\displaystyle\int_{0}^{v_{\text{esc}}(r)}$ $\displaystyle
d^{3}\vec{v}_{1}f_{c}(\vec{v}_{1})\int_{0}^{v_{\text{esc}}(r)}d^{3}\vec{v}_{2}f_{c}(\vec{v}_{2})=$
$\displaystyle=\int_{0}^{\infty}$ $\displaystyle
d^{3}\vec{v}_{1}\Theta(v_{\text{esc}}(r)-|\vec{v}_{1}|)f_{c}(\vec{v}_{1})\int_{0}^{\infty}d^{3}\vec{v}_{2}\Theta(v_{\text{esc}}(r)-|\vec{v}_{2}|)f_{c}(\vec{v}_{2})=$
$\displaystyle=\int_{0}^{\infty}$ $\displaystyle
d^{3}\vec{v}_{CM}\int_{0}^{\infty}d^{3}\vec{v}^{\prime}_{1}Jf_{c}(\vec{v}_{CM}+\vec{v}^{\prime}_{1})f_{c}(\vec{v}_{CM}-\vec{v}^{\prime}_{1})\times$
$\displaystyle\times\Theta(v_{\text{esc}}(r)-|\vec{v}_{CM}+\vec{v}^{\prime}_{1}|)\Theta(v_{\text{esc}}(r)-|\vec{v}_{CM}-\vec{v}^{\prime}_{1}|),$
(C.4)

where $J=8$ is the Jacobian of the transformation. First, we note the
following simplification for the exponentials appearing in the product of the
two Maxwell distributions:

$\displaystyle
e^{-\frac{M_{\text{DM}}}{2T_{\text{Sun}}}(v^{2}_{CM}+v^{{}^{\prime}2}_{1}+2\vec{v}_{CM}\cdot\vec{v}^{\prime}_{1})}e^{-\frac{M_{\text{DM}}}{2T_{\text{Sun}}}(v^{2}_{CM}+v^{{}^{\prime}2}_{1}-2\vec{v}_{CM}\cdot\vec{v}^{\prime}_{1})}=e^{-\frac{M_{\text{DM}}}{2T_{\text{Sun}}}(2v^{2}_{CM}+2v^{{}^{\prime}2}_{1})}.$
(C.5)

As before, we orient the $z$-axis of our frame along $\vec{v}^{\prime}_{1}$
such that the solid angle integration of $d^{3}\vec{v}^{\prime}_{1}$ proceeds
trivially. We let $\alpha$ and $\phi_{CM}$ specify the polar and azimuthal
angles of $\vec{v}_{CM}$ with respect to the $z$-axis. With this convention,
the two step functions
$\Theta(v_{\text{esc}}(r)-|\vec{v}_{CM}\pm\vec{v}^{\prime}_{1}|)$ in the
kinematic integrals of Eq. C.4 turn into bounds on $c\alpha\equiv\cos\alpha$:

$\displaystyle v^{2}_{\text{esc}}(r)\geq v^{2}_{CM}+v^{{}^{\prime}2}_{1}\pm
2v_{CM}v^{\prime}_{1}c\alpha,$ (C.6)

which we write compactly as

$\displaystyle|c\alpha|\leq\beta(r),$ (C.7)

with $\beta$ defined in Eq. 2.23. Collecting these calculations together,

$\displaystyle\frac{dC_{\text{sevap}}}{dV}=n^{2}_{c}(r)\int_{0}^{\infty}$
$\displaystyle
v^{2}_{CM}dv_{CM}\int_{0}^{\infty}v^{{}^{\prime}2}_{1}dv^{\prime}_{1}\int_{-1}^{1}dc\alpha\Theta(\beta(r)-|c\alpha|)\times$
$\displaystyle\times\int_{0}^{2\pi}d\phi_{CM}J4\pi
f_{c}(\sqrt{2}v_{CM})f_{c}(\sqrt{2}v^{\prime}_{1})2v^{\prime}_{1}\sigma_{\text{sevap}}.$
(C.8)

Next we calculate the self-evaporation cross-section in the non-relativistic
limit. For $\chi\chi$ and $\bar{\chi}\bar{\chi}$ scattering, both $t$\- and
$u$-channel diagrams contribute. The rate for these processes includes a
symmetry factor of $1/2$ to account for the identical initial state particles,
which for simplicity we incorporate into an effective matrix element. For
$\bar{\chi}\chi$ and $\chi\bar{\chi}$ self-evaporation only the $t$-channel
diagram contributes. Combining all four processes, we can write the effective
matrix element-squared as

$\displaystyle\quad\quad\quad|\mathcal{M}|^{2}=\frac{4g_{D}^{4}}{v^{{}^{\prime}4}_{1}}f_{\text{sevap}}(\cos\theta,\mu),\text{
where}$ (C.9) $\displaystyle f_{\text{sevap}}(\cos\theta,\mu)=$
$\displaystyle\frac{3}{(1-\cos\theta+\mu)^{2}}+\frac{1}{(1+\cos\theta+\mu)^{2}}+\frac{1}{(1-\cos\theta+\mu)(1+\cos\theta+\mu)}.$
(C.10)

Self-evaporation occurs when one of the particles has a final speed larger
than the local escape speed. This kinematic constraint is

$\displaystyle\Theta_{\text{sevap}}=\Theta(|\vec{v}^{lab}_{1,f}|-v_{\text{esc}}(r))=\Theta(\cos\gamma-\beta(r)),$
(C.11)

where the last equality is written in terms of the angle $\gamma$, defined in
Eq. 2.28. We convolute this step function with the matrix element above and
integrate over the solid angle of the final state particles to obtain the
self-evaporation cross-section,

$\displaystyle\sigma_{\text{sevap}}=\frac{1}{256\pi^{2}M^{2}_{\text{DM}}}$
$\displaystyle\int d\cos\theta
d\phi_{\chi}|\mathcal{M}|^{2}\Theta_{\text{sevap}}=$
$\displaystyle=\frac{\alpha^{2}_{D}}{4M^{2}_{\text{DM}}v^{{}^{\prime}4}_{1}}\int
d\cos\theta
d\phi_{\chi}f_{\text{sevap}}(\cos\theta,\mu)\Theta(\cos\gamma-\beta(r)).$
(C.12)

Inserting Eq. C.12 in Eq. C.8 and collecting terms gives the following result
for self-evaporation per unit volume:

$\displaystyle\frac{dC_{\text{sevap}}}{dV}=n^{2}_{c}(r)\int_{0}^{\infty}v^{2}_{CM}dv_{CM}$
$\displaystyle\int_{0}^{\infty}v^{{}^{\prime}2}_{1}dv^{\prime}_{1}\Theta(\beta(r))\Theta(1-\beta(r))\times$
$\displaystyle\times(4\pi)^{2}Jf_{c}(\sqrt{2}v_{CM})f_{c}(\sqrt{2}v^{\prime}_{1})2v^{\prime}_{1}\frac{\alpha^{2}_{D}}{4M^{2}_{\text{DM}}v^{{}^{\prime}4}_{1}}J_{\text{sevap}}(\beta(r),\mu),$
(C.13)

where we introduced the angular integral

$\displaystyle
J_{\text{sevap}}(\beta(r),\mu)\equiv\frac{1}{4\pi}\int_{0}^{2\pi}d\phi_{CM}$
$\displaystyle\int_{0}^{2\pi}d\phi_{\chi}\int_{-1}^{1}dc\alpha\int_{-1}^{1}d\cos\theta\times$
$\displaystyle\times
f_{\text{sevap}}(\cos\theta,\mu)\Theta(\beta(r)-|c\alpha|)\Theta(\cos\gamma-\beta(r)).$
(C.14)

The two step functions in Eq. C.14 impose $0\leq\beta(r)\leq 1$ for non-zero
$J_{\text{sevap}}$. For clarity, we further explicitly include this condition
in Eq. C.13. In what follows we consider the computation of
$J_{\text{sevap}}$. We apply the transformation
$\phi_{\pm}=\phi_{CM}\pm\phi_{\chi}$, outlined in Eqs. A.2–A.3, to obtain

$\displaystyle
J_{\text{sevap}}(\beta(r),\mu)=\frac{1}{4\pi}\int_{-1}^{1}dc\alpha$
$\displaystyle\int_{0}^{2\pi}d\phi_{-}(4\pi-2\phi_{-})\int_{-1}^{1}d\cos\theta\times$
$\displaystyle\times\Theta(\beta(r)-|c\alpha|)f_{\text{sevap}}(\cos\theta,\mu)\Theta(\cos\gamma-\beta(r)),$
(C.15)

followed by the symmetry argument used in Eq. A.4,

$\displaystyle J_{\text{sevap}}(\beta(r),\mu)=\int_{0}^{\pi/2}$ $\displaystyle
d\phi_{-}\int_{-1}^{1}dc\alpha\int_{-1}^{1}d\cos\theta\Theta(\beta(r)-|c\alpha|)\times$
$\displaystyle\times
f_{\text{sevap}}(\cos\theta,\mu)\big{(}\Theta(\cos\gamma-\beta(r))+\Theta(\cos\bar{\gamma}-\beta(r))\big{)},$
(C.16)

where we recall that
$\cos\bar{\gamma}\equiv-s\alpha\sin\theta\cos\phi_{-}+c\alpha\cos\theta$. We
note, however, that $\Theta(\cos\bar{\gamma}-\beta(r))$ is always zero in the
given volume of integration. To see this, note that
$\Theta(\beta(r)-|c\alpha|)$ imposes $c\alpha\cos\theta\leq\beta(r)$, where we
also recall the kinematic constraint $0\leq\beta\leq 1$. Therefore,
$\cos\bar{\gamma}=-s\alpha\sin\theta\cos\phi_{-}+c\alpha\cos\theta\leq\beta$
since the first term is non-positive in the integration region. We write the
$\Theta(\cos\gamma-\beta(r))$ step function as a quadratic inequality for
$c\alpha$:

$\displaystyle 0>($
$\displaystyle\cos^{2}\theta+\sin^{2}\theta\cos^{2}\phi_{-})c^{2}\alpha-2\beta(r)\cos\theta
c\alpha+\beta^{2}(r)-\sin^{2}\theta\cos^{2}\phi_{-},$ (C.17)

which is satisfied for $c\alpha\in(c\alpha_{-},c\alpha_{+})$, where

$\displaystyle
c\alpha_{\pm}=\frac{\beta(r)\cos\theta\pm\sqrt{\sin^{2}\theta\cos^{2}\phi_{-}(\cos^{2}\theta+\sin^{2}\theta\cos^{2}\phi_{-}-\beta^{2}(r))}}{\cos^{2}\theta+\sin^{2}\theta\cos^{2}\phi_{-}}.$
(C.18)

Figure 5: Left: Self-evaporation coefficient in the Sun as a function of
mediator mass $M_{\gamma_{D}}$. Dashed lines are included to highlight the
parametric scaling with $M_{\gamma_{D}}$ in the long- and short-range regimes.
Right: Self-evaporation coefficient in the Sun as a function of DM mass. We
fix $M_{\gamma_{D}}=10^{-6}$ GeV, corresponding to long-range scattering, and
$\alpha_{D}=10^{-3}$.

To ensure that $c\alpha_{\pm}$ are real we impose the positivity of the
discriminant with
$\Theta(\cos^{2}\theta+\sin^{2}\theta\cos^{2}\phi_{-}-\beta^{2}(r))$. The
$c\alpha$ integral can now be evaluated in closed form, producing four terms
which correspond to different orderings of upper and lower limits of
integration on $c\alpha$ imposed by the $\Theta(\beta(r)-|c\alpha|)$ and
$\Theta(c\alpha_{+}-c\alpha)\Theta(c\alpha-c\alpha_{-})$ step functions.
Inserting this calculation in Eq. C.16 gives

$\displaystyle J_{\text{sevap}}$
$\displaystyle(\beta(r),\mu)=\int_{0}^{\pi/2}d\phi_{-}\int_{-1}^{1}d\cos\theta
f_{\text{sevap}}(\cos\theta,\mu)\Theta(\cos^{2}\theta+\sin^{2}\theta\cos^{2}\phi_{-}-\beta^{2}(r))\times$
$\displaystyle\times\big{[}$
$\displaystyle(c\alpha_{+}-c\alpha_{-})\Theta(\beta(r)-c\alpha_{+})\Theta(c\alpha_{-}+\beta(r))+2\beta(r)\Theta(c\alpha_{+}-\beta(r))\Theta(-\beta(r)-c\alpha_{-})+$
$\displaystyle+(c\alpha_{+}+\beta(r))\Theta(\beta(r)-c\alpha_{+})\Theta(-\beta(r)-c\alpha_{-})\Theta(c\alpha_{+}+\beta(r))+$
$\displaystyle+(\beta(r)-c\alpha_{-})\Theta(c\alpha_{+}-\beta(r))\Theta(c\alpha_{-}+\beta(r))\Theta(\beta(r)-c\alpha_{-})\big{]}.$
(C.19)

We take advantage of the restricted functional dependence
$J_{\text{sevap}}(M_{\text{DM}},M_{\gamma_{D}},r,v_{CM},v^{\prime}_{1})=J_{\text{sevap}}(\beta(r),\mu)$,
with $\beta(r)$ in the finite range $0\leq\beta(r)\leq 1$, to precompute
$J_{\text{sevap}}$ numerically by placing it on the following grid:
$\beta\in[0,1]$ in 70 even steps and $\mu\in[10^{-15.5},10^{12.5}]$ in 28 log
steps. We collect all calculations to write the final expression for the self-
evaporation coefficient

$\displaystyle
C_{\text{sevap}}=(4\pi)^{3}\frac{\alpha^{2}_{D}}{M^{2}_{\text{DM}}}$
$\displaystyle\int_{0}^{R}drr^{2}n^{2}_{c}(r)\int_{0}^{v_{\text{esc}}(r)}dv_{CM}v^{2}_{CM}\times$
$\displaystyle\times\int_{v_{\text{esc}}(r)-v_{CM}}^{\sqrt{v^{2}_{\text{esc}}(r)-v^{2}_{CM}}}\frac{dv^{\prime}_{1}}{v^{\prime}_{1}}f_{c}(\sqrt{2}v_{CM})f_{c}(\sqrt{2}v^{\prime}_{1})J_{\text{sevap}}(\beta(r),\mu).$
(C.20)

We show numerical calculations of $C_{\text{sevap}}$ for the Sun in Fig. 5. We
note the strong exponential suppression of $C_{\text{sevap}}$ with increasing
$M_{\text{DM}}$, which can be traced to the Maxwell distribution for the
captured DM particles, given by Eq. 3.11. Similar to capture and self-capture,
self-evaporation features two distinct scattering regimes. The short-range
regime occurs for large $M_{\gamma_{D}}$, when the typical momentum transfer
is much smaller than the mediator mass. The interactions are effectively
pointlike and result in a $C_{\text{sevap}}\sim\mu^{-2}\sim
M^{-4}_{\gamma_{D}}$ scaling. Long-range scattering dominates in the opposite
regime, when $M_{\gamma_{D}}$ is small compared to the typical momentum
transfer, and result in a $C_{\text{sevap}}\sim\mu^{-1/2}\sim
M^{-1}_{\gamma_{D}}$ scaling with mediator mass.

## References

  * [1] J. Silk, K. A. Olive and M. Srednicki, _The Photino, the Sun and High-Energy Neutrinos_ , _Phys. Rev. Lett._ 55 (1985) 257–259.
  * [2] P. Schuster, N. Toro, N. Weiner and I. Yavin, _High Energy Electron Signals from Dark Matter Annihilation in the Sun_ , _Phys. Rev. D_ 82 (2010) 115012, [0910.1839].
  * [3] B. Batell, M. Pospelov, A. Ritz and Y. Shang, _Solar Gamma Rays Powered by Secluded Dark Matter_ , _Phys. Rev._ D81 (2010) 075004, [0910.1567].
  * [4] L. M. Krauss, M. Srednicki and F. Wilczek, _Solar System Constraints and Signatures for Dark Matter Candidates_ , _Phys. Rev. D_ 33 (1986) 2079–2083.
  * [5] C. Kouvaris, _WIMP Annihilation and Cooling of Neutron Stars_ , _Phys. Rev. D_ 77 (2008) 023006, [0708.2362].
  * [6] G. D. Mack, J. F. Beacom and G. Bertone, _Towards Closing the Window on Strongly Interacting Dark Matter: Far-Reaching Constraints from Earth’s Heat Flow_ , _Phys. Rev. D_ 76 (2007) 043523, [0705.4298].
  * [7] G. Bertone and M. Fairbairn, _Compact Stars as Dark Matter Probes_ , _Phys. Rev. D_ 77 (2008) 043515, [0709.1485].
  * [8] M. Baryakhtar, J. Bramante, S. W. Li, T. Linden and N. Raj, _Dark Kinetic Heating of Neutron Stars and An Infrared Window On WIMPs, SIMPs, and Pure Higgsinos_ , _Phys. Rev. Lett._ 119 (2017) 131801, [1704.01577].
  * [9] R. K. Leane and J. Smirnov, _Exoplanets as Sub-GeV Dark Matter Detectors_ , _Phys. Rev. Lett._ 126 (2021) 161101, [2010.00015].
  * [10] D. N. Spergel and W. H. Press, _Effect of hypothetical, weakly interacting, massive particles on energy transport in the solar interior_ , _Astrophys. J._ 294 (1985) 663–673.
  * [11] I. Goldman and S. Nussinov, _Weakly Interacting Massive Particles and Neutron Stars_ , _Phys. Rev. D_ 40 (1989) 3221–3230.
  * [12] I. P. Lopes, G. Bertone and J. Silk, _Solar seismic model as a new constraint on supersymmetric dark matter_ , _Mon. Not. Roy. Astron. Soc._ 337 (2002) 1179–1184, [astro-ph/0205066].
  * [13] A. de Lavallaz and M. Fairbairn, _Neutron Stars as Dark Matter Probes_ , _Phys. Rev. D_ 81 (2010) 123521, [1004.0629].
  * [14] M. T. Frandsen and S. Sarkar, _Asymmetric dark matter and the Sun_ , _Phys. Rev. Lett._ 105 (2010) 011301, [1003.4505].
  * [15] D. T. Cumberbatch, J. A. Guzik, J. Silk, L. S. Watson and S. M. West, _Light WIMPs in the Sun: Constraints from Helioseismology_ , _Phys. Rev. D_ 82 (2010) 103503, [1005.5102].
  * [16] M. Taoso, F. Iocco, G. Meynet, G. Bertone and P. Eggenberger, _Effect of low mass dark matter particles on the Sun_ , _Phys. Rev. D_ 82 (2010) 083509, [1005.5711].
  * [17] S. D. McDermott, H.-B. Yu and K. M. Zurek, _Constraints on Scalar Asymmetric Dark Matter from Black Hole Formation in Neutron Stars_ , _Phys. Rev. D_ 85 (2012) 023519, [1103.5472].
  * [18] A. Gould, _Resonant Enhancements in WIMP Capture by the Earth_ , _Astrophys.J._ 321 (1987) 571.
  * [19] G. Busoni, A. De Simone and W.-C. Huang, _On the Minimum Dark Matter Mass Testable by Neutrinos from the Sun_ , _JCAP_ 1307 (2013) 010, [1305.1817].
  * [20] J. Fan, A. Katz and J. Shelton, _Direct and indirect detection of dissipative dark matter_ , _JCAP_ 1406 (2014) 059, [1312.1336].
  * [21] B. Dasgupta, A. Gupta and A. Ray, _Dark matter capture in celestial objects: light mediators, self-interactions, and complementarity with direct detection_ , _JCAP_ 10 (2020) 023, [2006.10773].
  * [22] J. L. Feng, M. Kaplinghat, H. Tu and H.-B. Yu, _Hidden Charged Dark Matter_ , _JCAP_ 07 (2009) 004, [0905.3039].
  * [23] J. L. Feng, M. Kaplinghat and H.-B. Yu, _Halo Shape and Relic Density Exclusions of Sommerfeld-Enhanced Dark Matter Explanations of Cosmic Ray Excesses_ , _Phys. Rev. Lett._ 104 (2010) 151301, [0911.0422].
  * [24] A. Loeb and N. Weiner, _Cores in Dwarf Galaxies from Dark Matter with a Yukawa Potential_ , _Phys. Rev. Lett._ 106 (2011) 171302, [1011.6374].
  * [25] F.-Y. Cyr-Racine and K. Sigurdson, _Cosmology of atomic dark matter_ , _Phys. Rev. D_ 87 (2013) 103515, [1209.5752].
  * [26] J. Fan, A. Katz, L. Randall and M. Reece, _Double-Disk Dark Matter_ , _Phys. Dark Univ._ 2 (2013) 139–156, [1303.1521].
  * [27] D. Curtin and J. Setford, _Signatures of Mirror Stars_ , _JHEP_ 03 (2020) 041, [1909.04072].
  * [28] P. Schuster, N. Toro and I. Yavin, _Terrestrial and Solar Limits on Long-Lived Particles in a Dark Sector_ , _Phys. Rev._ D81 (2010) 016002, [0910.1602].
  * [29] P. Meade, S. Nussinov, M. Papucci and T. Volansky, _Searches for Long Lived Neutral Particles_ , _JHEP_ 06 (2010) 029, [0910.4160].
  * [30] N. F. Bell and K. Petraki, _Enhanced neutrino signals from dark matter annihilation in the Sun via metastable mediators_ , _JCAP_ 1104 (2011) 003, [1102.2958].
  * [31] The Fermi LAT collaboration, M. Ajello et al., _Constraints on dark matter models from a Fermi LAT search for high-energy cosmic-ray electrons from the Sun_ , _Phys. Rev._ D84 (2011) 032007, [1107.4272].
  * [32] J. L. Feng, J. Smolinsky and P. Tanedo, _Dark Photons from the Center of the Earth: Smoking-Gun Signals of Dark Matter_ , _Phys. Rev._ D93 (2016) 015014, [1509.07525].
  * [33] J. L. Feng, J. Smolinsky and P. Tanedo, _Detecting dark matter through dark photons from the Sun: Charged particle signatures_ , _Phys. Rev._ D93 (2016) 115036, [1602.01465].
  * [34] ANTARES collaboration, S. Adrián-Martínez et al., _A search for Secluded Dark Matter in the Sun with the ANTARES neutrino telescope_ , _JCAP_ 1605 (2016) 016, [1602.07000].
  * [35] R. K. Leane, K. C. Y. Ng and J. F. Beacom, _Powerful Solar Signatures of Long-Lived Dark Mediators_ , _Phys. Rev. D_ 95 (2017) 123016, [1703.04629].
  * [36] C. Arina, M. Backović, J. Heisig and M. Lucente, _Solar $\gamma$ rays as a complementary probe of dark matter_, _Phys. Rev._ D96 (2017) 063010, [1703.08087].
  * [37] M. Ardid, I. Felis, A. Herrero and J. A. Martínez-Mora, _Constraining Secluded Dark Matter models with the public data from the 79-string IceCube search for dark matter in the Sun_ , _JCAP_ 04 (2017) 010, [1701.08863].
  * [38] C. Niblaeus, A. Beniwal and J. Edsjo, _Neutrinos and gamma rays from long-lived mediator decays in the Sun_ , _JCAP_ 11 (2019) 011, [1903.11363].
  * [39] N. F. Bell, J. B. Dent and I. W. Sanderson, _Solar Gamma Ray Constraints on Dark Matter Annihilation to Secluded Mediators_ , 2103.16794.
  * [40] R. Garani and S. Palomares-Ruiz, _Dark matter in the Sun: scattering off electrons vs nucleons_ , _JCAP_ 05 (2017) 007, [1702.02768].
  * [41] J. Bovy, D. W. Hogg and H.-W. Rix, _Galactic masers and the milky way circular velocity_ , _The Astrophysical Journal_ 704 (2009) 1704.
  * [42] M. Reid, K. Menten, X. Zheng, A. Brunthaler, L. Moscadelli, Y. Xu et al., _Trigonometric parallaxes of massive star-forming regions. vi. galactic structure, fundamental parameters, and noncircular motions_ , _The Astrophysical Journal_ 700 (2009) 137.
  * [43] P. J. McMillan and J. J. Binney, _The uncertainty in galactic parameters_ , _Monthly Notices of the Royal Astronomical Society_ 402 (2010) 934–940.
  * [44] R. Schönrich, J. Binney and W. Dehnen, _Local kinematics and the local standard of rest_ , _Monthly Notices of the Royal Astronomical Society_ 403 (2010) 1829–1833.
  * [45] J. I. Read, _The Local Dark Matter Density_ , _J. Phys. G_ 41 (2014) 063101, [1404.1938].
  * [46] J. N. Bahcall, A. M. Serenelli and S. Basu, _New solar opacities, abundances, helioseismology, and neutrino fluxes_ , _Astrophys.J._ 621 (2005) L85–L88, [astro-ph/0412440].
  * [47] B. Holdom, _Two U(1)’s and Epsilon Charge Shifts_ , _Phys. Lett. B_ 166 (1986) 196–198.
  * [48] A. Gould, _WIMP Distribution in and Evaporation From the Sun_ , _Astrophys.J._ 321 (1987) 560.
  * [49] A. Gould, _Direct and indirect capture of weakly interacting massive particles by the earth_ , _The Astrophysical Journal_ 328 (1988) 919–939.
  * [50] J. Chen, Z.-L. Liang, Y.-L. Wu and Y.-F. Zhou, _Long-range self-interacting dark matter in the Sun_ , _JCAP_ 1512 (2015) 021, [1505.04031].
  * [51] G. Busoni, A. De Simone, P. Scott and A. C. Vincent, _Evaporation and scattering of momentum- and velocity-dependent dark matter in the Sun_ , _JCAP_ 10 (2017) 037, [1703.07784].
  * [52] P. Baratella, M. Cirelli, A. Hektor, J. Pata, M. Piibeleht and A. Strumia, _PPPC 4 DM $\nu$: a Poor Particle Physicist Cookbook for Neutrinos from Dark Matter annihilations in the Sun_, _JCAP_ 03 (2014) 053, [1312.6408].
  * [53] C. Gaidau and J. Shelton, _A Solar System Test of Self-Interacting Dark Matter_ , _JCAP_ 06 (2019) 022, [1811.00557].
  * [54] A. R. Zentner, _High-Energy Neutrinos From Dark Matter Particle Self-Capture Within the Sun_ , _Phys.Rev._ D80 (2009) 063501, [0907.3448].
  * [55] HAWC collaboration, A. Albert et al., _Constraints on Spin-Dependent Dark Matter Scattering with Long-Lived Mediators from TeV Observations of the Sun with HAWC_ , _Phys. Rev. D_ 98 (2018) 123012, [1808.05624].
  * [56] J. Redondo and G. Raffelt, _Solar constraints on hidden photons re-visited_ , _JCAP_ 08 (2013) 034, [1305.2920].
  * [57] H. An, M. Pospelov and J. Pradler, _New stellar constraints on dark photons_ , _Phys. Lett. B_ 725 (2013) 190–195, [1302.3884].
  * [58] H. An, M. Pospelov, J. Pradler and A. Ritz, _Direct Detection Constraints on Dark Photon Dark Matter_ , _Phys. Lett. B_ 747 (2015) 331–338, [1412.8378].
  * [59] J. H. Chang, R. Essig and S. D. McDermott, _Revisiting Supernova 1987A Constraints on Dark Photons_ , _JHEP_ 01 (2017) 107, [1611.03864].
  * [60] J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu et al., _Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump_ , _Phys. Rev._ D38 (1988) 3375.
  * [61] S. Andreas, C. Niebuhr and A. Ringwald, _New Limits on Hidden Photons from Past Electron Beam Dumps_ , _Phys. Rev._ D86 (2012) 095019, [1209.6083].
  * [62] LSND collaboration, C. Athanassopoulos et al., _Evidence for muon-neutrino $\rightarrow$ electron-neutrino oscillations from pion decay in flight neutrinos_, _Phys. Rev._ C58 (1998) 2489–2511, [nucl-ex/9706006].
  * [63] R. Essig, R. Harnik, J. Kaplan and N. Toro, _Discovering New Light States at Neutrino Experiments_ , _Phys. Rev._ D82 (2010) 113008, [1008.0636].
  * [64] J. Redondo and M. Postma, _Massive hidden photons as lukewarm dark matter_ , _JCAP_ 02 (2009) 005, [0811.0326].
  * [65] M. Pospelov, A. Ritz and M. B. Voloshin, _Bosonic super-WIMPs as keV-scale dark matter_ , _Phys. Rev. D_ 78 (2008) 115012, [0807.3279].
  * [66] W. DeRocco, P. W. Graham, D. Kasen, G. Marques-Tavares and S. Rajendran, _Observable signatures of dark photons from supernovae_ , _JHEP_ 02 (2019) 171, [1901.08596].
  * [67] E. Hardy and R. Lasenby, _Stellar cooling bounds on new light particles: plasma mixing effects_ , _JHEP_ 02 (2017) 033, [1611.05852].
  * [68] S. Knapen, T. Lin and K. M. Zurek, _Light Dark Matter: Models and Constraints_ , _Phys. Rev._ D96 (2017) 115021, [1709.07882].
  * [69] V. Bromm, _Formation of the First Stars_ , _Rept. Prog. Phys._ 76 (2013) 112901, [1305.5178].
  * [70] R. S. Klessen, _Formation of the first stars_ , 1807.06248.
  * [71] K. Freese, T. Rindler-Daller, D. Spolyar and M. Valluri, _Dark Stars: A Review_ , _Rept. Prog. Phys._ 79 (2016) 066902, [1501.02394].
  * [72] C. Ilie, C. Levy, J. Pilawa and S. Zhang, _Constraining Dark Matter properties with the first generation of stars_ , 2009.11474.
  * [73] H. An, M. Pospelov and J. Pradler, _Dark Matter Detectors as Dark Photon Helioscopes_ , _Phys. Rev. Lett._ 111 (2013) 041302, [1304.3461].
  * [74] A. Fradette, M. Pospelov, J. Pradler and A. Ritz, _Cosmological Constraints on Very Dark Photons_ , _Phys. Rev. D_ 90 (2014) 035022, [1407.0993].
  * [75] J. Jaeckel, J. Redondo and A. Ringwald, _Signatures of a hidden cosmic microwave background_ , _Phys. Rev. Lett._ 101 (2008) 131801, [0804.4157].
  * [76] D. F. Bartlett, P. E. Goldhagen and E. A. Phillips, _Experimental Test of Coulomb’s Law_ , _Phys. Rev. D_ 2 (1970) 483–487.
  * [77] D. F. Bartlett and S. Loegl, _LIMITS ON AN ELECTROMAGNETIC FIFTH FORCE_ , _Phys. Rev. Lett._ 61 (1988) 2285–2287.
  * [78] M. A. Buen-Abad, R. Essig, D. McKeen and Y.-M. Zhong, _Cosmological Constraints on Dark Matter Interactions with Ordinary Matter_ , 2107.12377.
  * [79] Z. Chacko, N. Craig, P. J. Fox and R. Harnik, _Cosmology in Mirror Twin Higgs and Neutrino Masses_ , _JHEP_ 07 (2017) 023, [1611.07975].
  * [80] N. Craig, S. Koren and T. Trott, _Cosmological Signals of a Mirror Twin Higgs_ , _JHEP_ 05 (2017) 038, [1611.07977].
  * [81] H. Vogel and J. Redondo, _Dark Radiation constraints on minicharged particles in models with a hidden photon_ , _JCAP_ 02 (2014) 029, [1311.2600].
  * [82] Z. Chacko, D. Curtin, M. Geller and Y. Tsai, _Cosmological Signatures of a Mirror Twin Higgs_ , _JHEP_ 09 (2018) 163, [1803.03263].