Dark matter bound-state formation in the Sun

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

Symmetric DM accumulating in the Sun or other celestial bodies 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. .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. 1   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 1 Particle-antiparticle BSF inside the Sun for symmetric DM was considered in [14].
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 χ, with a Yukawa coupling g s to a lighter scalar particle ϕ L ⊃ −g s ϕ χχ .
To assume a scalar mediator is convenient because it induces an attractive interaction. 2We 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 ⋆ ) [48,49], three additional terms come into play in the Boltzmann equation determining the DM particle number N χ , due to DM-ϕ interactions.First, depending on the strength of self-interaction, DM particles could efficiently form χ-χ bound-states by emitting the particle mediating DM self-interactions.The rate of the latter is denoted by A bsf .Since the number of free particles is reduced by two per process, the χ is introduced in the equation that describes the evolution of N χ , where the second factor, 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 buildup of DM bound-states within the celestial body, whose number N 2χ is determined by a second Boltzmann equation that simply involves a + 1 2 A bsf N 2 χ term.Second, there is a term coming from capture of Galactic DM on already-captured DM particles, whose rate is denoted by C χ .Finally, there is a term from capture of Galactic DM on the formed bound-states, denoted by C 2χ .The presence of both C χ and C 2χ terms increase DM accretion.In summary, both populations evolve according to the following set of differential equations, The various rates C ⋆ , C χ and C 2χ 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.We further assume that the scattering between a galactic DM and a captured bound-state is also elastic, neglecting the possible formation of three-body bound-states or bound-state dissociation.Also, in these Boltzmann equations we do not write down explicitly additional terms coming from possible evaporation 2 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 [46,47].
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 evap ≈ 5 GeV, which is about 50% larger than that in the standard scenario [35].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) 3 where the respective geometric rates on the DM and DMBS thermal spheres are denoted by C g χ and C g 2χ .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 χ particles, under the assumption that there exist bottlenecks to form heavier bound-states, such as fast decays of (3χ) → (2χ) + χ and (4χ) → (2χ) + (2χ).See e.g.[50][51][52][53] 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. 4  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 B.
A. Determination of the various rates a.Thermal radius and geometric rates C g ⋆ , C g χ and C g 2χ : 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 3 These 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 χ ) and on nucleons (C g ⋆ ) are never reached. 4The presence of stable many-body bound-states, as well as self-capture via inelastic scattering, may reduce DM evaporation efficiency by quickly capturing free χ 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.
Sun and lie within different spheres of thermal radius r th .Noting that the mass of a bound-state is approximately twice that of a DM particle, these radii are obtained [48,54] by equating the average thermal energy, 3T ⊙ /2, to the gravitational potential energy per particle, 2πGr 2 th ρ ⊙ nm χ /3, as , 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 √ 2. We take the solar core temperature, T ⊙ , to be 2.2 keV, and the core mass density to be ρ ⊙ ∼150 g/cm 3 .For example, the thermal radius of 3 GeV DM particles is about one-tenth of the Solar radius, R ⊙ ≃ 6.9 × 10 5 km [55].In addition, we will assume throughout that the DM radial distribution is isothermal, with its temperature T ≈ T ⊙ [48,56].
Upon taking into account the relative motion of the Sun with respect to the galactic DM halo, the geometric capture rate on nucleons is [57] Here 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.[56].Throughout this work we take the local DM density ρ χ = 0.3 GeV/cm 3 .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 th /R ⊙ ) 2 .This results in geometric self-capture rates which scale as m −2 χ , b. Capture on nucleons rate C ⋆ : The capture rate of DM particles in the Sun has been well studied in the literature.It can be written as [26] C where f ⊙ (u χ ) 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 [55].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, ω(r) = u 2 χ + v e (r) 2 .The differential cross section, dσ i /dE R , encodes the energy dependence of elastic scattering between DM and nucleus i in the nonrelativistic limit.The recoil energy in the solar frame is given by E R = µ 2 χN /m N ×w 2 (1−cos θ CM ) , where µ χN = m χ m N /(m χ + 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 For definiteness we assume that the dark sector communicates with the SM particles through the Higgs portal, and that the light mediator ϕ mixes with the Higgs boson with a mixing angle θ ϕ .
The differential cross section in the Born approximation is then where 1) number (see e.g.[58,59]).
The DM capture rates on nucleons are shown in Fig. 1, corresponding to a Higgs mixing angle sin θ ϕ = 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 e , which is above MeV for m χ ≥ m evap .This leads to a t-channel enhancement, which in the m χ ≫ m N limit gives a capture cross section scaling as m −1 χ , leading to a capture rate scaling as m −2 χ .While, for mediator masses larger than ∼ 100 MeV and m χ ≲ m N , the capture rate in turn scales as m χ .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.c.Bound-state formation rate A bsf As stated above, we consider DM to be asymmetric and made of Dirac fermions (χ).The interaction between two (identical) DM particles is attractive when the mediator is a scalar (ϕ), with a classical Yukawa potential V (r) = − α r e −m ϕ r .The binding energy of the bound-state is E bind ≃ mχα 2  4 − αm ϕ , with the dark coupling α ≡ g 2 s /4π.From the merging of two identical fermions, the radiative BSF rate has been estimated in [60], yielding in the limit of α/v rel ≫ 1 and m ϕ → 0, where the DM relative velocity is v rel ≃ 2v χ,⊙ .The recoil velocity of the DMBS is α 2 /8 in this limit, which must be below the escape velocity from the Solar   center, resulting in the requirement α ≲ 0.18. 5 Finally, the total annihilation rate that enters the 5 Due 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.[53,60] for related discussions.
number evolution equation is given by The normalized DM radial distribution is which corresponds to an isothermal sphere, with a radial dependence set by the gravitational potential ϕ(r) = r 0 GM ⊙ (r ′ )/r ′ 2 dr ′ , with G the gravitational constant and M ⊙ (r ′ ) the mass inside a sphere of radius r ′ , and N χ (t) is the total population of DM particles at a given time t.
Once DM particles are thermalized its average velocity is v χ ≃ 2T /m χ .The smallness of this value considerably boosts the BSF in Eq. ( 12): 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 χ u 2 χ .The maximum scattering angle for scattering of identical particles is π/2.Therefore cos θ min = 1 − 2 u 2 χ ω 2 and cos θ 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 with cos We have used m 2χ ≈ 2m χ and set that the coupling to bound-state is 2α (as the vertex is of the scalar type).We assume that the form factor has the form of exp(E r /Q χ ).The typical size of the bound-state is set by the Bohr radius, hence To a good approximation F χ → 1 for a typical value of α = 0.1.The radial number density of DMBS (n 2χ ) is assumed to be isothermal, analogous to n χ .
Bound-state formation suggests that non-perturbative effects are non-negligible for the calculation of the DM self-scattering cross section (see e.g.[61,62]).This is indeed our case.For the pa- we are never in the Born regime of self-scattering, for which m χ < m ϕ /α or α ≲ v [63].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 χ v/m ϕ ≳ 1.In this case classical mechanics can be employed to estimate the scattering cross section. 6ore 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 (θ 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.[64,65] 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 χ for m χ ≫ m ϕ .For moderate values of DM masses the scaling is less steep.Note that the capture rate on boundstates C 2χ is approximately larger by a factor of two with respect to C χ , due to larger Yukawa coupling induced by bound-states and larger maximal allowed scattering angle that keep both particles gravitationally captured.

B. 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 χ (N 2χ ) 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 χ and C 2χ terms (black), or only the C χ term (blue) or only the C 2χ term (green).At early times the number of free particles N χ grows as C ⋆ t (i.e.red and black solid lines coincide).As N χ increases, the BSF starts to occur efficiently and the corresponding term (quadratic in N χ ) quickly catches up with the constant term associated with capture on nucleons, so that in absence of the C 2χ 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 χ one gets that the associated time scale is τ 0 ≡ (C ⋆ A bsf ) −1/2 if one drops the C χ term -which has a subleading effect (black solid curve)-or τ s ≡ (C ⋆ A bsf + C 2 χ /4) −1/2 if one includes it (green solid curve).For more details, see the end of this section and Appendix A.

N(t)
Another consequence of efficient BSF is that the number of bound-states N 2χ quickly catches up to the number of DM particles N χ , 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 boundstates, C 2χ N 2χ , takes over and leads to an exponential growth of N χ (red and blue solid curves in Fig. 3) and of N 2χ (dashed curves with same colors).The exponential growth become significant at t ≳ τ s +C −1 2χ , or at τ 0 +C −1 2χ , depending on whether one includes the C χ 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χ N 2χ reaches C g 2χ .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 χ N χ term as N χ ≪ N 2χ , thus the Boltzmann equation takes the form, 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 χ /dt ≃ 0 gives the maximum final value We denote as τ g the time when N χ freezes in such a way.An analytic approximation is given in Appendix A.
Note that, contrary to the C 2χ N 2χ term -which saturates the geometric rate C g 2χ -the average DM density within the thermal radius r χ th is smaller than the DMBS density within r 2χ th , so the self-capture term C χ N χ does not reach C g χ at t = τ g , and will never reach it after, as N χ stops increasing.As for N 2χ , it keeps increasing in time forever.For t > τ g , it increases linearly in time and the rate of BSF is half the capture rate, inducing a flux of mediators given by the same rate.
The reason why the C χ term has a subleading effect with respect to the effect of the C 2χ term stems from the fact that these two terms are very different in nature.If we switch off the C 2χ term (green curves), the C χ term also leads to an exponential grow, starting slightly before t = τ s , but it is much less important than the exponential growth from the C 2χ term.This is analogous to what happens for the symmetric DM case with BSF playing the role of annihilation if C 2χ = 0.Moreover, the argument of the exponential from the C 2χ term is larger than the one from the C χ term, because C 2χ is about a factor of two times C χ , see Fig. 2.
In presence of both C χ and C 2χ terms, the C χ term has only a moderate impact on the evolution of N χ and N 2χ .As a comparison of the full evolution (red curves in Fig. 3) and the C χ = 0 evolution (blue curves) shows, including C χ reduces the timescale when the exponential growth starts, from t ∼ τ 0 to t ∼ τ s .Thus, due to the C χ term, the exponential effect of the C 2χ term starts somewhat earlier, and the geometric rate within the bound-state thermal sphere is also reached somewhat earlier.The values of N χ and N 2χ are insensitive to the C χ term at t ≫ τ 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 τ g , at which this happens, can be easily shorter than the age of the Sun, as solved in Eqs.(A10) and (A11).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 bind , the differential flux at a terrestrial detector is where d ⊙ is the Earth-Sun distance (AU), with ϑ being the angular sensitivity of detector (as the apparent angular diameter of the Sun ϑ ⊙ ≈ 0.5 • ).Since the light mediator is boosted, the neutrino energy spectrum this decay leads to is not monochromatic but has a characteristic box shape [70] In particular, E − ≃ 0 and The differential flux of neutrinos at the detector (ignoring oscillations) is The resulting neutrino fluxes are shown in Fig. 4  green points [77].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 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 [72].Atmospheric low energy neutrino measurements from Super-K are shown by the green dots [71].Thin green line denotes the predicted atmospheric neutrino flux for 30 • sky [76].
to C g 2χ ∝ m −2 χ .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.
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.

A. 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.[78][79][80]) Here we do not specify the scalar potential V (ϕ, H), and simply assume that it induces a ϕ-H mixing angle θ ϕ (which can be achieved in various ways), so that, upon rotation to the mass basis, the following interactions are obtained: 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 ϕ.Thus one has the extra interactions In the physical ν ′ , N ′ and ϕ ′ , h ′ mass eigenstate basis, this leads to the following relevant Yukawa interactions for the light scalar eigenstate ϕ ′ , up to second order in the neutrino mixing angle, Thus, both the Y ϕ and Y ν interactions induce a decay of the light mediator into a pair of light neutrinos.In practice the decay induced by Y ϕ will be dominant because various constraints require Y ν sin θ ϕ ≪ Y ϕ sin θ ν .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 Likewise for the associated Yukawa coupling, the typical seesaw gives with Y ν ∼ 2m ν m N /v 2 h .If m N > m ϕ the light mediator cannot decay into N N or N ν and the νν is the only possible channel induced by the seesaw interactions.The decay width is In addition to this neutrino-mixing channel, ϕ can also decay, through the ϕ-H scalar mixing, into pairs of charged leptons or quarks (or SM bosons for large ϕ mass).As will be discussed below, this needs to be suppressed.

B. 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: 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., [81][82][83][84][85] σ SI m χ ≲ 0.5 cm 2 /g , for DM velocities around v ∼ 1000 km/s.On the other hand, the non-observation of gravothermal collapse in dwarf-sized halos sets an upper bound of about 100 cm 2 /g at small scales, for which we take v ∼ 25 km/s [86,87].Concretely, we derive the bounds taking σ 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.This case is disfavored by BBN observables, given its sizeable coupling to neutrinos.For m ϕ = 100 MeV, the dwarf galaxy scale bound is weaker and basically irrelevant for BSF.For α 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 [53].
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 Thus a fast enough decay is obtained provided that Y ϕ is not too small and m N is not too large.
BBN and perturbative couplings, . Thus low scale seesaw is favored along this scenario.This BBN bound can be relaxed if the ϕ 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).
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 [92,93] set an upper bound on g s sin θ ϕ , 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 θ ϕ by a factor of 2 -3, depending on the mediator mass [94].
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 ν ≈ E bind /2 is equivalent to symmetric DM annihilation with a mass of E bind /2, if the latter only makes up a fraction of the observed DM abundance.Therefore, existing bounds can be rescaled as follows, where (σ 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 This is similar to boosted DM case [95], which can vary mildly due to the uncertainty the Galactic J-factor [96].The corresponding upper bound on α is shown in the right panel of Fig. 6, obtained from re-scaling the indirect bounds holding in the symmetric DM case [97].As also illustrated in Fig. 5, this excludes values of m χ between 3.0 and 11.4 GeV for α = 0.15, and gives irrelevant bound for α = 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 Φ ⊙ ν /(GeV/m χ ) 2 ∼ O(1) cm −2 s −1 when the DM capture on DMBS is saturated by its geometric rate.Using Eq. ( 12) results in α 5/2 /m χ ≲ 10 −4 /GeV, which in turn requires m χ ≳ 5.6 GeV (87.1 GeV) for α = 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 [98], Here, f 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 χ as in the case of DM annihilation, but is given by the binding energy ≈ α 2 m χ /4, times the electromagnetic branching ratio of ϕ decay, Br ϕ→EM .For m ϕ ≪ E 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 Br 1/3 ϕ→EM α 7/3 ≲ 10 −5 m χ /GeV.Combining this condition with Eq. ( 27) gives the constraint Y ϕ θ 2 ν /θ ϕ ≳ 10 2.2 α 7/2 (GeV/m χ ) 3/2 for m ϕ below twice the muon mass.For typical seesaw values of the neutrino mixing angle, Eq. ( 26), this translates to the condition Mediator decay into right-handed neutrinos: If 2m N < m ϕ the dominant decay is not anymore into a pair of light neutrinos but into a pair of right-handed neutrinos (or into νN for m ϕ /2 < m N < m ϕ ), with their decay widths given by, Using the typical seesaw expectation of Eq. ( 26), one gets for 2m N < m ϕ and m ϕ /2 < m N < m ϕ , respectively.This can be compatible with the extra radiation constraint τ ϕ ≲ 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 ∼ 10 −3 , α ∼ 0.1, m ϕ < GeV and m χ ∼ 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 −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.
Here again, the BSF term equilibrates with capture in Eq. ( 2), as it is negative and quadratic in N χ , whereas the C ⋆ (C χ ) term is constant (linear) in it.Thus, N χ saturates when the right-hand side of Eq. ( 2) vanishes, ).The equilibrium time scale τ s is smaller than that without the C χ term, τ 0 , because the C χ term increases the capture, so that BSF becomes important earlier.Note that before equilibrium is reached the solution of Eq. (A3) 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. 8Note also that the total number of particles grows linearly, Eq. (A4), except for a logarithmic correction, particularly for t ≫ τ s .More importantly, for realistic values of C ⋆ and C χ , the rate associated with the saturation value always lies well below the geometric rate C g χ of thermalized DM particles in the Sun.

Case with self-capture on DM bound-states
The appearance of a non-vanishing C 2χ term drastically changes the physics in several ways.
First of all, the C 2χ 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 χ and N 2χ 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, with boundary conditions: = C ⋆ and N χ (t = 0) = 0. Formally, Eq. (A5) is of the form , known as a Liénard equation, with no known general analytical solutions.
Second, the self-capture associated with C 2χ is larger than that induced by C χ .This and the fact that N 2χ never saturates -as Eq. ( 3) shows-imply that the term C 2χ N 2χ in Eq. ( 2) becomes much larger than C χ N χ .In fact, N χ saturates to values much larger than in the case of C χ ̸ = 0, C 2χ = 0, because the right-hand side of Eq. ( 2) vanishes for A bsf N 2 χ ≃ C 2χ N 2χ ≫ C χ N χ .Thus, unlike in the case of C χ ̸ = 0, C 2χ = 0, where the exponential growth induced by the C χ term does not last 8 Actually the solution of Eq. (A3) can be rewritten in exponential form as 2C⋆τs(e 2t/τs − 1)/(A + Be 2t/τs ) with A = −2 − Cχτs and B = −2 + CχτS.In practice if C⋆A bsf < C 2 χ /4 (as will be the case for our scenario), then τs ≃ 2/Cχ which gives B < A so that Nχ ≃ 2C⋆τs(e 2t/τs − 1)/A, which is exponentially growing until it reaches the equilibrium plateau when Be 2t/τs becomes larger than A. The amount of exponential grows is limited because in practice B is not much smaller than A.
long, the C 2χ term induces an exponential growth for N χ which is both faster (since the argument of the exponential will be proportional to C 2χ rather than to C χ ) and much larger.Therefore, N χ 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 χ = 0, C 2χ ̸ = 0 : For the parameters of interest in this work, C 2χ ≪ 100/τ 0 , we numerically find that the following exponential ansatz provides a good approximation from t = τ 0 to the point when C 2χ N 2χ saturates the geometric capture rate.That is, in practice This is also suggested by the quasi-static solution of dN χ /dt ∼ 0 in Eqs. ( 2) and ( 3), where For the opposite case, C 2χ ≥ 100/τ 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 ⊙ − τ 0 ≫ C −1 2χ is satisfied.The fact that the exponential growth starts when N 2χ becomes of order N χ can also be seen in the numerical examples shown in Fig. 3.
BSF in the full general case: C χ ̸ = 0, C 2χ ̸ = 0 : For the reasons explained above, switching on the C χ term -in addition to the C 2χ term-does not drastically change the result.It induces an additional moderate exponential growth that makes the contribution of C 2χ important slightly earlier (i.e.around t = τ s rather than at τ 0 ), see Fig. 3.An approximate solution to this general case is obtained in the same way as Eq.(A6), by matching Eq. (A3) at τ s rather than at τ 0 .
Interestingly, since τ s < 2/C χ ∼ 4/C 2χ , we have τ s C 2χ ≲ 4. That is, the condition for the validity of Eq. (A6), C 2χ ≪ 100/τ 0 , in which τ 0 is now replaced by τ s , is automatically satisfied.This suggests where N s χ (τ s ) is well approximated by Eq. (A3).This exponential growth lasts until the geometric capture rate within the DMBS thermal sphere is saturated.Quickly after t = τ g , N χ stops increasing when the BSF term (quadratic in N χ ) compensates the constant capture rate. 9The quasi-static equilibrium solution, obtained from The associated time τ g , where such an equilibrium is reached, is determined as (A10) The condition that τ g ≤ t ⊙ , so that the DM number can be maximized at present, is satisfied if The latter case has τ 0 ≫ τ s ≃ 2/C χ ∼ 4/C 2χ , so the saturation should happen within a few τ s .
We can now analytically match the above two solutions at times ∼ τ s and τ g .The final approximate solution for species i has the parametric form where and N t>τg i are given by Eqs.(A3), (A8) and (A9), respectively.We have checked that these analytical expressions agree well with our numerical results presented in the various figures from solving the Boltzmann equations.by a free DM particle scattered by the corresponding Yukawa potential V (r) = −α ′ e −m ϕ r /r, with α ′ = g s y N sin θ ϕ /4π, α or 2α 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 where ℓ is the orbital angular momentum, P ℓ (cos θ) is the corresponding Legendre polynomial and δ ℓ 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 [99] ℓ (kr) Here, h 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. [100], the gravothermal evolution of such halos is best characterized by the so-called viscosity cross section, 11 The parameter regimes where the Born and semi-classical approximations apply are given in left panel of Fig. 7.
For α ≲ v or m χ ≲ m ϕ /α, the first-order Born approximation is justified [63].Note that µ = m χ /2 and α ′ = α in this case.Moreover, accounting for the indistinguishability of the DM particles in χ − χ scattering, we find Where we have introduced β = 2m ϕ α/(m χ v 2 ) for later convenience.Note that if one adopts the modified transfer cross section, defined as σ

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 ≫ m −1 ϕ .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. (B2) can be performed using the stationary phase approximation [63], which, up to an inconsequential global phase, gives , where with β ≡ µm ϕ α ′ /k 2 .The integral in Eq. (B7) gives the scattering angle according to classical mechanics.In this work, we follow the numerical approach presented in Ref. [64,65] to solve Eq. (B7) and obtain the scattering rates.
As shown by Ref. [65], the calculation can be further simplified for β > 13.In turn, the rate of energy loss is given by Φ = 2 d 3 k ′ (2π) 3 S(q 0 , q) × (E i − E f ) , (C4) where k ′ 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 c starting with an initial kinetic energy E i = m χ v 2 esc /2: Integrating Eq. (C3) over p ′ , the response function becomes where E k = E k ′ = m χ and E p = E p ′ = 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 = m N + p ′ 2 /2m N .Moreover, the δ 3 enforces momentum conservation, q = p ′ − p, and the energy delta-function is recast in terms of angle between the incoming nuclei and the momentum transfer as shown in Ref. [107], yielding S(q 0 , q) ≈ dp d cos θ (2π) Taking the Higgs-mixing model introduced in the main text, the squared amplitude for scattering on nucleon in the non-relativistic limit reads which, in the limit of small energy transfer, corresponds to DM-nuclei scattering cross section of The k ′ 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 E k S(q 0 , q)q 0 .(C12) 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 [107].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.

1 . 2 . 27 C
parametric solution to number evolution of DM particles 21 Case with no self-capture on DM bound-states 21 Case with self-capture on DM bound-states 22 B. (Non-)perturbative treatment of elastic DM scattering 24 1.Integrated DM self-scattering in DM halos 25 2. Differential DM-DM(BS) scattering for DM capture

FIG. 1 .
FIG. 1. Left: DM capture rate on nucleons for several mediator masses for a Higgs mixing angle sin θ ϕ = 10 −10 .The geometric capture on nucleons C g ⋆ and C g χ self-capture rates are shown by the dashed and dashed-dotted black curves.The C g 2χ self-capture rate line (not shown) is a factor 2 below the C g χ line.Right: The DM BSF rate inside the Sun, for several values of α.
be approximately independent of the DM mass in practice, as shown in right panel of Fig. 1. d. DM self capture rates C χ and C 2χ : Analogous to DM-nucleon capture rate, we estimate the self-capture rate per target DM particle in the Sun as follows 10 12 10 13 10 14 10 15 10 16 10 17 10 18 10 19 t [s] FIG. 3. Number evolution of DM and DM bound-states for two parameter sets, as a function of time.Here τ s and τ 0 are characteristic times explained in the main text (and exactly defined in Appendix A), t ⊙ 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.

7
InEq. (3) the C χ term (linear in N χ ) is quickly counter-balanced by formation of bound-states from the A bsf term (quadratic in N χ ): after a short period of exponential growth, N χ ∼ C ⋆ (e Cχt −1)/C χ , the green curve (C 2χ = 0), goes to a plateau from an equilibration of capture and BSF.In contrast, if C 2χ ̸ = 0, the capture rate on bound-state grows as N 2χ .Thus, due to capture on bound-state term C 2χ , if it were not for the geometric rate upper floor, both N χ and N 2χ would grow forever.
for a fixed Higgs mixing angle sin θ ϕ = 10 −12 , for α = 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 ν e neutrino background (DSNB) flux from Super-K runs III, IV and KamLAND are shown in blue, red and black points [71-73].Also shown are predicted atmospheric fluxes of ν e down to 10 MeV in thin green line, assuming 30 • angular resolution in the sky, adapted from [74-76].Current measurements of atmospheric neutrino fluxes are shown in FIG. 4. Left: Flux of light mediators emitted by the BSF process for various parameter choices.Right:

FIG. 5 .
FIG. 5. Constraints for two example values of the dark coupling α.Contour lines give the values of σ SI /m χ in dwarf-sized halos.Black (gray) hatched regions are excluded because σ SI /m χ ≥ 0.5 cm 2 /g at cluster scales (or 100 cm 2 /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.

Fig. 5 α
Fig.5shows the resulting values for this cross section as a function of m χ and m ϕ /m χ for two values of α.The cluster bound excludes χ masses below ∼ 10-1000 GeV depending on m ϕ /m χ and α.This figure also shows the region where m ϕ 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. σ SI /m χ ∈ [1, 100] cm 2 /g at typical DM velocities of the order of 25 km/s in dwarf galaxies[88-  91].In the m χ -α plane, Fig.6(right) shows the corresponding upper bound on α.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 ℓ = j ℓ + i n ℓ is the spherical Hankel function of the first kind.Nevertheless, solving Eq. (B3) for all values of ℓ 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.1.Integrated DM self-scattering in DM halos