# Dark matter bound-state formation in the Sun

Xiaoyong Chu<EMAIL_ADDRESS>Institute of High Energy Physics,
Austrian Academy of Sciences, Nikolsdorfer Gasse 18, 1050 Vienna, Austria
Raghuveer Garani<EMAIL_ADDRESS>INFN Sezione di Firenze, Via G. Sansone 1,
I-50019 Sesto Fiorentino, Italy  Camilo García-Cely<EMAIL_ADDRESS>Instituto de Física Corpuscular (IFIC), Universitat de València-CSIC, Parc
Científic UV, C/ Catedrático José Beltrán 2, E-46980 Paterna, Spain  Thomas
Hambye<EMAIL_ADDRESS>Service de Physique Théorique Université Libre de
Bruxelles, Boulevard du Triomphe, CP225, 1050 Brussels, Belgium  CERN,
Theoretical Physics Department, Geneva, Switzerland

###### Abstract

The Sun may capture asymmetric dark matter (DM), which can subsequently form
bound-states through the radiative emission of a sub-GeV scalar. This process
enables generation of scalars without requiring DM annihilation. In addition
to DM capture on nucleons, the DM-scalar coupling responsible for bound-state
formation also induces capture from self-scatterings of ambient DM particles
with DM particles already captured, as well as with DM bound-states formed in-
situ within the Sun. This scenario is studied in detail by solving Boltzmann
equations numerically and analytically. In particular, we take into
consideration that the DM self-capture rates require a treatment beyond the
conventional Born approximation. We show that, thanks to DM scatterings on
bound-states, the number of DM particles captured increases exponentially,
leading to enhanced emission of relativistic scalars through bound-state
formation, whose final decay products could be observable. We explore
phenomenological signatures with the example that the scalar mediator decays
to neutrinos. We find that the neutrino flux emitted can be comparable to
atmospheric neutrino fluxes within the range of energies below one hundred
MeV. Future facilities like Hyper-K, and direct DM detection experiments can
further test such scenario.

††preprint: ULB-TH/24-01

###### Contents

  1. I Introduction
  2. II Number evolution of DM particles
    1. II.1 Determination of the various rates
    2. II.2 Integrating the Boltzmann equations
  3. III Neutrino flux and Terrestrial detection
  4. IV Results and constraints for an explicit model
    1. IV.1 Model
    2. IV.2 Constraints and results
  5. V Summary
  6. A Approximate parametric solution to number evolution of DM particles
    1. A.1 Case with no self-capture on DM bound-states
    2. A.2 Case with self-capture on DM bound-states
  7. B (Non-)perturbative treatment of elastic DM scattering
    1. B.1 Integrated DM self-scattering in DM halos
    2. B.2 Differential DM-DM(BS) scattering for DM capture
  8. C DM thermalization

##  I Introduction

Wherever dark matter (DM) particles are numerous, it is possible that two or
more of them form bound-states March-Russell and West (2009); Pospelov and
Ritz (2009); Shepherd  _et al._ (2009); Laha and Braaten (2014); von Harling
and Petraki (2014); Petraki  _et al._ (2014, 2015); Pearce  _et al._ (2015);
Laha (2015); An  _et al._ (2016a, b, 2017); Bi  _et al._ (2017); Kouvaris  _et
al._ (2016); Cirelli  _et al._ (2017); Asadi  _et al._ (2017); Petraki  _et
al._ (2017); Johnson  _et al._ (2016); Mitridate  _et al._ (2017); Baldes  _et
al._ (2018). This could occur if DM undergoes attractive self-interactions
mediated by a scalar or vector boson. In this case bound-state formation (BSF)
can occur radiatively via the emission of this mediator, leading to an
observable flux of particles. Indirect signals of BSF in the center of the
Milky Way has been studied for DM with no particle-antiparticle asymmetry, as
well as for asymmetric DM Pearce and Kusenko (2013); Baldes  _et al._ (2018);
Mahbubani  _et al._ (2020, 2021). In this work we consider the possibility
that the BSF process occurs in the Sun and study the corresponding DM indirect
detection signals.

Symmetric DM accumulating in the Sun from its capture on nucleons and their
corresponding indirect detection signals due to emitted meta-stable mediators
have been studied in numerous works, see e.g. Chen  _et al._ (2015); Feng  _et
al._ (2016); Arina  _et al._ (2017); Leane  _et al._ (2017); Lucente  _et al._
(2017); Niblaeus  _et al._ (2019); Cuoco  _et al._ (2020); Mazziotta  _et al._
(2020); Dasgupta  _et al._ (2020); Bell  _et al._ (2021a); Garani and
Palomares-Ruiz (2022); Abud  _et al._ (2021); Leane  _et al._ (2021); Abbasi
_et al._ (2022); Bell  _et al._ (2021b); Zakeri and Zhou (2022); Bose  _et
al._ (2022); Maity  _et al._ (2023). In these scenarios, the accumulation of
DM particles in the Sun reaches a saturation point when the annihilation rate
matches the capture rate. This is not the case of asymmetric DM scenarios as
DM cannot annihilate. Interestingly, this absence of annihilation permits a
greater accumulation of DM particles. However, it comes with the drawback of
not generating any indirect detection signal.

The BSF of asymmetric DM particles in the Sun from radiative emission of a
light scalar has the interesting property of allowing both DM indirect
detection and large accumulation of DM particles in the Sun. A crucial aspect
of this scenario is that the DM-scalar interaction that is needed for BSF
inherently implies that DM capture results not only from interactions with
nucleons but also from interactions with previously accumulated DM particles
and DM bound-states (DMBS). This means that the capture rate is larger than
that of the usual symmetric or asymmetric scenarios in which the capture only
arises from DM-nucleon scatterings. This allows for a larger accumulation and
an enhanced flux of emitted particles. To quantify this effect, it will be
necessary to compute the rates for DM-DM and DM-DMBS scatterings, which –as we
will show– receive non-perturbative contributions which can be calculated in
the semi-classical approximation. To our knowledge this possibility that
asymmetric DM particles form bound-states in the Sun and that DM particles are
captured by scattering off DM bound-states has not been considered
before.111Particle-antiparticle BSF inside the Sun for symmetric DM was
considered in Kouvaris  _et al._ (2016).

For concreteness, in this work we assume that the associated light scalars,
once emitted when the bound-states form, decay into SM particles, in
particular to neutrinos. As is well known, unlike other SM particles, low
energy neutrinos can escape the Sun leading to observable signatures even if
the decay takes place inside the Sun. After solving the set of Boltzmann
equations that describe the dynamics of DM accumulation in the Sun in section
II, we show in section III that an observable neutrino flux could arise from
the decays of the emitted mediator via the BSF of those accumulated DM.
Finally, in section IV we discuss relevant constraints from both astrophysical
observations and terrestrial experiments, such as DM self-interactions, BBN,
CMB, and direct/indirect searches.

##  II Number evolution of DM particles

Before considering quantitatively any concrete model, this section introduces
the basic relevant processes and related Boltzmann equations determining the
number of DM particles captured in the Sun. The framework of our interest is
based on an interaction between the DM particle and a lighter particle. For
definiteness in the following we consider that DM is made of Dirac fermions
$\chi$, with a Yukawa coupling $g_{s}$ to a lighter scalar particle $\phi$

${\cal L}\supset-g_{s}\phi\bar{\chi}\chi\,.$ (1)

To assume a scalar mediator is convenient because it induces an attractive
interaction.222 Such interactions do not result in the formation of mini black
hole in the Sun, i.e. the Chandrasekhar limit for fermions is not modified
Gresham and Zurek (2019); Garani  _et al._ (2022). We also assume that DM is
asymmetric, thus it does not annihilate in the present epoch. Nevertheless, in
addition to the usual capture due to scattering on target nuclei ($C_{\star}$)
Press and Spergel (1985); Griest and Seckel (1987), three additional terms
come into play in the Boltzmann equation determining the DM particle number
$N_{\chi}$, due to DM-$\phi$ interactions. First, depending on the strength of
self-interaction, DM particles could efficiently form $\chi$-$\chi$ bound-
states by emitting the particle mediating DM self-interactions. The rate of
the latter is denoted by $A_{\rm bsf}$. Since the number of free particles is
reduced by two per process, the term $-2\cdot\frac{1}{2}A_{\rm
bsf}N^{2}_{\chi}$ is introduced in the equation that describes the evolution
of $N_{\chi}$, where the second factor, $\frac{1}{2}$, counts for double-
counting of identical initial states. In contrast to the case of annihilation,
in this scenario DM particles are not lost by the system. In effect, there is
a build-up of DM bound-states within the celestial body, whose number
$N_{2\chi}$ is determined by a second Boltzmann equation that simply involves
a $+\frac{1}{2}A_{\rm bsf}N^{2}_{\chi}$ term. Second, there is a term coming
from capture of Galactic DM on already-captured DM particles, whose rate is
denoted by $C_{\chi}$. Finally, there is a term from capture of Galactic DM on
the formed bound-states, denoted by $C_{2\chi}$. The presence of both
$C_{\chi}$ and $C_{2\chi}$ terms increase DM accretion. In summary, both
populations evolve according to the following set of differential equations,

$\displaystyle\frac{{\rm d}N_{\chi}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle C_{\star}-A_{\rm
bsf}N^{2}_{\chi}+C_{\chi}N_{\chi}+C_{2\chi}N_{2\chi}\,,$ (2)
$\displaystyle\frac{{\rm d}N_{2\chi}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle\frac{1}{2}A_{\rm bsf}N^{2}_{\chi}~{}.$ (3)

The various rates $C_{\star}$, $C_{\chi}$ and $C_{2\chi}$ take into account
the fact that the three kinds of DM capture occur within different spheres
around the solar center. Here we have neglected self-capture via direct BSF
between a galactic DM and a captured DM particle, as such inelastic scattering
is typically much weaker than the elastic scattering with two DM initial
states. Also, in these Boltzmann equations we do not write down explicitly
additional terms coming from possible evaporation of the captured DM
particles. We have checked that the DM mass above which evaporation is
negligible is not significantly different than in the standard scenario (no-
self interaction effects). The value is found to be $m_{\rm evap}\approx 5$
GeV, which is about 50% larger than that in the standard scenario Garani and
Palomares-Ruiz (2022). Therefore, in the following phenomenological discussion
we consider only DM masses above $5$ GeV. Once DM particles are captured they
thermalize with the solar material, leading to efficient formation of bound-
states. We always assume that the thermalization happens quickly, and refer to
Appendix C for detailed calculation of the thermalization process.

The number of DM particles captured in any of the three ways cannot be larger
than the corresponding geometric rates, as the latter assume that all DM
particles crossing the corresponding thermal spheres are captured. This is
taken into account through the following conservative matching for the regimes
of small optical depth (thin) to large (thick)333These equations are accurate
if only one geometric rate can be saturated. This will be the case studied
here, for which the geometric capture rates on free DM ($C^{g}_{\chi}$) and on
nucleons ($C^{g}_{\star}$) are never reached.

$\displaystyle C_{\chi}N_{\chi}$ $\displaystyle\rightarrow$ $\displaystyle
C_{\chi}N_{\chi}\Theta(C^{g}_{\chi}-C_{\chi}N_{\chi})+C^{g}_{\chi}\Theta(C_{\chi}N_{\chi}-C^{g}_{\chi})~{},$
(4) $\displaystyle C_{2\chi}N_{2\chi}$ $\displaystyle\rightarrow$
$\displaystyle
C_{2\chi}N_{2\chi}\Theta(C^{g}_{2\chi}-C_{2\chi}N_{2\chi})+C^{g}_{2\chi}\Theta(C_{2\chi}N_{2\chi}-C^{g}_{2\chi})~{}.$
(5)

where the respective geometric rates on the DM and DMBS thermal spheres are
denoted by $C^{g}_{\chi}$ and $C^{g}_{2\chi}$. These geometric rates determine
the maximal possible capture rates, independently of the underlying particle
model.

Note that we do not consider bound-states containing more than two $\chi$
particles, under the assumption that there exist bottlenecks to form heavier
bound-states, such as fast decays of $(3\chi)\to(2\chi)+\chi$ and
$(4\chi)\to(2\chi)+(2\chi)$. See e.g. Wise and Zhang (2015); Hardy  _et al._
(2015); Gresham  _et al._ (2018, 2017) for further discussions. In practice,
if there is formation of many-body bound-states, each captured DM particle may
cause the emission of a few more mediator particles while thermal radii for
heavier bound-states shrink by a factor of few. The total effect at most
modifies our results mildly.444The presence of stable many-body bound-states,
as well as self-capture via inelastic scattering, may reduce DM evaporation
efficiency by quickly capturing free $\chi$ particles into heavier bound-
states. This can have qualitative consequences for DM candidates below GeV or
those captured by Earth, which is left for a future study.

We begin in the next subsection by determining the various rates. The reader
seeking immediate understanding of their interplay in the Boltzmann equations
can directly refer to section II.2.

###  II.1 Determination of the various rates

#### Thermal radius and geometric rates $C^{g}_{\star}$, $C^{g}_{\chi}$ and
$C^{g}_{2\chi}$:

Once a DM particle has been captured in the Sun, or a DMBS has formed, these
particles will thermalize with the SM material of the Sun and lie within
different spheres of thermal radius $r_{\rm th}$. Noting that the mass of a
bound-state is approximately twice that of a DM particle, these radii are
obtained Spergel and Press (1985); Press and Spergel (1985) by equating the
average thermal energy, $3T_{\odot}/2$, to the gravitational potential energy
per particle, $2\pi Gr^{2}_{\rm th}\rho_{\odot}nm_{\chi}/3$, as

$\displaystyle\\!\\!r_{\rm th}\\!=\\!\left(\frac{9T_{\odot}}{4\pi
G\rho_{\odot}nm_{\chi}}\right)^{\frac{1}{2}}\\!\\!=0.03R_{\odot}\\!\\!\left(\frac{T_{\odot}}{2.2\,{\rm
keV}}\frac{150\,{\rm\frac{g}{cm^{3}}}}{\rho_{\odot}}\frac{\rm
10\,GeV}{n\,m_{\chi}}\right)^{\frac{1}{2}}\\!\\!,$ $\displaystyle\\!\text{with
$n=1(2)$ for DM (DMBS)}.$ (6)

That is, the DM bound-state thermal radius is smaller than the DM thermal
radius by a factor of $\sqrt{2}$. We take the solar core temperature,
$T_{\odot}$, to be $2.2\,$keV, and the core mass density to be
$\rho_{\odot}\sim$150 g/cm3. For example, the thermal radius of $3$ GeV DM
particles is about one-tenth of the Solar radius, $R_{\odot}\simeq 6.9\times
10^{5}$ km Asplund  _et al._ (2009). In addition, we will assume throughout
that the DM radial distribution is isothermal, with its temperature $T\approx
T_{\odot}$ Press and Spergel (1985); Garani and Palomares-Ruiz (2017).

Upon taking into account the relative motion of the Sun with respect to the
galactic DM halo, the geometric capture rate on nucleons is Bottino  _et al._
(2002)

$C^{g}_{\star}\,\,{=\frac{\rho_{\chi}}{m_{\chi}}\,\pi
R_{\odot}^{2}\,\bar{v}}\simeq{6\times
10^{28}}\,{\mathrm{s}}^{-1}\left(\frac{\rho_{\chi}}{0.3\,\mathrm{GeV/cm^{3}}}\cdot\frac{10\,\mathrm{GeV}}{m_{\chi}}\right)\,.$
(7)

Here $\bar{v}$ is a factor with units of velocity, accounting for the relative
motion of the Sun and the velocity distribution of DM, see e.g. Garani and
Palomares-Ruiz (2017). Throughout this work we take the local DM density
$\rho_{\chi}=0.3$ GeV/cm3. Apparently, the geometric capture rate is
proportional to the corresponding cross-sectional area. So this expression can
be rescaled to obtain the geometric rate for DM self-capture on DM and on DMBS
(i.e. when the mean free path of DM is smaller than the corresponding thermal
sphere) by adding the corresponding factor $(r_{\rm th}/R_{\odot})^{2}$. This
results in geometric self-capture rates which scale as $m^{-2}_{\chi}$,

$\displaystyle C^{g}_{\chi}$ $\displaystyle=$ $\displaystyle 5\times
10^{25}\,{\rm s^{-1}}\left(\frac{10\,{\rm
GeV}}{m_{\chi}}\right)^{2}\,,\quad\quad C^{g}_{2\chi}\approx 2.6\times
10^{25}\,{\rm s^{-1}}\left(\frac{10\,{\rm GeV}}{m_{\chi}}\right)^{2}\,.$ (8)

####  Capture on nucleons rate $C_{\star}$:

The capture rate of DM particles in the Sun has been well studied in the
literature. It can be written as Feng  _et al._ (2016)

$\displaystyle C_{\star}$ $\displaystyle\approx$
$\displaystyle\sum_{i}\int_{0}^{R_{\odot}}4\pi r^{2}n_{i}(r){\rm
d}r\int_{0}^{\infty}{\rm
d}u_{\chi}\,\left(\frac{\rho_{\chi}}{m_{\chi}}\right)\times\,u_{\chi}\,\omega^{2}(r)f_{\odot}(u_{\chi})\int_{E_{R}^{min}}^{E_{R}^{max}}\frac{{\rm
d}\sigma_{i}}{{\rm d}E_{R}}\,{\rm d}E_{R}\,,$ (9)

where $f_{\odot}(u_{\chi})$ is the normalized asymptotic DM velocity
distribution far from the Sun in the solar frame. The values of radial number
density distribution of each element $i$, denoted as $n_{i}(r)$ above, are
adopted from the AGSS09 Solar model Asplund  _et al._ (2009). Together with
the escape velocity at a distance $r$ away from the centre of the Sun,
$v_{e}(r)$, it provides the relative velocity of a DM particle when it
scatters with the nucleus, $\omega(r)=\sqrt{u_{\chi}^{2}+v_{e}(r)^{2}}$. The
differential cross section, ${{\rm d}\sigma_{i}}/{{\rm d}E_{R}}$, encodes the
energy dependence of elastic scattering between DM and nucleus $i$ in the non-
relativistic limit. The recoil energy in the solar frame is given by
$E_{R}=\mu_{\chi N}^{2}/m_{N}\times w^{2}(1-\cos\theta_{CM})\,,$ where
$\mu_{\chi N}=m_{\chi}m_{N}/(m_{\chi}+m_{N})$ is the reduced mass of the DM-
nucleus system. For the Solar capture to actually occur, we need to make sure
that the recoil energy lies within the range of
$E_{R}^{min}<E_{R}<E_{R}^{max}$, with

$E_{R}^{min}=\frac{1}{2}m_{\chi}u_{\chi}^{2},\,\quad
E_{R}^{max}=\frac{2\mu_{\chi N}^{2}\omega^{2}}{m_{N}}\,.$ (10)

For definiteness we assume that the dark sector communicates with the SM
particles through the Higgs portal, and that the light mediator $\phi$ mixes
with the Higgs boson with a mixing angle $\theta_{\phi}$. The differential
cross section in the Born approximation is then

$\frac{{\rm d}\sigma}{{\rm
d}E_{R}}=\frac{g^{2}_{s}\cos^{2}{\theta_{\phi}}\sin^{2}{\theta_{\phi}}}{2\pi}\frac{m^{2}_{N}f^{2}_{N}}{v^{2}_{H}}\frac{m_{N}}{\omega^{2}(2m_{N}E_{R}+m^{2}_{\phi})^{2}}F^{2}\left(\frac{E_{R}}{Q_{0}}\right)~{}\,,$
(11)

where $F(E_{R}/Q_{0})$ is the Helm form factor and
$f_{N}=f^{N}_{u}+f^{N}_{d}+f^{N}_{s}+6/27f_{G}$, which is an $\mathcal{O}(1)$
number (see e.g. Lewin and Smith (1996); Duda  _et al._ (2007)).

The DM capture rates on nucleons are shown in Fig. 2, corresponding to a Higgs
mixing angle $\sin\theta_{\phi}=10^{-10}$ and $g_{s}=1$, for several mediator
masses (labeled in red). The geometric capture and self-capture rates, Eq. (7)
and (8), are shown as black dashed and dotted curves. For mediator mass below
MeV the differential cross section is independent of the mediator mass since
$m_{N}E^{max}_{R}\sim\mu_{\chi N}v^{2}_{e}$, which is above MeV for
$m_{\chi}\geq m_{\rm evap}$. This leads to a $t$-channel enhancement, which in
the $m_{\chi}\gg m_{N}$ limit gives a capture cross section scaling as
$m^{-1}_{\chi}$, leading to a capture rate scaling as $m^{-2}_{\chi}$. While,
for mediator masses larger than $\sim 100$ MeV and $m_{\chi}\lesssim m_{N}$,
the capture rate in turn scales as $m_{\chi}$. Note that in Eq. (11) we have
neglected the effect of DM-electron scattering, which is subleading due to
lower energy loss and tiny electron Yukawa coupling.

Figure 1: _Left_ : DM capture rate on nucleons for several mediator masses for
a Higgs mixing angle $\sin\theta_{\phi}=10^{-10}$. The geometric capture on
nucleons $C^{g}_{\star}$ and $C_{\chi}^{g}$ self-capture rates are shown by
the dashed and dashed-dotted black curves. The $C^{g}_{2\chi}$ self-capture
rate line (not shown) is a factor 2 below the $C_{\chi}^{g}$ line. _Right_ :
The DM BSF rate inside the Sun, for several values of $\alpha$.

Figure 2: _Left_ (_right_): DM self-capture rates due to DM-DM (DM-DMBS)
scattering. See text for details.

#### Bound-state formation rate $A_{bsf}$

As stated above, we consider DM to be asymmetric and made of Dirac fermions
($\chi$). The interaction between two (identical) DM particles is attractive
when the mediator is a scalar ($\phi$), with a classical Yukawa potential
$V(r)=-{\alpha\over r}e^{-m_{\phi}r}\,.$ The binding energy of the bound-state
is $E_{\rm bind}\simeq\frac{m_{\chi}\alpha^{2}}{4}-\alpha m_{\phi}$, with the
dark coupling $\alpha\equiv g_{s}^{2}/4\pi$. From the merging of two identical
fermions, the radiative BSF rate has been estimated in Wise and Zhang (2014),
yielding

$\sigma_{\rm BSF}v\simeq\frac{256\pi^{2}\alpha^{5}}{5e^{4}m_{\chi}^{2}v_{\rm
rel}}\,,$ (12)

in the limit of $\alpha/v_{\rm rel}\gg 1$ and $m_{\phi}\to 0$, where the DM
relative velocity is $v_{\rm rel}\simeq 2v_{\chi,\odot}$. The recoil velocity
of the DMBS is $\alpha^{2}/8$ in this limit, which must be below the escape
velocity from the Solar center, resulting in the requirement $\alpha\lesssim
0.18$.555Due to velocity enhancement, the BSF process in the early Universe
could be efficient such that very few free DM particles could exist in the
current epoch. Then qualitative changes in some regions of parameter space of
interest are expected. However, this depends on the assumed cosmological
history, and it is not explored in this work. See e.g. Wise and Zhang (2014);
Gresham  _et al._ (2017) for related discussions.  Finally, the total
annihilation rate that enters the number evolution equation is given by

$A_{\rm bsf}=\frac{\int\,n_{\chi}^{2}\,\sigma_{\rm BSF}v\,dV}{\left(\int
n_{\chi}dV\right)^{2}}~{}.$ (13)

The normalized DM radial distribution is

$n_{\chi}(r,t)=N_{\chi}(t)\,\frac{e^{-m_{\chi}\phi(r)/T}}{\int_{0}^{R_{\odot}}e^{-m_{\chi}\phi(r)/T}\,4\pi
r^{2}\,{\rm d}r}~{},$ (14)

which corresponds to an isothermal sphere, with a radial dependence set by the
gravitational potential
$\phi(r)=\int_{0}^{r}GM_{\odot}(r^{\prime})/{r^{\prime}}^{2}\,{\rm
d}r^{\prime}$, with $G$ the gravitational constant and $M_{\odot}(r^{\prime})$
the mass inside a sphere of radius $r^{\prime}$, and $N_{\chi}(t)$ is the
total population of DM particles at a given time $t$.

Once DM particles are thermalized its average velocity is
$v_{\chi}\simeq{}\sqrt{2T/m_{\chi}}$. The smallness of this value considerably
boosts the BSF in Eq. (12): $v_{rel}\simeq 4.2\times
10^{-3}\sqrt{1\,\hbox{GeV}/m_{\chi}}$. Consequently, $A_{\rm bsf}$ scales as
$\alpha^{5}(m_{\chi}^{2}v_{\rm rel})^{-1}(r^{\chi}_{\rm th})^{-3}$, be
approximately independent of the DM mass in practice, as shown in right panel
of Fig. 2.

####  DM self capture rates $C_{\chi}$ and $C_{2\chi}$:

Analogous to DM-nucleon capture rate, we estimate the self-capture rate per
target DM particle in the Sun as follows

$C_{\chi}=\int_{0}^{r^{\chi}_{th}}{\rm d}r4\pi
r^{2}n_{\chi}(r)\int_{0}^{\infty}{\rm
d}u_{\chi}\,\left(\frac{\rho_{\chi}}{m_{\chi}}\right)\,u_{\chi}\,f_{\odot}(u_{\chi})\,w^{2}\int_{\cos\theta_{\rm
min}}^{\cos\theta_{\rm max}}\frac{{\rm d}\sigma_{\chi-\chi}}{{\rm
d}\cos\theta}\,{\rm d}\cos\theta~{}.$ (15)

The minimum scattering angle is set by the requirement of a minimum energy
that has to be lost in a single scattering event to be captured, i.e. DM
kinetic energy at infinity must be smaller than $1/2\,m_{\chi}u^{2}_{\chi}$.
The maximum scattering angle for scattering of identical particles is $\pi/2$.
Therefore $\cos\theta_{\rm min}=1-2\frac{u^{2}_{\chi}}{\omega^{2}}$ and
$\cos\theta_{\rm max}=0$ . Similar to the DM-DM self-capture, the capture rate
due to DM-DMBS scattering, per target DMBS in the Sun, is given by

$C_{2\chi}=\int_{0}^{r^{2\chi}_{th}}{\rm d}r4\pi
r^{2}n_{2\chi}(r)\int_{0}^{\infty}{\rm
d}u_{\chi}\,\left(\frac{\rho_{\chi}}{m_{\chi}}\right)\,u_{\chi}\,f_{\odot}(u_{\chi})\,\omega^{2}\int_{\cos\theta_{\rm
min}}^{\cos\theta_{\rm max}}\frac{{\rm d}\sigma_{2\chi-\chi}}{{\rm
d}\cos\theta}\,F_{\chi}^{2}\left(\frac{E_{r}}{Q_{\chi}}\right){\rm
d}\cos\theta~{},$ (16)

with $\cos\theta_{\rm min}=1-\frac{9}{8}\frac{u_{\chi}^{2}}{\omega^{2}}$ and
$\cos\theta_{\rm max}=-1$. We have used $m_{2\chi}\approx 2m_{\chi}$ and set
that the coupling to bound-state is $2\alpha$ (as the vertex is of the scalar
type). We assume that the form factor has the form of $\exp(E_{r}/Q_{\chi})$.
The typical size of the bound-state is set by the Bohr radius, hence
$Q_{\chi}=m_{\chi}\alpha$. To a good approximation $F_{\chi}\rightarrow 1$ for
a typical value of $\alpha=0.1$. The radial number density of DMBS
($n_{2\chi}$) is assumed to be isothermal, analogous to $n_{\chi}$.

Bound-state formation suggests that non-perturbative effects are non-
negligible for the calculation of the DM self-scattering cross section (see
e.g. Tulin  _et al._ (2013); Colquhoun  _et al._ (2021)). This is indeed our
case. For the parameter values of interest $\alpha\in[0.02,0.2]$, $m_{\chi}>5$
GeV, ${\rm MeV}\lesssim m_{\phi}\ll m_{\chi}$, and $v\sim v_{esc}\sim 4\times
10^{-3}$, we are never in the Born regime of self-scattering, for which
$m_{\chi}<m_{\phi}/\alpha$ or $\alpha\lesssim v$ Landau and Lifshits (1991).
As a result, a perturbative expansion of the scattering cross section is not
justified here, see left panel of Fig. 7 in Appendix B.1. However, for these
parameters above, it is not necessary to solve the Schrödinger equation
associated with the non-relativistic DM collisions because the scattering
process is typically semi-classical, for which the range of the potential is
larger than the DM de Broglie wavelength, that is, $m_{\chi}v/m_{\phi}\gtrsim
1$. In this case classical mechanics can be employed to estimate the
scattering cross section.666A precise treatment would require to solve the
phase shift for a large number of partial waves, see e.g. Tulin  _et al._
(2013); Chen  _et al._ (2015); Colquhoun  _et al._ (2021).

More importantly, unlike phenomenological studies of DM self-scattering in
galaxies and galaxy clusters, where the transfer and viscosity cross sections
are relevant, here the differential cross section is needed to calculate the
integrated self-capture rates introduced above. This is because DM capture
needs that enough initial kinetic energy is lost in a single scattering. This
condition is imposed by integrating the differential cross section from
minimum possible scattering angle ($\theta_{\rm min}$, set by required energy
loss) to the maximum one where the scattered particle still does not gain
enough energy to escape. These limits are explicitly indicated in Eqs. (15)
and (16). To calculate the differential cross sections, we follow Refs.
Khrapak  _et al._ (2003, 2004) and solve the elastic scattering with a Yukawa
potential classically. Further discussions and full expressions are presented
in Appendix B.

We present the results for DM self capture rates through scattering on free DM
particles (DM bound-states) in the left (right) panel of Fig. 2. The rates
scale proportionally as $m^{-2}_{\chi}$ for $m_{\chi}\gg m_{\phi}$. For
moderate values of DM masses the scaling is less steep. Note that the capture
rate on bound-states $C_{2\chi}$ is approximately larger by a factor of two
with respect to $C_{\chi}$, due to larger Yukawa coupling induced by bound-
states and larger maximal allowed scattering angle that keep both particles
gravitationally captured.

###  II.2 Integrating the Boltzmann equations

The set of coupled Boltzmann equations of Eqs. (2) and (3) has no closed form
analytical solution. In Appendix A we describe at length how they can be
solved approximately. We summarize in this subsection the main outcome of this
discussion. To this end, in Fig. 3 we present the evolution of both
populations for two parameter sets. The evolution of $N_{\chi}$ ($N_{2\chi}$)
from the numerical integration of the full Boltzmann Eqs. (2) and (3) are
given by the solid (dashed) red curves. Curves with other colors correspond to
what is obtained switching off both the $C_{\chi}$ and $C_{2\chi}$ terms
(black), or only the $C_{\chi}$ term (blue) or only the $C_{2\chi}$ term
(green).

Figure 3: Number evolution of DM and DM bound-states for two parameter sets,
as a function of time. Here $\tau_{s}$ and $\tau_{0}$ are characteristic times
explained in the main text (and exactly defined in Appendix A), $t_{\odot}$ is
the age of the Sun. The color code of the curves in the left panel are the
same as that in the right panel.

At early times the number of free particles $N_{\chi}$ grows as $C_{\star}\,t$
(i.e. red and black solid lines coincide). As $N_{\chi}$ increases, the BSF
starts to occur efficiently and the corresponding term (quadratic in
$N_{\chi}$) quickly catches up with the constant term associated with capture
on nucleons, so that in absence of the $C_{2\chi}$ term, a quasi-static
equilibrium between both terms is reached (see black and green solid curves in
Fig. 3). Solving the Boltzmann equation for $N_{\chi}$ one gets that the
associated time scale is $\tau_{0}\equiv(C_{\star}A_{\rm bsf})^{-1/2}$ if one
drops the $C_{\chi}$ term –which has a subleading effect (black solid curve)–
or $\tau_{s}\equiv(C_{\star}A_{\rm bsf}+C^{2}_{\chi}/4)^{-1/2}$ if one
includes it (green solid curve). For more details, see the end of this section
and Appendix A.

Another consequence of efficient BSF is that the number of bound-states
$N_{2\chi}$ quickly catches up to the number of DM particles $N_{\chi}$, and
never stops growing thereafter, as it is not counterbalanced by any other term
in Eq. (3). From this point, the term associated with the capture on bound-
states, $C_{2\chi}N_{2\chi}$, takes over and leads to an exponential growth of
$N_{\chi}$ (red and blue solid curves in Fig. 3) and of $N_{2\chi}$ (dashed
curves with same colors). The exponential growth become significant at
$t\gtrsim\tau_{s}+C_{2\chi}^{-1}$, or at $\tau_{0}+C_{2\chi}^{-1}$, depending
on whether one includes the $C_{\chi}$ term. This growth lasts until the
capture rate saturates the geometric rate within the bound-state thermal
sphere, i.e. when the term $C_{2\chi}N_{2\chi}$ reaches $C^{g}_{2\chi}$. This
happens when the mean free path of DM particles becomes much smaller than the
bound-state thermal radius. At this moment, we can safely neglect the
$C_{\chi}N_{\chi}$ term as $N_{\chi}\ll N_{2\chi}$, thus the Boltzmann
equation takes the form, $dN_{\chi}/dt\simeq C_{\star}+C^{g}_{2\chi}-A_{\rm
bsf}N_{\chi}^{2}$. This means that, quickly after the geometric capture rate
is saturated, the BSF term compensates the constant capture rate from both
nucleon and bound-state terms. Setting $dN_{\chi}/dt\simeq 0$ gives the
maximum final value

$N_{\chi,{\rm eq}}\simeq\left(\frac{C_{\star}+C_{2\chi}^{g}}{A_{\rm
bsf}}\right)^{1/2}~{}.$ (17)

We denote as $\tau_{g}$ the time when $N_{\chi}$ freezes in such a way. An
analytic approximation is given in Appendix A.

Note that, contrary to the $C_{2\chi}N_{2\chi}$ term –which saturates the
geometric rate $C_{2\chi}^{g}$– the average DM density within the thermal
radius $r^{\chi}_{\rm th}$ is smaller than the DMBS density within
$r^{2\chi}_{\rm th}$, so the self-capture term $C_{\chi}N_{\chi}$ does not
reach $C^{g}_{\chi}$ at $t=\tau_{g}$, and will never reach it after, as
$N_{\chi}$ stops increasing. As for $N_{2\chi}$, it keeps increasing in time
forever. For $t>\tau_{g}$, it increases linearly in time and the rate of BSF
is half the capture rate,

$\Gamma(t)=\frac{dN_{2\chi}}{dt}=\frac{1}{2}A_{\rm bsf}N^{2}_{\chi,{\rm
eq}}\simeq\frac{C_{\star}+C^{g}_{2\chi}}{2}\,,$ (18)

inducing a flux of mediators given by the same rate.

The reason why the $C_{\chi}$ term has a subleading effect with respect to the
effect of the $C_{2\chi}$ term stems from the fact that these two terms are
very different in nature. If we switch off the $C_{2\chi}$ term (green
curves), the $C_{\chi}$ term also leads to an exponential grow, starting
slightly before $t=\tau_{s}$, but it is much less important than the
exponential growth from the $C_{2\chi}$ term. This is analogous to what
happens for the symmetric DM case with BSF playing the role of annihilation if
$C_{2\chi}=0$.777In the limit of $C_{2\chi}\rightarrow 0$ and $A_{\rm bsf}=0$
(and with an additional DM annihilation term), the solution has an analytical
form, first presented in Zentner (2009).  In Eq. (3) the $C_{\chi}$ term
(linear in $N_{\chi}$) is quickly counter-balanced by formation of bound-
states from the $A_{\rm bsf}$ term (quadratic in $N_{\chi}$): after a short
period of exponential growth, $N_{\chi}\sim
C_{\star}(e^{C_{\chi}t}-1)/C_{\chi}$, the green curve ($C_{2\chi}=0$), goes to
a plateau from an equilibration of capture and BSF. In contrast, if
$C_{2\chi}\neq 0$, the capture rate on bound-state grows as $N_{2\chi}$. Thus,
due to capture on bound-state term $C_{2\chi}$, if it were not for the
geometric rate upper floor, both $N_{\chi}$ and $N_{2\chi}$ would grow
forever. Moreover, the argument of the exponential from the $C_{2\chi}$ term
is larger than the one from the $C_{\chi}$ term, because $C_{2\chi}$ is about
a factor of two times $C_{\chi}$, see Fig. 2.

In presence of both $C_{\chi}$ and $C_{2\chi}$ terms, the $C_{\chi}$ term has
only a moderate impact on the evolution of $N_{\chi}$ and $N_{2\chi}$. As a
comparison of the full evolution (red curves in Fig. 3) and the $C_{\chi}=0$
evolution (blue curves) shows, including $C_{\chi}$ reduces the timescale when
the exponential growth starts, from $t\sim\tau_{0}$ to $t\sim\tau_{s}$. Thus,
due to the $C_{\chi}$ term, the exponential effect of the $C_{2\chi}$ term
starts somewhat earlier, and the geometric rate within the bound-state thermal
sphere is also reached somewhat earlier. The values of $N_{\chi}$ and
$N_{2\chi}$ are insensitive to the $C_{\chi}$ term at $t\gg\tau_{s}$.

##  III Neutrino flux and Terrestrial detection

As explained above, thanks to the capture on DMBS, the capture rate will
saturate the geometric rate within the DMBS thermal sphere, leading to an
equilibrium between the capture and the BSF processes. The time $\tau_{g}$, at
which this happens, can be easily shorter than the age of the Sun, as solved
in Eqs. (44) and (45). At this point, the flux of mediators emitted from BSF
becomes constant and is approximately given by Eq. (18) above.

For concreteness (and other reasons explained below) we will assume in the
following that the mediator decays dominantly into a pair of neutrinos. Under
the assumption that the mediator mass is much smaller than the binding energy,
$E_{\rm bind}$, the differential flux at a terrestrial detector is

${d^{2}\Phi_{\nu}\over dE_{\nu}d\Omega}={A_{\rm
bsf}N_{\chi}^{2}(t_{\odot})\over 2\Delta\Omega}{dN_{\nu}\over
dE_{\nu}}={A_{\rm bsf}N_{\chi}^{2}(t_{\odot})\over 8\pi
d_{\odot}^{2}\left(1-\cos\vartheta\right)}{dN_{\nu}\over dE_{\nu}}\,,$ (19)

where $d_{\odot}$ is the Earth-Sun distance (AU), with $\vartheta$ being the
angular sensitivity of detector (as the apparent angular diameter of the Sun
$\vartheta_{\odot}\approx 0.5^{\circ}$). Since the light mediator is boosted,
the neutrino energy spectrum this decay leads to is not monochromatic but has
a characteristic box shape Ibarra  _et al._ (2012)

$\displaystyle{dN_{\nu}\over dE_{\nu}}\simeq{2\over E_{\rm
bind}}\text{~{}~{}~{}for~{}}E_{-}\leq E_{\nu}\leq E_{+}\,,$ with
$\displaystyle E_{\pm}=\frac{E_{\rm
bind}}{2}\Big{(}1\pm\sqrt{1-\frac{4m^{2}_{\phi}}{E^{2}_{bind}}}\Big{)}\,.$
(20)

In particular, $E_{-}\simeq 0$ and $E_{+}\simeq E_{\rm bind}$ if $m_{\phi}\ll
E_{\rm bind}$. The differential flux of neutrinos at the detector (ignoring
oscillations) is

$\displaystyle\frac{d\Phi_{\nu}}{dE_{\nu}}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi
d^{2}_{\odot}}\Gamma(m_{\chi},m_{\phi},\alpha_{\chi},s_{\theta})\frac{2}{E_{\rm
bind}}\Theta(E_{+}-E_{\nu})\Theta(E_{\nu}-E_{-})~{}.$ (21)

The resulting neutrino fluxes are shown in Fig. 4 for a fixed Higgs mixing
angle $\sin\theta_{\phi}=10^{-12}$, for $\alpha=0.05,\,0.1$ and $0.15$ (in
blue, red and maroon colors), respectively. For comparison, the current limits
(90% C.L.) on diffuse supernova $\nu_{e}$ neutrino background (DSNB) flux from
Super-K runs III, IV and KamLAND are shown in blue, red and black points
Richard  _et al._ (2016); Abe  _et al._ (2021, 2022). Also shown are predicted
atmospheric fluxes of $\nu_{e}$ down to 10 MeV in thin green line, assuming
$30^{\circ}$ angular resolution in the sky, adapted from Battistoni  _et al._
(2005); Zhou and Beacom (2023); Suliga and Beacom (2023). Current measurements
of atmospheric neutrino fluxes are shown in green points Honda  _et al._
(2015). As this figure suggests, the neutrino flux induced by BSF could be
observable. When the geometric limit for DMBS sphere is reached, the emitted
flux of mediators is proportional to $C^{g}_{2\chi}\propto m^{-2}_{\chi}$.
This scaling is seen in Fig. 4. If the geometric limit is not reached, heavier
mediators would lead to smaller self-capture rates, as shown by the falling
tails of the neutrino flux.

Figure 4: _Left_ : Flux of light mediators emitted by the BSF process for
various parameter choices. _Right_ : Neutrino flux as a function of neutrino
energy, coming from the Sun for a detector placed on the surface of the Earth.
Shown in red, blue and black scattered points are the current limits on
diffuse supernova neutrinos adapted from Abe  _et al._ (2021). Atmospheric low
energy neutrino measurements from Super-K are shown by the green dots Richard
_et al._ (2016). Thin green line denotes the predicted atmospheric neutrino
flux for $30^{\circ}$ sky Suliga and Beacom (2023).

Finally, note that our results can be easily generalized, e.g. if the mediator
decays electromagnetically. Nevertheless, the experimental constraints, as
studied below, also become more stringent for electromagnetic decays, further
narrowing down allowed regions of parameter space.

##  IV Results and constraints for an explicit model

###  IV.1 Model

As mentioned above, we consider DM in the form of a (vector-like) Dirac
fermion self-interacting through the exchange of a light scalar. Concretely
the Lagrangian is (see e.g. Krnjaic (2016); Winkler (2019); Matsumoto  _et
al._ (2019))

$\displaystyle{\cal L}$ $\displaystyle\supset$ $\displaystyle{\cal L}_{\rm
SM}+{\cal
L}_{N_{R}}+\bar{\chi}\left(i\not{\partial}-m_{\chi}\right)\chi+\frac{1}{2}\left(\partial\phi\right)^{2}-g_{s}\phi\bar{\chi}\chi-V(\phi,H)\,.$
(22)

Here we do not specify the scalar potential $V(\phi,H)$, and simply assume
that it induces a $\phi$-$H$ mixing angle $\theta_{\phi}$ (which can be
achieved in various ways), so that, upon rotation to the mass basis, the
following interactions are obtained:

$\displaystyle{\cal L}_{\rm int}$ $\displaystyle=$ $\displaystyle
g_{s}\cos\theta_{\phi}\phi^{\prime}\bar{\chi}\chi+g_{s}\sin\theta_{\phi}h^{\prime}\bar{\chi}\chi+\left(\cos\theta_{\phi}h^{\prime}+\sin\theta_{\phi}\phi^{\prime}\right)\frac{m_{f}}{v_{H}}\bar{f}_{L}f_{R}\,.$
(23)

If we assume that neutrino masses are generated through the usual type-I
seesaw mechanism, nothing prevents the right-handed neutrinos to couple in
pairs to the light mediator $\phi$. Thus one has the extra interactions

$\displaystyle-{\cal L}_{N_{R}}$ $\displaystyle=$ $\displaystyle
Y_{\nu}\bar{N}L\cdot
H+\frac{m_{N}}{2}\overline{N^{c}}N+Y_{\phi}\phi\overline{N^{c}}N+{\rm
h.c.}~{}.$ (24)

In the physical $\nu^{\prime},N^{\prime}$ and $\phi^{\prime},h^{\prime}$ mass
eigenstate basis, this leads to the following relevant Yukawa interactions for
the light scalar eigenstate $\phi^{\prime}$, up to second order in the
neutrino mixing angle,

$\displaystyle-{\cal L}_{\rm N_{R}}$ $\displaystyle\supset$ $\displaystyle
Y_{\phi}\phi^{\prime}\cos\theta_{\phi}(\overline{N^{\prime
c}}N^{\prime}\cos^{2}\theta_{\nu}+\overline{N^{\prime c}}\nu^{\prime
c}\sin\theta_{\nu}+\overline{\nu^{\prime}}N^{\prime}\sin\theta_{\nu}+\overline{\nu^{\prime}}\nu^{\prime
c}\sin^{2}\theta_{\nu})+h.c.$ (25) $\displaystyle-$ $\displaystyle
Y_{\nu}\phi^{\prime}\sin\theta_{\phi}(\bar{N^{\prime}}\nu^{\prime}\cos^{2}\theta_{\nu}+\overline{\nu^{\prime
c}}\nu^{\prime}\sin\theta_{\nu}-\bar{N^{\prime}}N^{\prime
c}\sin\theta_{\nu}-\overline{\nu^{\prime c}}N^{\prime
c}\sin^{2}\theta_{\nu})+h.c.\,.$

Thus, both the $Y_{\phi}$ and $Y_{\nu}$ interactions induce a decay of the
light mediator into a pair of light neutrinos. In practice the decay induced
by $Y_{\phi}$ will be dominant because various constraints require
$Y_{\nu}\sin\theta_{\phi}\ll Y_{\phi}\sin\theta_{\nu}$. In the following, for
simplicity, we will consider only one left-handed neutrino and one right-
handed Majorana neutrino, with the neutrino mixing angle
$\sin\theta_{\nu}\approx Y_{\nu}v_{h}/\sqrt{2}M_{N}$. Likewise for the
associated Yukawa coupling, the typical seesaw gives

$\sin^{2}\theta_{\nu}\sim\frac{m_{\nu}}{m_{N}}=10^{-13}\cdot\frac{m_{\nu}}{0.1\hbox{eV}}\cdot\frac{1\,\hbox{TeV}}{m_{N}}$
(26)

with $Y_{\nu}\sim\sqrt{2m_{\nu}m_{N}/v_{h}^{2}}$.

If $m_{N}>m_{\phi}$ the light mediator cannot decay into $NN$ or $N\nu$ and
the $\nu\nu$ is the only possible channel induced by the seesaw interactions.
The decay width is

$\Gamma(\phi\rightarrow\nu\nu)=\frac{1}{16\pi}|Y_{\phi}|^{2}\,m_{\phi}\,\sin^{4}\theta_{\nu}\,.$
(27)

In addition to this neutrino-mixing channel, $\phi$ can also decay, through
the $\phi$-$H$ scalar mixing, into pairs of charged leptons or quarks (or SM
bosons for large $\phi$ mass). As will be discussed below, this needs to be
suppressed.

###  IV.2 Constraints and results

An observable neutrino flux requires sufficiently large capture and BSF rates,
eventually leading to an equilibrium between both processes. This scenario
must also fulfill additional phenomenological constraints:

Figure 5: Constraints for two example values of the dark coupling $\alpha$.
Contour lines give the values of $\sigma_{\rm SI}/m_{\chi}$ in dwarf-sized
halos. Black (gray) hatched regions are excluded because $\sigma_{\rm
SI}/m_{\chi}\geq 0.5$ cm2/g at cluster scales (or 100 cm2/g at dwarfs scales,
see contour lines). The resonant regime (with peaks) will become smaller in
the Milky Way halo, where the average DM velocity is much larger than that at
dwarfs scale.

Self-interactions: On the one hand, the corresponding self-scattering cross
section is bounded from above by the observation of galaxy cluster collisions,
e.g., Harvey  _et al._ (2015); Bondarenko  _et al._ (2018); Harvey  _et al._
(2019); Sagunski  _et al._ (2021); Cross  _et al._ (2023)

${\sigma_{\rm SI}\over m_{\chi}}\,\lesssim\,0.5\,\text{cm}^{2}/\text{g}~{},$
(28)

for DM velocities around $v\sim 1000\,$km/s. On the other hand, the non-
observation of gravothermal collapse in dwarf-sized halos sets an upper bound
of about $100\,\text{cm}^{2}/\text{g}$ at small scales, for which we take
$v\sim 25\,$km/s Tulin and Yu (2018); Adhikari  _et al._ (2022). Concretely,
we derive the bounds taking $\sigma_{\rm SI}$ equal to the viscosity cross
section of DM particles. Adopting the modified transfer cross section barely
changes the results, see Appendix B.1 for detailed discussions of these
effective cross sections.

Fig. 5 shows the resulting values for this cross section as a function of
$m_{\chi}$ and $m_{\phi}/m_{\chi}$ for two values of $\alpha$. The cluster
bound excludes $\chi$ masses below $\sim 10$-$1000$ GeV depending on
$m_{\phi}/m_{\chi}$ and $\alpha$. This figure also shows the region where
$m_{\phi}$ is larger than the binding energy where BSF is not possible. Values
of the cross section that could address the small scale anomalies are allowed,
i.e. $\sigma_{\rm SI}/m_{\chi}\in[1,\,100]\,\text{cm}^{2}/\text{g}$ at typical
DM velocities of the order of 25 km/s in dwarf galaxies Vogelsberger  _et al._
(2012); Zavala  _et al._ (2013); Elbert  _et al._ (2015); Ren  _et al._
(2019). In the $m_{\chi}$-$\alpha$ plane, Fig. 6 (right) shows the
corresponding upper bound on $\alpha$.

Figure 6: _Left_ : Upper bounds on the DM-SM interaction strength from XENON1T
Aprile  _et al._ (2018), depending on $m_{\phi}$. _Right_ : Upper bounds on DM
coupling $\alpha$ from indirect search for neutrino flux (purple shaded area)
and DM self-scattering at dwarfs scales. The bounds from the Bullet cluster,
Eq. (28), are given for $m_{\phi}=1\,$MeV and $10\,$MeV. The conservative
dwarf scale upper bound, ${\sigma_{\rm
SI}/m_{\chi}}\,\lesssim\,100\,\text{cm}^{2}/\text{g}$ is given for
$m_{\phi}=1$ MeV. This case is disfavored by BBN observables, given its
sizeable coupling to neutrinos. For $m_{\phi}=100\,$MeV, the dwarf galaxy
scale bound is weaker and basically irrelevant for BSF. For $\alpha$ above the
gray solid line, one expects a galactic flux from BSF in the galactic halo
larger than the flux from BSF in the Sun. Efficient BSF in the early Universe
is possible for parameters above the dot-dashed gray line Gresham  _et al._
(2017).

BBN: In order to observe the neutrino flux in the detector, the in-flight
lifetime of the mediator must be smaller than eight minutes. This constraint
is typically less stringent than that resulting from upper bounds on extra
radiation after neutrino decoupling, which requires its lifetime to be shorter
than one second. Using the typical seesaw expectation of Eq. (26), one gets

$|Y_{\phi}|^{2}\sim
0.03\,\,\frac{100\,\hbox{MeV}}{m_{\phi}}\frac{1\,\hbox{s}}{\tau_{\phi}}\,\Big{(}\frac{0.1\hbox{eV}}{m_{\nu}}\frac{m_{N}}{1\hbox{GeV}}\Big{)}^{2}\,.$
(29)

Thus a fast enough decay is obtained provided that $Y_{\phi}$ is not too small
and $m_{N}$ is not too large. BBN and perturbative couplings,
$Y_{\phi}^{2}\lesssim 4\pi$, require $m_{N}\lesssim
20\,\hbox{GeV}\cdot(m_{\nu}/0.1\,\hbox{eV})\cdot(m_{\phi}/100\,\hbox{MeV})^{1/2}$.
Thus low scale seesaw is favored along this scenario. This BBN bound can be
relaxed if the $\phi$ number density gets suppressed before neutrino
decoupling.

Direct detection: As already mentioned above, due to evaporation, DM masses
below a few GeV are not relevant for our purpose. Thus, non-observation of
spin-independent nucleon recoil signals provides the best direct-detection
bounds. For purely scalar-mediated interactions, the differential elastic
scattering cross section with nuclei has been given by Eq. (11). If
$m_{\phi}\lesssim\sqrt{2m_{N}E_{R}}\lesssim 1\,\mathrm{MeV}$ for typical keV-
scale recoil energies, direct detection rates, albeit in the Born regime, are
boosted due to $t$-channel exchange of the light mediator, leading to more
stringent bounds than for heavier mediators. The corresponding measurements by
the Xenon1T experiment Aprile  _et al._ (2018); Athron  _et al._ (2019) set an
upper bound on $g_{s}\sin\theta_{\phi}$, which we show in Fig. 6 (left) for
different masses of mediator and DM. This obviously translates into an upper
bound on the capture rate on nucleons, as both processes involve the same DM-
nucleon cross section. Recent LZ data can improve the limit on
$g_{s}\sin\theta_{\phi}$ by a factor of $2$ – $3$, depending on the mediator
mass Aalbers  _et al._ (2023).

Indirect detection from galactic center emission: If DM BSF occurs in the Sun,
it is reasonable to anticipate its occurrence in the galactic center of the
Milky Way as well. Consequently, we would expect corresponding emission of
energetic neutrinos originating from the galactic center, also with a box-
shaped energy spectrum. Here we estimate the indirect search limit by re-
scaling the current bounds on symmetric DM from neutrino telescope
observations. Regarding indirect signals, a BSF process that generates one
mediator and eventually two neutrinos with $E_{\nu}\approx E_{\rm bind}/2$ is
equivalent to symmetric DM annihilation with a mass of $E_{\rm bind}/2$, if
the latter only makes up a fraction of the observed DM abundance. Therefore,
existing bounds can be rescaled as follows,

$(\sigma_{\rm symm.}v)|_{{m_{\chi}\to E_{\rm bind}/2}}\rightarrow\left({E_{\rm
bind}/2\over m_{\chi}}\right)^{2}(\sigma_{\rm BSF}v)\,,$ (30)

where $(\sigma_{\rm symm.}v)$ denotes the known bounds on symmetric DM as a
function of the DM mass.

Quantitatively, the neutrino flux generated from the BSF process of halo DM
particles is estimated to be

$\Phi^{\rm D}_{\nu}\sim 1\,\text{cm${}^{-2}$s${}^{-1}$}\left({\text{GeV}\over
m_{\chi}}\right)^{2}\left({\sigma_{\rm BSF}v\over 3\times 10^{-22}\,{\rm
cm}^{3}/{\rm s}}\right)\,.$ (31)

This is similar to boosted DM case Agashe  _et al._ (2014), which can vary
mildly due to the uncertainty of the Galactic J-factor Cirelli  _et al._
(2011). The corresponding upper bound on $\alpha$ is shown in the right panel
of Fig. 6, obtained from re-scaling the indirect bounds holding in the
symmetric DM case Argüelles  _et al._ (2021). As also illustrated in Fig. 5,
this excludes values of $m_{\chi}$ between $3.0$ and $11.4$ GeV for
$\alpha=0.15$, and gives irrelevant bound for $\alpha=0.05$.

We can further speculate when this galactic neutrino flux is larger than the
one induced by BSF inside the Sun. The latter is approximately
$\Phi^{\odot}_{\nu}/({\rm GeV}/m_{\chi})^{2}\sim\mathcal{O}(1)$ cm-2s-1 when
the DM capture on DMBS is saturated by its geometric rate. Using Eq. (12)
results in $\alpha^{5/2}/m_{\chi}\lesssim 10^{-4}/{\rm GeV}$, which in turn
requires $m_{\chi}\gtrsim 5.6\,$GeV (87.1 GeV) for $\alpha=0.05$ (0.15), see
also the right panel of Fig. 6. In practice, the non-vanishing mass of the
mediator suppresses the BSF cross section, reducing the neutrino flux, but
barely affects the BSF-induced flux from the Sun.

CMB: Since the CMB observables are mostly sensitive to the total
electromagnetic energy injected by extra processes in the high-redshift
Universe, the associated constraint on BSF can be obtained by re-scaling the
CMB bound on the annihilation of a symmetric DM candidate Elor  _et al._
(2016),

$f_{\rm eff}{\sigma_{\rm BSF}v\over m_{\chi}}\,\frac{{\rm
Br}_{\phi\to\hbox{\tiny EM}}E_{\rm bind}}{2m_{\chi}}\lesssim 4.1\times
10^{-28}\,\frac{\hbox{cm}^{3}/\hbox{s}}{\hbox{GeV}}\,.$ (32)

Here, $f_{\rm eff}$ is an efficiency factor that depends on the spectrum of
injected electrons and photons, and we have taken into account that the energy
injected per process is not about $2m_{\chi}$ as in the case of DM
annihilation, but is given by the binding energy
$\approx\alpha^{2}m_{\chi}/4$, times the electromagnetic branching ratio of
$\phi$ decay, ${\rm Br}_{\phi\to\text{\tiny EM}}$. For $m_{\phi}\ll E_{\rm
bind}$, the mass of the light mediator does not play any role in the CMB
bound. For decay into SM particles other than neutrinos, the fact that BSF is
strongly enhanced when the DM relative velocity is small gives a CMB bound
which is much stronger than other bounds derived from local cosmic-ray
observations such as Fermi-Lat and AMS experiments. For instance, taking the
DM velocity at CMB to be a few km/s leads to ${\rm
Br}^{1/3}_{\phi\to\text{\tiny EM}}\alpha^{7/3}\lesssim
10^{-5}m_{\chi}/\text{GeV}$. Combining this condition with Eq. (27) gives the
constraint $Y_{\phi}\theta_{\nu}^{2}/\theta_{\phi}\gtrsim
10^{2.2}\alpha^{7/2}(\text{GeV}/m_{\chi})^{3/2}$ for $m_{\phi}$ below twice
the muon mass. For typical seesaw values of the neutrino mixing angle, Eq.
(26), this translates to the condition $m_{N}\lesssim 22\,\hbox{TeV}\cdot
Y_{\phi}\,(m_{\chi}/5\,\hbox{GeV})^{3/2}(10^{-12}/\theta_{\phi})(0.1/\alpha)^{7/2}(m_{\nu}/0.1\,\hbox{eV})$.

Mediator decay into right-handed neutrinos: If $2m_{N}<m_{\phi}$ the dominant
decay is not anymore into a pair of light neutrinos but into a pair of right-
handed neutrinos (or into $\nu N$ for $m_{\phi}/2<m_{N}<m_{\phi}$), with their
decay widths given by,

$\displaystyle\Gamma(\phi\rightarrow
NN)=\frac{1}{16\pi}|Y_{\phi}|^{2}\,m_{\phi}\,,$
$\displaystyle\Gamma(\phi\rightarrow
N\nu)=\frac{1}{16\pi}|Y_{\phi}|^{2}\,m_{\phi}\,\sin^{2}\theta_{\nu}\,.$ (33)

Using the typical seesaw expectation of Eq. (26), one gets

$\displaystyle|Y_{\phi}|^{2}\sim
10^{-22}\,\,\frac{1\,\hbox{GeV}}{m_{\phi}}\frac{1\,\hbox{s}}{\tau_{\phi}}\,,$
$\displaystyle|Y_{\phi}|^{2}\sim
10^{-12}\,\,\frac{1\,\hbox{GeV}}{m_{\phi}}\frac{1\,\hbox{s}}{\tau_{\phi}}\,\frac{0.1\hbox{eV}}{m_{\nu}}\frac{m_{N}}{1\hbox{GeV}}$
(34)

for $2m_{N}<m_{\phi}$ and $m_{\phi}/2<m_{N}<m_{\phi}$, respectively. This can
be compatible with the extra radiation constraint $\tau_{\phi}\lesssim 1$ s
but is basically excluded by the CMB bound because the right-handed neutrino
decay product contains a non-negligible amount of electromagnetic material.

##  V Summary

We have considered the possibility that asymmetric DM forms bound-states in
the Sun, and showed that this leads to novel phenomenology. BSF in the Sun can
proceed via emission of light scalar particles that carry energy roughly equal
to the binding energy. Their decays to neutrinos lead to potentially testable
low energy signals at neutrino detectors.

Unlike for annihilating DM, BSF produces a flux of particles without reducing
the number of DM particles in the Sun. We point out that on top of the DM
particles captured in the Sun, the bound-states piling up in this way become
additional scattering targets through which DM from the galactic halo could be
captured. We have determined the associated DM accretion rates on DM and DM
bound-states by evaluating the differential cross section, taking into account
that for typical parameters, $v\sim 10^{-3}$, $\alpha\sim 0.1$, $m_{\phi}<$
GeV and $m_{\chi}\sim 100$ GeV, the DM-DM and DM-DMBS scattering processes
proceed in the semi-classical regime. As soon as these self capture rates are
larger than $t_{\odot}^{-1}$, they can become phenomenologically relevant. In
particular, we have shown that, thanks to the self-capture on bound-states,
the number of DM particles in the Sun can exponentially increase, so much that
the capture rate can reach the geometric rate, i.e. all the DM particles
intercepting the DM bound-state thermal sphere are captured as the mean free
path becomes smaller than the sphere. As a result, this exponential effect
also considerably boosts the BSF and thus the associated flux of light
mediators. In an example model, where DM is a Dirac fermion which self-
interacts through exchange of light scalar that mixes with the Higgs boson,
with the scalar decaying into two neutrinos through seesaw interactions, this
leads to a neutrino flux which reach the predicted atmospheric neutrino fluxes
at energies below hundred MeV. Near future experiments such as Hyper-K, as
well as direct detection experiments, will be able to probe further this
scenario.

#### Acknowledgments:

We thank Sergio Palomares-Ruiz and Hai-bo Yu for discussions, and Sebastian
Wild for help with the implementation of direct detection rates in DDCalc
Athron  _et al._ (2019). X.C. is supported by the Research Network Quantum
Aspects of Spacetime (TURIS) and co-funded by the European Union (ERC, NLO-DM,
101044443). R.G. is supported by MIUR grants PRIN 2017FMJFMW and 2017L5W2PT,
and thanks Galileo Galilei Institute for hospitality. C.G.C. is supported by a
Ramón y Cajal contract with Ref. RYC2020-029248-I, the Spanish National Grant
PID2022-137268NA-C55 and Generalitat Valenciana through the grant
CIPROM/22/69. The work of T.H. is supported by the Excellence of Science (EoS)
project No. 30820817 - be.h “The H boson gateway to physics beyond the
Standard Model”, by the IISN convention 4.4503.15, and by CERN where part of
this work has been done. The authors acknowledge the workshop on Self-
Interacting Dark Matter: Models, Simulations and Signals 2023, Pollica, for
hospitality.

##  Appendix A Approximate parametric solution to number evolution of DM
particles

The set of coupled Boltzmann equations for $N_{\chi}$ and $N_{2\chi}$, Eqs.
(2)-(3), has no closed-form analytical solution. However, we can construct an
approximate solution and better understand the underlying physics by
proceeding step by step as follows.

###  A.1 Case with no self-capture on DM bound-states

Without the $C_{2\chi}$ term, the equation for $N_{\chi}$ does not depend on
$N_{2\chi}$

$\displaystyle\frac{{\rm d}N_{\chi}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle C_{\star}-A_{\rm bsf}N^{2}_{\chi}+C_{\chi}N_{\chi}\,,$ (35)

which can be solved easily as follows.

####  BSF with capture only on nucleons: $C_{\chi}=C_{2\chi}=0$ :

Without any capture due to DM self-interactions, the above equation has the
solution

$N^{0}_{\chi}(t)=C_{\star}\tau_{0}\tanh\left(\frac{t}{\tau_{0}}\right)\quad\quad,\quad\quad
N^{0}_{2\chi}(t)={1\over
2}C_{\star}\left(t-\tau_{0}\tanh\left(\frac{t}{\tau_{0}}\right)\right)~{},$
(36)

with the time constant given by $\tau_{0}=(C_{\star}A_{\rm bsf})^{-1/2}$.
Without BSF, i.e. when $\tau_{0}\to\infty$, this solution gives a $N_{\chi}$
linearly increasing in time as expected,
$N^{0}_{\chi}(t\lesssim\tau_{0})=C_{\star}t$. With BSF, such an increase
occurs until there are enough DM particles for BSF to proceed. For
$t\gtrsim\tau_{0}$, an equilibrium between BSF and capture on nucleons is
reached so that $N_{\chi}$ saturates becoming constant:
$N_{\chi}=C_{\star}\tau_{0}=(C_{\star}/A_{\rm bsf})^{1/2}$. This value can be
obtained also directly from the Boltzmann equation, Eq. (2), whose right-hand
side vanishes when $C_{\star}-A_{\rm bsf}N_{\chi}^{2}=0$. As expected, the
larger $A_{\rm bsf}$, the sooner the saturation is reached and the smaller the
asymptotic value of $N_{\chi}$ is. Conversely, the larger $C_{\star}$, the
larger $N_{\chi}$ must be for the term associated with BSF to compensate that
of $C_{\star}$. Note, nevertheless, that the larger $C_{\star}$ and the larger
$A_{\rm bsf}$, the smaller the saturation time scale $\tau_{0}$. This stems
from the fact that the number of particles captured increases more rapidly,
and BSF starts to become important earlier. Consequently, BSF catches up with
the capture process sooner.

####  BSF with capture on nucleons and on free DM particles: $C_{\chi}\neq
0,\,C_{2\chi}=0$ :

Early in the evolution the contribution of $C_{2\chi}N_{2\chi}$ is negligible.
The analytical solution when $C_{2\chi}=0$ and $C_{\chi}\neq 0$ is

$\displaystyle N^{s}_{\chi}(t)$ $\displaystyle=$
$\displaystyle\frac{C_{\chi}}{2A_{\rm bsf}}+\frac{1}{A_{\rm
bsf}\tau_{s}}\tanh\left(\frac{t}{\tau_{s}}-{\rm
arctanh}\left(\frac{C_{\chi}\tau_{s}}{2}\right)\right)\,,$ (37) $\displaystyle
2N^{s}_{2\chi}(t)+N^{s}_{\chi}(t)$ $\displaystyle=$
$\displaystyle\left(C_{\star}+\frac{C_{\chi}^{2}}{2A_{\rm
bsf}}+\frac{C_{\chi}}{A_{\rm bsf}\tau_{s}}\right)t+\frac{C_{\chi}}{A_{\rm
bsf}}\log\left(\frac{2-C_{\chi}\tau_{s}+e^{-\frac{2t}{\tau_{s}}}\left(2+C_{\chi}\tau_{s}\right)}{4}\right)\,,$
(38)

with $\tau_{s}=(C_{\star}A_{\rm bsf}+C^{2}_{\chi}/4)^{-1/2}$. As it must be,
for $C_{\chi}=0$ these equations reduce to Eq. (36) and $\tau_{s}$ reduces to
$\tau_{0}$.

Here again, the BSF term equilibrates with capture in Eq. (2), as it is
negative and quadratic in $N_{\chi}$, whereas the $C_{\star}$ ($C_{\chi}$)
term is constant (linear) in it. Thus, $N_{\chi}$ saturates when the right-
hand side of Eq. (2) vanishes, $C_{\star}-A_{\rm
bsf}N_{\chi}^{2}+C_{\chi}N_{\chi}=0$, which gives $N^{s,\rm
eq}_{\chi}=C_{\chi}/(2A_{\rm bsf})+1/(\tau_{s}A_{\rm
bsf})=C_{\star}/(\tau_{s}^{-1}-C_{\chi}/2)$. The equilibrium time scale
$\tau_{s}$ is smaller than that without the $C_{\chi}$ term, $\tau_{0}$,
because the $C_{\chi}$ term increases the capture, so that BSF becomes
important earlier. Note that before equilibrium is reached the solution of Eq.
(37) is exponential. However the effect of this exponential is very limited as
the equilibrium is reached soon, as can be seen in Fig. (3), for various
examples of parameter sets.888Actually the solution of Eq. (37) can be
rewritten in exponential form as
$2C_{\star}\tau_{s}(e^{2t/\tau_{s}}-1)/(A+Be^{2t/\tau_{s}})$ with
$A=-2-C_{\chi}\tau_{s}$ and $B=-2+C_{\chi}\tau_{S}$. In practice if
$C_{\star}A_{\rm bsf}<C^{2}_{\chi}/4$ (as will be the case for our scenario),
then $\tau_{s}\simeq 2/C_{\chi}$ which gives $B<A$ so that $N_{\chi}\simeq
2C_{\star}\tau_{s}(e^{2t/\tau_{s}}-1)/A$, which is exponentially growing until
it reaches the equilibrium plateau when $Be^{2t/\tau_{s}}$ becomes larger than
$A$. The amount of exponential grows is limited because in practice $B$ is not
much smaller than $A$. Note also that the total number of particles grows
linearly, Eq. (38), except for a logarithmic correction, particularly for
$t\gg\tau_{s}$. More importantly, for realistic values of $C_{\star}$ and
$C_{\chi}$, the rate associated with the saturation value always lies well
below the geometric rate $C^{g}_{\chi}$ of thermalized DM particles in the
Sun.

###  A.2 Case with self-capture on DM bound-states

The appearance of a non-vanishing $C_{2\chi}$ term drastically changes the
physics in several ways. First of all, the $C_{2\chi}$ term implies that the
number of free DM particles in the Sun depends on the number of DM bound-
states so that the Boltzmann equations for $N_{\chi}$ and $N_{2\chi}$ are
coupled. Differentiating Eq. (2) with respect to time and using the equation
for the bound-state, Eq. (3), results in the following second-order
differential equation,

$\displaystyle\frac{{\rm d}^{2}N_{\chi}}{{\rm d}^{2}t}$ $\displaystyle=$
$\displaystyle\left(-2A_{\rm bsf}N_{\chi}+C_{\chi}\right)\frac{{\rm
d}N_{\chi}}{{\rm d}t}+\frac{1}{2}C_{2\chi}A_{\rm bsf}N^{2}_{\chi}\,,$ (39)

with boundary conditions: $\frac{{\rm d}N_{\chi}(t=0)}{{\rm d}t}=C_{\star}$
and $N_{\chi}(t=0)=0$. Formally, Eq. (39) is of the form
$N^{\prime\prime}=f(N)N^{\prime}+g(N)$, known as a Liénard equation, with no
known general analytical solutions.

Second, the self-capture associated with $C_{2\chi}$ is larger than that
induced by $C_{\chi}$. This and the fact that $N_{2\chi}$ never saturates –as
Eq. (3) shows– imply that the term $C_{2\chi}N_{2\chi}$ in Eq. (2) becomes
much larger than $C_{\chi}N_{\chi}$. In fact, $N_{\chi}$ saturates to values
much larger than in the case of $C_{\chi}\neq 0,C_{2\chi}=0$, because the
right-hand side of Eq. (2) vanishes for $A_{\rm bsf}N_{\chi}^{2}\simeq
C_{2\chi}N_{2\chi}\gg C_{\chi}N_{\chi}$. Thus, unlike in the case of
$C_{\chi}\neq 0,\,C_{2\chi}=0$, where the exponential growth induced by the
$C_{\chi}$ term does not last long, the $C_{2\chi}$ term induces an
exponential growth for $N_{\chi}$ which is both faster (since the argument of
the exponential will be proportional to $C_{2\chi}$ rather than to $C_{\chi}$)
and much larger. Therefore, $N_{\chi}$ can reach in this way much higher
values, large enough for the capture rate to reach the geometric limit within
the DMBS thermal sphere. Approximate solutions can be obtained as follows.

####  BSF with capture on nucleons and on DM bound-states:
$C_{\chi}=0,\,C_{2\chi}\neq 0$ :

For the parameters of interest in this work, $C_{2\chi}\ll 100/\tau_{0}$, we
numerically find that the following exponential ansatz provides a good
approximation from $t=\tau_{0}$ to the point when $C_{2\chi}N_{2\chi}$
saturates the geometric capture rate. That is, in practice

$N_{\chi}(\tau_{0}\lesssim
t\lesssim\tau_{g})=N_{\chi}^{0}(\tau_{0})\exp\Bigg{[}{C_{2\chi}\over
4}(t-\tau_{0})\Bigg{]}~{}.$ (40)

This is also suggested by the quasi-static solution of $dN_{\chi}/dt\sim 0$ in
Eqs. (2) and (3), where

$\frac{dN_{2\chi}}{dt}\simeq\frac{C_{\star}}{2}+\frac{C_{2\chi}}{2}N_{2\chi}\,~{}~{}\text{and~{}~{}~{}}N_{\chi}\simeq\left(\frac{C_{\star}+C_{2\chi}N_{2\chi}}{A_{\rm
bsf}}\right)^{1/2}\,.$ (41)

For the opposite case, $C_{2\chi}\geq 100/\tau_{0}$, the divergence happens
even faster.

Therefore, capture on bound-states exponentially increases the number of DM
particles within the Solar lifetime if the condition, $t_{\odot}-\tau_{0}\gg
C_{2\chi}^{-1}$ is satisfied. The fact that the exponential growth starts when
$N_{2\chi}$ becomes of order $N_{\chi}$ can also be seen in the numerical
examples shown in Fig. 3.

####  BSF in the full general case: $C_{\chi}\neq 0,\,C_{2\chi}\neq 0$ :

For the reasons explained above, switching on the $C_{\chi}$ term –in addition
to the $C_{2\chi}$ term– does not drastically change the result. It induces an
additional moderate exponential growth that makes the contribution of
$C_{2\chi}$ important slightly earlier (i.e. around $t=\tau_{s}$ rather than
at $\tau_{0}$), see Fig. 3. An approximate solution to this general case is
obtained in the same way as Eq. (40), by matching Eq. (37) at $\tau_{s}$
rather than at $\tau_{0}$. Interestingly, since $\tau_{s}<2/C_{\chi}\sim
4/C_{2\chi}$, we have $\tau_{s}C_{2\chi}\lesssim 4$. That is, the condition
for the validity of Eq. (40), $C_{2\chi}\ll 100/\tau_{0}$, in which $\tau_{0}$
is now replaced by $\tau_{s}$, is automatically satisfied. This suggests

$N_{\chi}(\tau_{s}\lesssim
t\lesssim\tau_{g})=N^{s}_{\chi}(\tau_{s})\exp\Bigg{[}{{C_{2\chi}\over
4}\left(t-\tau_{s}\right)}\Bigg{]}~{},$ (42)

where $N^{s}_{\chi}(\tau_{s})$ is well approximated by Eq. (37).

This exponential growth lasts until the geometric capture rate within the DMBS
thermal sphere is saturated. Quickly after $t=\tau_{g}$, $N_{\chi}$ stops
increasing when the BSF term (quadratic in $N_{\chi}$) compensates the
constant capture rate.999Once the capture rate reaches the DMBS sphere
geometric rate, the time it takes for $N_{\chi}$ to become constant is small,
$\Delta t\simeq 1/\sqrt{A_{\rm bsf}C^{g}_{2\chi}}$. The quasi-static
equilibrium solution, obtained from $dN_{\chi}/dt\simeq
C_{\star}+C_{2\chi}^{g}-A_{bsf}N_{\chi}^{2}\simeq 0$, yields the final
particle number101010Here we neglect the subleading capture contribution, from
the $C_{\chi}$ term, within the DM thermal sphere outside the DMBS thermal
sphere. To be accurate the $C_{\star}$ term in Eq. (43) has to be taken into
account only outside the DMBS thermal sphere as the $C_{2\chi}^{g}$ term
already accounts for the maximal total capture rate within this sphere, but
this concerns a negligible effect.

$N_{\chi}(t>\tau_{g})=\left(\frac{C^{g}_{2\chi}+C_{\star}}{A_{\rm
bsf}}\right)^{1/2}\,.$ (43)

The associated time $\tau_{g}$, where such an equilibrium is reached, is
determined as

$\tau_{g}=\tau_{s}+\frac{2}{C_{2\chi}}\log\left(\frac{C^{g}_{2\chi}+C_{\star}}{A_{\rm
bsf}(N^{s}_{\chi}(\tau_{s}))^{2}}\right)\,.$ (44)

The condition that $\tau_{g}\leq t_{\odot}$, so that the DM number can be
maximized at present, is satisfied if

$\,\frac{C_{2\chi}}{2}\left(t_{\odot}-\tau_{s}\right)\geq\log\left(\frac{C^{g}_{2\chi}+C_{\star}}{A_{\rm
bsf}(N^{s}_{\chi}(\tau_{s}))^{2}}\right)=\begin{cases}\log\left(\frac{C^{g}_{2\chi}/C_{\star}+1}{\tanh[1]^{2}}\right)&\text{if~{}}C_{\chi}^{2}\leq
A_{\rm bsf}C_{\star}\\\
\log\left(\frac{C^{g}_{2\chi}/C_{\star}+1}{(e^{2}-1)^{2}}\,\frac{C_{\chi}^{2}}{A_{\rm
bsf}C_{\star}}\right)&\text{if~{}}C_{\chi}^{2}\gg A_{\rm
bsf}C_{\star}\end{cases}\,.$ (45)

The latter case has $\tau_{0}\gg\tau_{s}\simeq 2/C_{\chi}\sim 4/C_{2\chi}$, so
the saturation should happen within a few $\tau_{s}$.

We can now analytically match the above two solutions at times $\sim\tau_{s}$
and $\tau_{g}$. The final approximate solution for species $i$ has the
parametric form

$N_{i}(t)~{}=~{}N_{i}^{t<\tau_{s}}\Theta(\tau_{s}-t)+N_{i}^{\tau_{s}<t<\tau_{g}}\Theta(t-\tau_{s})\Theta(\tau_{g}-t)+N_{i}^{t>\tau_{g}}\Theta(t-\tau_{g})\,,$
(46)

where $N_{i}^{t<\tau_{s}}$, $N_{i}^{\tau_{s}<t<\tau_{g}}$ and
$N_{i}^{t>\tau_{g}}$ are given by Eqs. (37), (42) and (43), respectively. We
have checked that these analytical expressions agree well with our numerical
results presented in the various figures from solving the Boltzmann equations.

##  Appendix B (Non-)perturbative treatment of elastic DM scattering

Suppose that DM scatters off a target with a mass $m_{T}$ with a negligible
form factor, the cross section differential in the recoil energy is

$\displaystyle\frac{{\rm d}\sigma}{{\rm
d}E_{R}}=\frac{2\pi\,m_{T}}{k^{2}}\big{|}{\cal M}(k,\theta)\big{|}^{2}\,,$
with $\displaystyle\cos\theta=1-\frac{m_{T}E_{R}}{k^{2}}\,.$ (47)

Here $\theta$ and $k=\mu v$ are respectively the scattering angle and the
incoming momentum in the center of mass frame, while the reduced mass of the
DM-target system $\mu=m_{\chi}m_{T}/(m_{\chi}+m_{T})$. We aim to calculate
this when the scattering is triggered by the exchange of a scalar mediator,
particularly when non-perturbative effects cannot be disregarded. In the non-
relativistic limit, this is described by a free DM particle scattered by the
corresponding Yukawa potential $V(r)=-\alpha^{\prime}e^{-m_{\phi}r}/r$, with
$\alpha^{\prime}=g_{s}y_{N}\sin\theta_{\phi}/4\pi,\alpha$ or $2\alpha$ when
target is a nucleon, a DM particle or a DMBS, respectively.

The amplitude of such quantum scattering can be obtained from a partial-wave
expansion

$\displaystyle{{\cal M}}(k,\theta)$
$\displaystyle=\frac{1}{2ik}\sum^{\infty}_{\ell=0}(2l+1)P_{\ell}\left(\cos\theta\right)\left(e^{2i\delta_{\ell}}-1\right)\,,$
(48)

where $\ell$ is the orbital angular momentum, $P_{\ell}(\cos\theta)$ is the
corresponding Legendre polynomial and $\delta_{\ell}$ is the phase-shift. In
general, the phase-shift must be obtained by solving the radial part of the
Schrödinger equation describing the collision, which is equivalent to solving
Chu  _et al._ (2020)

$\displaystyle\delta_{\ell,k}^{\prime}(r)=-2k\mu\,r^{2}V(r){\text{
Re}}\left[e^{i\delta_{\ell,k}(r)}h^{(1)}_{\ell}(kr)\right]^{2}\,\,\text{~{}with~{}}\delta_{\ell,k}(0)=0\,\text{~{}and~{}~{}}\delta_{\ell,k}(r)|_{r\to\infty}\to\delta_{\ell}\,.$
(49)

Here, $h^{(1)}_{\ell}=j_{\ell}+i\,n_{\ell}$ is the spherical Hankel function
of the first kind. Nevertheless, solving Eq. (49) for all values of $\ell$ is
impractical, and thus simplifications are necessary for different parameter
regimes, as shown in left panel of Fig. 7. This will be discussed below for
integrated and differential DM-DM(BS) scattering, separately. The Born
approximation is always justified for the very weak scattering between DM and
nucleon.

###  B.1 Integrated DM self-scattering in DM halos

Sufficiently large DM self-scattering can alter the evolution of halos, which
gives rise to constraints on the corresponding scattering rates. As argued in
Ref. Yang and Yu (2022), the gravothermal evolution of such halos is best
characterized by the so-called viscosity cross section,
$\sigma_{V}\equiv{3\over
2}\int^{1}_{-1}d\cos\theta\,\sin^{2}\theta{d\sigma\over
d\cos\theta}$.111111Following Ref. Yang and Yu (2022), we introduce a
prefactor of $3/2$ to make $\sigma_{V}=\sigma$ for isotropic scatterings. The
parameter regimes where the Born and semi-classical approximations apply are
given in left panel of Fig. 7.

For $\alpha\lesssim v$ or $m_{\chi}\lesssim m_{\phi}/\alpha$, the first-order
Born approximation is justified Landau and Lifshits (1991). Note that
$\mu=m_{\chi}/2$ and $\alpha^{\prime}=\alpha$ in this case. Moreover,
accounting for the indistinguishability of the DM particles in $\chi-\chi$
scattering, we find

$\sigma_{V}=\frac{6\pi\beta^{2}}{m_{\phi}^{2}}\left(\frac{\left(\frac{m_{\chi}v}{m_{\phi}}\right)^{4}+5\left(\frac{m_{\chi}v}{m_{\phi}}\right)^{2}+5}{\left(\frac{m_{\chi}v}{m_{\phi}}\right)^{4}+2\left(\frac{m_{\chi}v}{m_{\phi}}\right)^{2}}\log\left[\left(\frac{m_{\chi}v}{m_{\phi}}\right)^{2}+1\right]-\frac{5}{2}\right)\,.$
(50)

Where we have introduced $\beta=2m_{\phi}\alpha/(m_{\chi}v^{2})$ for later
convenience. Note that if one adopts the modified transfer cross section,
defined as
$\sigma_{|T|}\equiv\int^{1}_{-1}d\cos\theta\,(1-|\cos\theta|){d\sigma\over
d\cos\theta}$, approximately $\sigma_{V}\simeq 2\sigma_{|T|}$ if
$m_{\chi}v/m_{\phi}\ll 1$ as well as $\sigma_{V}\simeq 3\sigma_{|T|}$ if
$m_{\chi}v/m_{\phi}\gg 1$.

Outside the Born regime, non-perturbative effects are non-negligible,
especially in the resonant regime ($\alpha m_{\chi}/m_{\phi}\gtrsim 1$ and
$m_{\chi}v/m_{\phi}\lesssim 1$). In spite of this, in this parameter region
the scattering is generally dominated by the $s$-wave, $\ell=0$ Tulin  _et
al._ (2013), for which we can adopt the analytical results obtained by solving
for the phase shift, $\delta_{0}$ of $\ell=0$, under the assumption that the
Yukawa potential can be approximated by a Hulthén potential Cassel (2010). For
identical fermions, after adding the symmetric factor $1/2$ for final states,
this leads to

$\sigma_{V}={4\pi\over 2}{d\sigma\over d\Omega}\bigg{|}_{\ell=0}={2\pi\over
k^{2}}\Bigg{|}e^{i\delta_{0}}\Bigl{(}P_{0}(\cos\theta)+P_{0}(-\cos\theta)\Bigr{)}\sin\delta_{0}\Bigg{|}^{2}={8\pi\over
k^{2}}\Bigg{|}e^{i\delta_{0}}\sin\delta_{0}\Bigg{|}^{2}\,,$ (51)

where the scattering is isotropic, $k=m_{\chi}v/2$, and $P_{0}$ is the $0$th-
order Legendre polynomial. Note that this cross section only includes an even
spatial wave-function and an anti-symmetric spin wave-function. We take the
probability of such anti-parallel spin alignment to be 1/4. That is, its
corresponding cross section for identical fermion is reduced by half with
respect to that of non-identical fermions. The characteristic peaks of this
regime are visible in Figs. 2 and 5.

Finally, with the resulting expressions of the viscosity cross section for all
regimes, we impose bounds on DM self-interaction at both cluster and dwarf
scales, obtaining the observationally allowed parameter region shown in Fig.
5, see also section IV.

Figure 7: _Left_ : Dark matter self-scattering regimes. For self-capture of
dark matter by the Sun, the relevant parameters roughly lie within the dotted
line, and correspond to $v\sim 4\cdot 10^{-3}$, $m_{\chi}\gtrsim 5$ GeV,
$10^{-4}$ GeV $\lesssim m_{\phi}\lesssim 10^{-2}$ GeV, and
$0.02\lesssim\alpha<0.2$. See text for details. _Right_ : Curves above which
the captured DM particles thermalize with the nucleons in the Sun within
$t_{\odot}$.

###  B.2 Differential DM-DM(BS) scattering for DM capture

As explained in the main text and above, while the impact of DM self-
interactions in galaxies and galaxy clusters can be characterized by the
integrated transfer or viscosity total cross sections, for DM self-capture in
the Sun, it is crucial to obtain the differential scattering cross section. In
the latter case, left panel of Fig. 7 (dotted contour) illustrates the
parameter region of interest for the DM capture via scatterings off free and
bound DM particles already captured in the Sun.

The relevant parameters largely lie in the semi-classical regime, where $k\gg
m^{-1}_{\phi}$. That is, the de Broglie wavelength is much larger than the
range of the Yukawa potential and the scattering amplitude can be estimated
semi-classically. More precisely, the resulting sum in Eq. (48) can be
performed using the stationary phase approximation Landau and Lifshits (1991),
which, up to an inconsequential global phase, gives

$\displaystyle{\cal
M}=\frac{1}{k}\left(\frac{\ell_{0}(\theta)}{\sin\theta}\Big{|}\frac{{\rm
d}\ell_{0}(\theta)}{{\rm d}\theta}\Big{|}\right)^{\frac{1}{2}}\,,$ where
$\displaystyle\int^{\infty}_{r_{0}}\frac{\ell_{0}(\theta)+\frac{1}{2}}{r^{2}\sqrt{k^{2}-\left(\frac{\ell_{0}(\theta)+\frac{1}{2}}{r}\right)^{2}-2\mu
V(r)}}{\rm d}r=\frac{\pi+\theta}{2}\,.$ (52)

Here $r_{0}$ is the closest distance from the potential centre, and
$\ell_{0}(\theta)$ is the angular-momentum minimizing the phase in Eq. (48)
for a given scattering angle $\theta$. All this applies as long as
$\ell_{0}(\theta)\gg 1$, which justifies neglecting $1/2$ in Eq. (52). To see
how Eq. (52) resembles the classical equations of motion, we define the impact
parameter as $\rho\equiv\ell_{0}(\theta)/k$, which together with Eq. (47),
leads to

$\displaystyle\frac{{\rm d}\sigma}{{\rm d}\cos\theta}=2\pi\rho\frac{{\rm
d}\rho}{{\rm d}\cos\theta}\,,$ where
$\displaystyle\int^{\infty}_{r_{0}}\frac{\rho\,{\rm
d}r}{r^{2}\sqrt{1-\frac{\rho^{2}}{r^{2}}+2\beta\left(\frac{e^{-m_{\phi}r}}{m_{\phi}r}\right)}}=\frac{\pi+\theta}{2}\,,$
(53)

with $\beta\equiv\mu m_{\phi}\alpha^{\prime}/k^{2}$. The integral in Eq. (53)
gives the scattering angle according to classical mechanics. In this work, we
follow the numerical approach presented in Ref. Khrapak  _et al._ (2004, 2003)
to solve Eq. (53) and obtain the scattering rates.

As shown by Ref. Khrapak  _et al._ (2004), the calculation can be further
simplified for $\beta>13.2$, where the angular dependence of ${{\rm
d}\sigma}/{{\rm d}\cos\theta}$ becomes insensitive to $\beta$. Performing a
numerical fit to the cross sections of non-identical scattering and defining
${{\rm d}\sigma}/{{\rm
d}\cos\theta}\equiv\pi\rho^{2}_{*}(\beta,m_{\phi}){\mathcal{N}}(\cos\theta)$,
we find

$\mathcal{N}(c_{\theta})\simeq
0.35+0.05c_{\theta}+2.09c_{\theta}^{3}+3.25c_{\theta}^{4}-5.51c_{\theta}^{5}-8.53c_{\theta}^{6}+5.06c_{\theta}^{7}+7.22c_{\theta}^{8}\,.$
(54)

Note that the normalization,
$\int^{1}_{-1}dc_{\theta}(1-c_{\theta})\mathcal{N}(c_{\theta})\simeq 1$, is
chosen in a way we can directly adopt the overall prefactor from Khrapak  _et
al._ (2003) as

$\rho^{2}_{\ast}(\beta,m_{\phi})=\begin{cases}\frac{1}{m^{2}_{\phi}}\left(\log\beta+1-\frac{1}{2\,\log\beta}\right)^{2}&{\rm
for}\quad\quad\beta>1000~{},\\\
\frac{8}{m^{2}_{\phi}}\left(\frac{\beta^{2}}{1+1.5\,\beta^{1.65}}\right)&{\rm
for}\quad\quad 1000>\beta>13.2~{}.\\\ \end{cases}$ (55)

The amplitudes at the vicinity of $\cos\theta=1$, which may approach infinity,
are not included in the numerical fit above, given that small-angle
scatterings barely contribute to the self-capture rate as a result of
$\theta_{\rm min}$ in Eqs. (15) and (16).

Besides, under the assumption that the interference is subleading in the
(semi-)classical regime Landau and Lifshitz (1976), the previous approach can
also be applied to the scattering of identical particles by symmetrizing the
differential scattering amplitude $\theta\to\pi-\theta$, as

$\frac{d\sigma}{d\cos\theta}\bigg{|}_{\rm\chi\chi}=\frac{d\sigma}{d\cos\theta}+\frac{d\sigma}{d\cos\theta}\bigg{|}_{\cos\theta\to-\cos\theta}\,,$
(56)

where the range of the scattering angle $\theta$ reduces to $[0,\,\pi/2]$.

##  Appendix C DM thermalization

We summarize below the general formalism for DM thermalization with the Solar
medium, and present the results applicable to the scenario presented in the
main text. After being captured, DM particle orbits become smaller than the
size of the celestial body. For contact interactions, the time taken for the
orbits to shrink to its thermal radius is given by Kouvaris and Tinyakov
(2011)

$t_{\rm therm}=\frac{m_{\chi}}{\left(\sum_{N}\rho_{N}\sigma_{\chi
N}\right)}\sqrt{\frac{m_{\chi}}{2\,T}}~{}.$ (57)

For the Sun, the time scale is estimated via

$\displaystyle t^{\odot}_{2}$ $\displaystyle=$ $\displaystyle 1.5\,{\rm
yrs}\left(\frac{m_{\chi}}{10\,\rm GeV}\right)^{3/2}\left(\frac{2.2\,\rm
KeV}{T_{c}}\right)^{1/2}\left(\frac{150\,\rm
g/cm^{3}}{\rho_{c}}\right)\left(\frac{10^{-40}\,\rm cm^{2}}{\sigma_{\chi
N}}\right)~{}.$ (58)

Thermalization time is dominated by the last stages of the process with the
typical energy and momentum transfer of the order of temperature Tinyakov  _et
al._ (2021). Formally, thermalization time is computed through the energy loss
rate. Consider the elastic scattering of a DM particle $\chi$ off a
distribution of target particles $T$:
$\chi(k)+T(p)\rightarrow\chi(k^{\prime})+T(p^{\prime})$. The interaction rate
per DM particle reads Bertoni  _et al._ (2013); Garani  _et al._ (2021):

$\displaystyle d\Gamma=2\frac{{\rm
d}^{3}k^{\prime}}{(2\pi)^{3}}S(q_{0},q)~{},$ $\displaystyle
S(q_{0},q)=\iint\frac{{\rm d}^{3}p^{\prime}{\rm
d}^{3}p}{(2\pi)^{6}2E_{p^{\prime}}2E_{k^{\prime}}2E_{p}2E_{k}}(2\pi)^{4}\delta^{4}\left(k+p-k^{\prime}-p^{\prime}\right)|\mathcal{M}|^{2}f(E_{p})\left(1-f(E_{p^{\prime}})\right)\,,$
(59)

where the second line above is the response function, and
$f(E_{p,p^{\prime}})$ being the Maxwell-Boltzmann distributions w.r.t. the
core temperature $T_{c}$. For non-degenerate medium, there exists little
Pauli-blocking from the final states, so $1-f(E_{p^{\prime}})\simeq 1$. In
turn, the rate of energy loss is given by

$\Phi=2\int\frac{{\rm
d}^{3}k^{\prime}}{(2\pi)^{3}}S(q_{0},q)\times(E_{i}-E_{f})\,,$ (60)

where $k^{\prime}$ is integrated from 0 to $k$.

Using this we can write down the time taken to thermalize with the celestial
body, i.e. to reach a final energy $E_{f}=3/2\,T_{\rm c}$ starting with an
initial kinetic energy $E_{i}=m_{\chi}v_{\rm esc}^{2}/2$:

$\displaystyle\tau_{\rm therm}=-\int^{E_{f}}_{E_{i}}\frac{{\rm
d}E_{i}}{\Phi}.$ (61)

Integrating Eq. (C) over $p^{\prime}$, the response function becomes

$\displaystyle S(q_{0},q)$ $\displaystyle=$ $\displaystyle\int\frac{{\rm
d}^{3}p}{(2\pi)^{2}2}\frac{|\mathcal{M}|^{2}}{16m_{N}^{2}m^{2}_{\chi}}\delta^{3}\delta\left(q_{0}-E_{p}+E_{p^{\prime}}\right)f(E_{p}),$
(62)

where $E_{k}=E_{k^{\prime}}=m_{\chi}$ and $E_{p}=E_{p^{\prime}}=m_{N}$ have
been taken in the denominator, while for the numerator we use
$E_{p}=m_{N}+p^{2}/2m_{N}$ and $E_{p}^{\prime}=m_{N}+p^{\prime\,2}/2m_{N}$.
Moreover, the $\delta^{3}$ enforces momentum conservation,
$\mathbf{q}=\mathbf{p^{\prime}}-\mathbf{p}$, and the energy delta-function is
recast in terms of angle between the incoming nuclei and the momentum transfer
as shown in Ref. Reddy  _et al._ (1998), yielding

$S(q_{0},q)\approx\int\frac{{\rm d}p\,{\rm d}\cos\theta}{(2\pi)^{2}2}2\pi
p^{2}\frac{|\mathcal{M}|^{2}}{16m^{2}_{N}m^{2}_{\chi}}e^{\bar{\mu}/T_{c}}e^{-p^{2}/2m_{N}T_{c}}\frac{m_{N}}{|p||q|}\delta\left(\cos\theta-\cos\theta_{0}\right)\Theta(p^{2}-p^{2}_{-})~{},$
(63)

along with

$\cos\theta_{0}=\frac{m_{N}}{|p||q|}\left(q_{0}-\frac{q^{2}}{2m_{N}}\right)\quad\quad{\rm
and}\quad\quad
p^{2}_{-}=\frac{m^{2}_{N}}{q^{2}}\left(q_{0}-\frac{q^{2}}{2m_{N}}\right)^{2}~{}.$
(64)

Integrating the above equation we get

$S(q_{0},q)\approx\frac{|\mathcal{M}|^{2}}{64\pi
m^{2}_{\chi}}\frac{T_{c}}{|q|}e^{\bar{\mu}/T_{c}}e^{-p^{2}_{-}/2m_{N}T_{c}}\,.$
(65)

Taking the Higgs-mixing model introduced in the main text, the squared
amplitude for scattering on nucleon in the non-relativistic limit reads

$|\mathcal{M}|^{2}=g^{2}_{s}\cos^{2}_{\theta_{\phi}}\sin^{2}_{\theta_{\phi}}\frac{m^{2}_{N}f^{2}_{N}}{v^{2}_{H}}\frac{m^{2}_{\chi}m^{2}_{N}}{\left(q_{0}^{2}-q^{2}-m^{2}_{\phi}\right)^{2}}~{},$
(66)

which, in the limit of small energy transfer, corresponds to DM-nuclei
scattering cross section of

$\sigma_{\chi
N}=g^{2}_{s}\cos^{2}_{\theta_{\phi}}\sin^{2}_{\theta_{\phi}}\frac{1}{16\pi(m_{\chi}+m_{N})^{2}}\frac{m^{2}_{N}f^{2}_{N}}{v^{2}_{H}}\frac{m^{2}_{\chi}m^{2}_{N}}{m^{4}_{\phi}}\,.$
(67)

The $k^{\prime}$ integral is recast in terms of momentum transfer and energy
transfer. The limits of integration are $0<q_{0}<E_{k}$ and
$q_{0}<q<2E_{k}-q_{0}$. Therefore the energy loss rate for a particle with
initial momenta $k$ is

$\frac{d\Gamma}{dq_{0}}\Delta
E=\int\frac{dq}{(2\pi)^{2}}q\frac{E_{k^{\prime}}}{E_{k}}S(q_{0},q)(E_{k}-E_{k^{\prime}})=\int\frac{dq}{(2\pi)^{2}}q\frac{E_{k}-q_{0}}{E_{k}}S(q_{0},q)q_{0}~{}.$
(68)

In the non-degenerate limit the above equation factorizes. Hence the total
rate is proportional to cross section on one target particle times their
number density Reddy  _et al._ (1998). The corresponding numerical results are
shown in the right panel of Fig. 7, showing that the thermalization can be
reached within the lifetime of the Sun, even with very tiny mixing angles.
Including the DM self-interactions would further reduce the time scale of
thermalization.

## References

  * March-Russell and West (2009) J. D. March-Russell and S. M. West, Phys. Lett. B 676, 133 (2009), arXiv:0812.0559 [astro-ph] .
  * Pospelov and Ritz (2009) M. Pospelov and A. Ritz, Phys. Lett. B 671, 391 (2009), arXiv:0810.1502 [hep-ph] .
  * Shepherd  _et al._ (2009) W. Shepherd, T. M. P. Tait, and G. Zaharijas, Phys. Rev. D 79, 055022 (2009), arXiv:0901.2125 [hep-ph] .
  * Laha and Braaten (2014) R. Laha and E. Braaten, Phys. Rev. D 89, 103510 (2014), arXiv:1311.6386 [hep-ph] .
  * von Harling and Petraki (2014) B. von Harling and K. Petraki, JCAP 12, 033 (2014), arXiv:1407.7874 [hep-ph] .
  * Petraki  _et al._ (2014) K. Petraki, L. Pearce, and A. Kusenko, JCAP 07, 039 (2014), arXiv:1403.1077 [hep-ph] .
  * Petraki  _et al._ (2015) K. Petraki, M. Postma, and M. Wiechers, JHEP 06, 128 (2015), arXiv:1505.00109 [hep-ph] .
  * Pearce  _et al._ (2015) L. Pearce, K. Petraki, and A. Kusenko, Phys. Rev. D 91, 083532 (2015), arXiv:1502.01755 [hep-ph] .
  * Laha (2015) R. Laha, Phys. Rev. D 92, 083509 (2015), arXiv:1505.02772 [hep-ph] .
  * An  _et al._ (2016a) H. An, B. Echenard, M. Pospelov, and Y. Zhang, Phys. Rev. Lett. 116, 151801 (2016a), arXiv:1510.05020 [hep-ph] .
  * An  _et al._ (2016b) H. An, M. B. Wise, and Y. Zhang, Phys. Rev. D 93, 115020 (2016b), arXiv:1604.01776 [hep-ph] .
  * An  _et al._ (2017) H. An, M. B. Wise, and Y. Zhang, Phys. Lett. B 773, 121 (2017), arXiv:1606.02305 [hep-ph] .
  * Bi  _et al._ (2017) X.-J. Bi, Z. Kang, P. Ko, J. Li, and T. Li, Phys. Rev. D 95, 043540 (2017), arXiv:1602.08816 [hep-ph] .
  * Kouvaris  _et al._ (2016) C. Kouvaris, K. Langæble, and N. G. Nielsen, JCAP 10, 012 (2016), arXiv:1607.00374 [hep-ph] .
  * Cirelli  _et al._ (2017) M. Cirelli, P. Panci, K. Petraki, F. Sala, and M. Taoso, JCAP 05, 036 (2017), arXiv:1612.07295 [hep-ph] .
  * Asadi  _et al._ (2017) P. Asadi, M. Baumgart, P. J. Fitzpatrick, E. Krupczak, and T. R. Slatyer, JCAP 02, 005 (2017), arXiv:1610.07617 [hep-ph] .
  * Petraki  _et al._ (2017) K. Petraki, M. Postma, and J. de Vries, JHEP 04, 077 (2017), arXiv:1611.01394 [hep-ph] .
  * Johnson  _et al._ (2016) E. Johnson, E. Braaten, and H. Zhang, PoS ICHEP2016, 210 (2016), arXiv:1611.06212 [hep-ph] .
  * Mitridate  _et al._ (2017) A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, JCAP 05, 006 (2017), arXiv:1702.01141 [hep-ph] .
  * Baldes  _et al._ (2018) I. Baldes, M. Cirelli, P. Panci, K. Petraki, F. Sala, and M. Taoso, SciPost Phys. 4, 041 (2018), arXiv:1712.07489 [hep-ph] .
  * Pearce and Kusenko (2013) L. Pearce and A. Kusenko, Phys. Rev. D 87, 123531 (2013), arXiv:1303.7294 [hep-ph] .
  * Mahbubani  _et al._ (2020) R. Mahbubani, M. Redi, and A. Tesi, Phys. Rev. D 101, 103037 (2020), arXiv:1908.00538 [hep-ph] .
  * Mahbubani  _et al._ (2021) R. Mahbubani, M. Redi, and A. Tesi, JCAP 02, 039 (2021), arXiv:2007.07231 [hep-ph] .
  * Chen  _et al._ (2015) J. Chen, Z.-L. Liang, Y.-L. Wu, and Y.-F. Zhou, JCAP 12, 021 (2015), arXiv:1505.04031 [hep-ph] .
  * Feng  _et al._ (2016) J. L. Feng, J. Smolinsky, and P. Tanedo, Phys. Rev. D 93, 115036 (2016), [Erratum: Phys.Rev.D 96, 099903 (2017)], arXiv:1602.01465 [hep-ph] .
  * Arina  _et al._ (2017) C. Arina, M. Backović, J. Heisig, and M. Lucente, Phys. Rev. D 96, 063010 (2017), arXiv:1703.08087 [astro-ph.HE] .
  * Leane  _et al._ (2017) R. K. Leane, K. C. Y. Ng, and J. F. Beacom, Phys. Rev. D 95, 123016 (2017), arXiv:1703.04629 [astro-ph.HE] .
  * Lucente  _et al._ (2017) M. Lucente, C. Arina, M. Backović, and J. Heisig, PoS EPS-HEP2017, 628 (2017), arXiv:1710.03947 [hep-ph] .
  * Niblaeus  _et al._ (2019) C. Niblaeus, A. Beniwal, and J. Edsjo, JCAP 11, 011 (2019), arXiv:1903.11363 [astro-ph.HE] .
  * Cuoco  _et al._ (2020) A. Cuoco, P. De La Torre Luque, F. Gargano, M. Gustafsson, F. Loparco, M. N. Mazziotta, and D. Serini, Phys. Rev. D 101, 022002 (2020), arXiv:1912.09373 [astro-ph.HE] .
  * Mazziotta  _et al._ (2020) M. N. Mazziotta, F. Loparco, D. Serini, A. Cuoco, P. De La Torre Luque, F. Gargano, and M. Gustafsson, Phys. Rev. D 102, 022003 (2020), arXiv:2006.04114 [astro-ph.HE] .
  * Dasgupta  _et al._ (2020) B. Dasgupta, A. Gupta, and A. Ray, JCAP 10, 023 (2020), arXiv:2006.10773 [hep-ph] .
  * Bell  _et al._ (2021a) N. F. Bell, J. B. Dent, and I. W. Sanderson, Phys. Rev. D 104, 023024 (2021a), arXiv:2103.16794 [hep-ph] .
  * Garani and Palomares-Ruiz (2022) R. Garani and S. Palomares-Ruiz, JCAP 05, 042 (2022), arXiv:2104.12757 [hep-ph] .
  * Abud  _et al._ (2021) A. A. Abud _et al._ (DUNE), JCAP 10, 065 (2021), arXiv:2107.09109 [hep-ex] .
  * Leane  _et al._ (2021) R. K. Leane, T. Linden, P. Mukhopadhyay, and N. Toro, Phys. Rev. D 103, 075030 (2021), arXiv:2101.12213 [astro-ph.HE] .
  * Abbasi  _et al._ (2022) R. Abbasi _et al._ (IceCube), Phys. Rev. D 105, 062004 (2022), arXiv:2111.09970 [astro-ph.HE] .
  * Bell  _et al._ (2021b) N. F. Bell, M. J. Dolan, and S. Robles, JCAP 11, 004 (2021b), arXiv:2107.04216 [hep-ph] .
  * Zakeri and Zhou (2022) M. Zakeri and Y.-F. Zhou, JCAP 04, 026 (2022), arXiv:2109.11662 [hep-ph] .
  * Bose  _et al._ (2022) D. Bose, T. N. Maity, and T. S. Ray, Phys. Rev. D 105, 123013 (2022), arXiv:2112.08286 [hep-ph] .
  * Maity  _et al._ (2023) T. N. Maity, A. K. Saha, S. Mondal, and R. Laha, (2023), arXiv:2308.12336 [hep-ph] .
  * Gresham and Zurek (2019) M. I. Gresham and K. M. Zurek, Phys. Rev. D 99, 083008 (2019), arXiv:1809.08254 [astro-ph.CO] .
  * Garani  _et al._ (2022) R. Garani, M. H. G. Tytgat, and J. Vandecasteele, Phys. Rev. D 106, 116003 (2022), arXiv:2207.06928 [hep-ph] .
  * Press and Spergel (1985) W. H. Press and D. N. Spergel, Astrophys. J. 296, 679 (1985).
  * Griest and Seckel (1987) K. Griest and D. Seckel, Nucl. Phys. B 283, 681 (1987), [Erratum: Nucl.Phys.B 296, 1034–1036 (1988)].
  * Wise and Zhang (2015) M. B. Wise and Y. Zhang, JHEP 02, 023 (2015), [Erratum: JHEP 10, 165 (2015)], arXiv:1411.1772 [hep-ph] .
  * Hardy  _et al._ (2015) E. Hardy, R. Lasenby, J. March-Russell, and S. M. West, JHEP 06, 011 (2015), arXiv:1411.3739 [hep-ph] .
  * Gresham  _et al._ (2018) M. I. Gresham, H. K. Lou, and K. M. Zurek, Phys. Rev. D 97, 036003 (2018), arXiv:1707.02316 [hep-ph] .
  * Gresham  _et al._ (2017) M. I. Gresham, H. K. Lou, and K. M. Zurek, Phys. Rev. D 96, 096012 (2017), arXiv:1707.02313 [hep-ph] .
  * Spergel and Press (1985) D. N. Spergel and W. H. Press, Astrophys. J. 294, 663 (1985).
  * Asplund  _et al._ (2009) M. Asplund, N. Grevesse, A. Sauval, and P. Scott, Ann. Rev. Astron. Astrophys. 47, 481 (2009), arXiv:0909.0948 [astro-ph.SR] .
  * Garani and Palomares-Ruiz (2017) R. Garani and S. Palomares-Ruiz, JCAP 05, 007 (2017), arXiv:1702.02768 [hep-ph] .
  * Bottino  _et al._ (2002) A. Bottino, G. Fiorentini, N. Fornengo, B. Ricci, S. Scopel, and F. L. Villante, Phys. Rev. D 66, 053005 (2002), arXiv:hep-ph/0206211 .
  * Lewin and Smith (1996) J. D. Lewin and P. F. Smith, Astropart. Phys. 6, 87 (1996).
  * Duda  _et al._ (2007) G. Duda, A. Kemper, and P. Gondolo, JCAP 04, 012 (2007), arXiv:hep-ph/0608035 .
  * Wise and Zhang (2014) M. B. Wise and Y. Zhang, Phys. Rev. D 90, 055030 (2014), [Erratum: Phys.Rev.D 91, 039907 (2015)], arXiv:1407.4121 [hep-ph] .
  * Tulin  _et al._ (2013) S. Tulin, H.-B. Yu, and K. M. Zurek, Phys. Rev. D 87, 115007 (2013), arXiv:1302.3898 [hep-ph] .
  * Colquhoun  _et al._ (2021) B. Colquhoun, S. Heeba, F. Kahlhoefer, L. Sagunski, and S. Tulin, Phys. Rev. D 103, 035006 (2021), arXiv:2011.04679 [hep-ph] .
  * Landau and Lifshits (1991) L. D. Landau and E. M. Lifshits,  _Quantum Mechanics_ , Course of Theoretical Physics, Vol. v.3 (Butterworth-Heinemann, Oxford, 1991).
  * Khrapak  _et al._ (2003) S. A. Khrapak, A. V. Ivlev, G. E. Morfill, and S. K. Zhdanov, Phys. Rev. Lett. 90, 225002 (2003).
  * Khrapak  _et al._ (2004) S. Khrapak, A. Ivlev, G. Morfill, S. Zhdanov, and H. Thomas, IEEE Transactions on Plasma Science 32, 555 (2004).
  * Zentner (2009) A. R. Zentner, Phys. Rev. D 80, 063501 (2009), arXiv:0907.3448 [astro-ph.HE] .
  * Ibarra  _et al._ (2012) A. Ibarra, S. Lopez Gehler, and M. Pato, JCAP 07, 043 (2012), arXiv:1205.0007 [hep-ph] .
  * Richard  _et al._ (2016) E. Richard _et al._ (Super-Kamiokande), Phys. Rev. D 94, 052001 (2016), arXiv:1510.08127 [hep-ex] .
  * Abe  _et al._ (2021) K. Abe _et al._ (Super-Kamiokande), Phys. Rev. D 104, 122002 (2021), arXiv:2109.11174 [astro-ph.HE] .
  * Abe  _et al._ (2022) S. Abe _et al._ (KamLAND), Astrophys. J. 925, 14 (2022), arXiv:2108.08527 [astro-ph.HE] .
  * Battistoni  _et al._ (2005) G. Battistoni, A. Ferrari, T. Montaruli, and P. R. Sala, Astropart. Phys. 23, 526 (2005).
  * Zhou and Beacom (2023) B. Zhou and J. F. Beacom, (2023), arXiv:2311.05675 [hep-ph] .
  * Suliga and Beacom (2023) A. M. Suliga and J. F. Beacom, Phys. Rev. D 108, 043035 (2023), arXiv:2306.11090 [hep-ph] .
  * Honda  _et al._ (2015) M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa, Phys. Rev. D 92, 023004 (2015), arXiv:1502.03916 [astro-ph.HE] .
  * Krnjaic (2016) G. Krnjaic, Phys. Rev. D 94, 073009 (2016), arXiv:1512.04119 [hep-ph] .
  * Winkler (2019) M. W. Winkler, Phys. Rev. D99, 015018 (2019), arXiv:1809.01876 [hep-ph] .
  * Matsumoto  _et al._ (2019) S. Matsumoto, Y.-L. S. Tsai, and P.-Y. Tseng, JHEP 07, 050 (2019), arXiv:1811.03292 [hep-ph] .
  * Harvey  _et al._ (2015) D. Harvey, R. Massey, T. Kitching, A. Taylor, and E. Tittley, Science 347, 1462 (2015), arXiv:1503.07675 [astro-ph.CO] .
  * Bondarenko  _et al._ (2018) K. Bondarenko, A. Boyarsky, T. Bringmann, and A. Sokolenko, JCAP 04, 049 (2018), arXiv:1712.06602 [astro-ph.CO] .
  * Harvey  _et al._ (2019) D. Harvey, A. Robertson, R. Massey, and I. G. McCarthy, Mon. Not. Roy. Astron. Soc. 488, 1572 (2019), arXiv:1812.06981 [astro-ph.CO] .
  * Sagunski  _et al._ (2021) L. Sagunski, S. Gad-Nasr, B. Colquhoun, A. Robertson, and S. Tulin, JCAP 01, 024 (2021), arXiv:2006.12515 [astro-ph.CO] .
  * Cross  _et al._ (2023) D. Cross _et al._ (DES), (2023), arXiv:2304.10128 [astro-ph.CO] .
  * Tulin and Yu (2018) S. Tulin and H.-B. Yu, Phys. Rept. 730, 1 (2018), arXiv:1705.02358 [hep-ph] .
  * Adhikari  _et al._ (2022) S. Adhikari _et al._ , (2022), arXiv:2207.10638 [astro-ph.CO] .
  * Vogelsberger  _et al._ (2012) M. Vogelsberger, J. Zavala, and A. Loeb, Mon. Not. Roy. Astron. Soc. 423, 3740 (2012), arXiv:1201.5892 [astro-ph.CO] .
  * Zavala  _et al._ (2013) J. Zavala, M. Vogelsberger, and M. G. Walker, Mon. Not. Roy. Astron. Soc. 431, L20 (2013), arXiv:1211.6426 [astro-ph.CO] .
  * Elbert  _et al._ (2015) O. D. Elbert, J. S. Bullock, S. Garrison-Kimmel, M. Rocha, J. Oñorbe, and A. H. G. Peter, Mon. Not. Roy. Astron. Soc. 453, 29 (2015), arXiv:1412.1477 [astro-ph.GA] .
  * Ren  _et al._ (2019) T. Ren, A. Kwa, M. Kaplinghat, and H.-B. Yu, Phys. Rev. X 9, 031020 (2019), arXiv:1808.05695 [astro-ph.GA] .
  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Athron  _et al._ (2019) P. Athron _et al._ (GAMBIT), Eur. Phys. J. C 79, 38 (2019), arXiv:1808.10465 [hep-ph] .
  * Aalbers  _et al._ (2023) J. Aalbers _et al._ (LZ), Phys. Rev. Lett. 131, 041002 (2023), arXiv:2207.03764 [hep-ex] .
  * Agashe  _et al._ (2014) K. Agashe, Y. Cui, L. Necib, and J. Thaler, JCAP 10, 062 (2014), arXiv:1405.7370 [hep-ph] .
  * Cirelli  _et al._ (2011) M. Cirelli, G. Corcella, A. Hektor, G. Hutsi, M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Strumia, JCAP 03, 051 (2011), [Erratum: JCAP 10, E01 (2012)], arXiv:1012.4515 [hep-ph] .
  * Argüelles  _et al._ (2021) C. A. Argüelles, A. Diaz, A. Kheirandish, A. Olivares-Del-Campo, I. Safa, and A. C. Vincent, Rev. Mod. Phys. 93, 035007 (2021), arXiv:1912.09486 [hep-ph] .
  * Elor  _et al._ (2016) G. Elor, N. L. Rodd, T. R. Slatyer, and W. Xue, JCAP 06, 024 (2016), arXiv:1511.08787 [hep-ph] .
  * Chu  _et al._ (2020) X. Chu, C. Garcia-Cely, and H. Murayama, JCAP 06, 043 (2020), arXiv:1908.06067 [hep-ph] .
  * Yang and Yu (2022) D. Yang and H.-B. Yu, JCAP 09, 077 (2022), arXiv:2205.03392 [astro-ph.CO] .
  * Cassel (2010) S. Cassel, J. Phys. G 37, 105009 (2010), arXiv:0903.5307 [hep-ph] .
  * Landau and Lifshitz (1976) L. Landau and E. Lifshitz,  _Mechanics: Volume 1_ , Course of theoretical physics (Elsevier Science, 1976).
  * Kouvaris and Tinyakov (2011) C. Kouvaris and P. Tinyakov, Phys. Rev. D 83, 083512 (2011), arXiv:1012.2039 [astro-ph.HE] .
  * Tinyakov  _et al._ (2021) P. Tinyakov, M. Pshirkov, and S. Popov, Universe 7, 401 (2021), arXiv:2110.12298 [astro-ph.HE] .
  * Bertoni  _et al._ (2013) B. Bertoni, A. E. Nelson, and S. Reddy, Phys. Rev. D 88, 123505 (2013), arXiv:1309.1721 [hep-ph] .
  * Garani  _et al._ (2021) R. Garani, A. Gupta, and N. Raj, Phys. Rev. D 103, 043019 (2021), arXiv:2009.10728 [hep-ph] .
  * Reddy  _et al._ (1998) S. Reddy, M. Prakash, and J. M. Lattimer, Phys. Rev. D 58, 013009 (1998), arXiv:astro-ph/9710115 .