Analytic Approximations for the Velocity Suppression of Dark Matter Capture

Dark matter (DM)—nonbaryonic, nonluminous matter that interacts predominantly gravitationally—has been a scientific puzzle since Zwicky coined the term Dunkele Materie (in translation "dark matter") in 1933 (Zwicky 1933). Today, there are numerous viable theories on the nature of DM. Some of the most notable include weakly interacting massive particles (WIMPs; see Roszkowski et al. 2018, and references therein), WIMPzillas (Kolb et al. 1999), and axions or axion-like particles (see Marsh 2016, and references therein), to name a few. For a while, massive astrophysical compact halo objects (MACHOs) were also popular candidates (for a review, see Evans & Belokurov 2005), and today a related version of this line of reasoning are primordial black holes as DM candidates (see Carr & Kuhnel 2020, and references therein). There are two distinct strategies for DM detection. One is direct detection, based on the interactions between DM and baryonic matter and the minute energy transferred to nuclei by collisions with the omnipresent sea of DM particles within our galaxy, through which the Earth and the Sun travel (for a recent review, see Schumann 2019). Any such DM signal should have a clear annual modulation, as predicted by Drukier et al. (1986). Intriguingly, such a signal has been detected by the DAMA experiment, starting in 1998 (Bernabei et al. 1998), and has persisted for more than two decades with an ever-increasing statistical significance (Bernabei et al. 2018). It is striking that none of the other direct-detection experiments have identified a similar signal. Recently, two experiments (ANAIS and COSINE) have been set up with the same detector technology (NaI) as the DAMA experiment, and, while preliminary, there is no indication of a statistically significant annual modulation in their data (Adhikari et al. 2019; Amaré et al. 2021). Other very sensitive DM direct-detection searches include the XENON1T (Aprile et al. 2012, 2018, 2019a, 2019b, 2020) experiment in Gran Sasso, Italy, the PICO experiment located at SNOLAB in Canada (Amole et al. 2019), and PandaX-II in China (Tan et al. 2016), among others. Despite the fact that these experiments have been running for some time, none of them have yet detected DM directly.

Rather than relying solely on direct detection, one can extract DM parameters for any model based on annihilation signals that could originate from DM-dense regions. This, in a nutshell, is the essence of indirect-detection techniques (for a review, see Feng et al. 2001). For instance, Freese et al. (2008), Iocco (2008), Ilie & Zhang (2019, 2020), and Ilie et al. (2020a, 2021) have discussed the impact of DM on a Population III star's luminosity. Most stars shine below the Eddington luminosity—that is, the brightest theoretically possible luminosity that preserves hydrostatic equilibrium. However, a star that has accreted enough DM may shine at the Eddington limit, and, as such, a limit on its mass can be placed if we know the DM–proton interaction cross section. Conversely, if Population III stars (zero-metallicity nuclear-powered first stars) are observed, their mere existence implies an upper bound on the cross section.

The role of astrophysical objects as potential DM laboratories has been recognized in the literature for a while. For instance, the pioneering work on DM capture (Faulkner & Gilliland 1985; Press & Spergel 1985; Spergel & Press 1985; Gould 1987, 1988) also deals with the potentially observable effects this phenomenon has on our Sun or the Earth. Via collisions with nuclei or electrons inside the dense environment of a compact object, such as a star, DM particles traversing it can be slowed down. Some of those will lose enough energy to become gravitationally trapped, and therefore captured, by the object. Subsequently, they will sink in toward the center of the star, where DM annihilations can produce energy (or secondary particles) that can have observable effects.

For dense capturing stars, and/or for heavy DM, typically the DM particle will experience more than one collision per crossing. Using the multiscatter-capture formalism (Gould 1992a; Bramante et al. 2017; Dasgupta et al. 2019; Bell et al. 2020; Ilie et al. 2020b), the potential of several classes of objects to constrain properties of DM has been explored in the literature. Below we include a nonexhaustive list of the more recent papers where such effects have been analyzed for Population III stars (Freese et al. 2008; Iocco et al. 2008; Taoso et al. 2008; Ilie & Zhang 2019; Ilie et al. 2020a, 2021), neutron stars (Gould et al. 1990; Bertone & Fairbairn 2008; Kouvaris 2008; Baryakhtar et al. 2017; Bramante et al. 2017; Bell et al. 2018; Chen & Lin 2018; Croon et al. 2018; Raj et al. 2018; Bell et al. 2019; Gresham & Zurek 2019; Hamaguchi et al. 2019; Leung et al. 2019; Joglekar et al. 2020a, 2020b; Acevedo et al. 2020; Bell et al. 2020; Génolini et al. 2020; Keung et al. 2020; Kumar et al. 2020; Leroy et al. 2020; Pérez-García & Silk 2020; Bell et al. 2021a, 2021b; Garani et al. 2021), white dwarfs (Moskalenko & Wai 2007; Bertone & Fairbairn 2008; Miller Bertolami et al. 2014; Bramante et al. 2017; Dasgupta et al. 2019; Horowitz2020; Panotopoulos & Lopes 2020), exoplanets (Leane & Smirnov 2021), and the Earth (Freese 1986; Krauss et al. 1986; Gould 1987, 1992b; Mack et al. 2007).

One common assumption made in most of the aforementioned studies is that the capturing object is at rest with respect to the DM halo. However, as shown by Gould (1987), for the case of single scattering, the capture rates are suppressed when the effects of stellar velocities are included. The aim of this paper is to generalize the result of Gould (1987) and provide an analytic estimation of the suppression coefficient, for the more general case of multiscatter capture of DM.

We end the introduction with a description of the structure of this paper. In Section 2 we consider the case of zero stellar velocity, and briefly review the DM multiscatter-capture formalism (Bramante et al. 2017; Ilie et al. 2020b) and the closed-form analytic approximation formulae of Ilie et al. (2020b, 2021). In Section 3 we present and validate our main result, the velocity-dependent suppression coefficients (see Equations (26)–(29) and Figure 1). These could prove to be useful for future research, as using them, in conjunction with the analytic estimates for the zero-velocity capture rates, can bypass the need for a full numeric, computationally expensive, calculation. Moreover, the suppression coefficients (Equations (26)–(29)) would allow one to quickly rescale, by simply dividing by the corresponding suppression factor, any bounds on DM–nucleon cross section (σ) obtained under the assumption of a capturing object at rest, once the velocity of the capturing star is known. Finally, in Section 4 we revisit the bounds obtained in Ilie et al. (2020a, 2021) by using Population III stars as DM probes. Namely, using the formalism developed in Section 3 we estimate the role of a possible stellar velocity on the projected bounds on σ and find that for Population III stars stellar velocity only weakens the bounds by at most a factor of a few. We end with Section 5, where our conclusions are presented.

Zoom In Zoom Out Reset image size

We want to emphasize that our main results, presented in Section 3, can be applied to any astrophysical object that is capturing DM, and are not restricted to the Sun or Population III stars, which were the focus of this paper. The main reason for us restricting our attention to Population III stars in Section 4, where we estimate projected bounds on σ, is that our main motivation for this work was to reevaluate, as explained above, the forecast bounds previously obtained by our group.

In this section we give a brief overview of the formalism necessary to predict the capture rates of DM by astrophysical compact objects, such as stars, planets, etc. In order for a DM particle to be captured by a star, its velocity must fall below the star's escape velocity. This can occur through collisions with baryonic nuclei in the star. Depending on the mass of the DM particle, this may happen after one (Press & Spergel 1985; Gould 1987, 1988) or more collisions (Gould 1992a; Bramante et al. 2017; Dasgupta et al. 2019; Bell et al. 2020; Ilie et al. 2020b), where very massive DM particles will need more collisions for capture than less massive DM particles. Additionally, the number of collisions that are likely to occur depends on the characteristics of the star and the cross section of interaction between DM and targets inside the star; this number is roughly equal to the optical depth, τ = n_T σ(2R_⋆), where n_T is the average number density of target particles in the star, σ is the cross section of interaction, and R_⋆ is the stellar radius. For all objects considered we will assume one constituent dominates over all others. For example, in the case of Population III stars, or the Sun, which are composed primarily of hydrogen, we assume the atomic nuclei to have the mass of one proton. The probability of capture after exactly N scatters may be represented as follows (Bramante et al. 2017):

$$\begin{eqnarray}&&{C}_{N}=\pi {R}_{\star }^{2}{p}_{N}(\tau ){\int }_{0}^{\infty }f(u)\displaystyle \frac{{du}}{u}{w}^{2}{g}_{N}(w),\end{eqnarray} \tag{ 1 }$$

where R_⋆ is the radius of the star, p_N (τ) is the probability of N collisions occurring, f(u) is the DM velocity distribution, and g_N (w) is the probability that the DM particle's velocity will fall below the escape velocity after N scatters. The quantity p_N (τ) may be represented as (Ilie & Zhang 2019):

$$\begin{eqnarray}&&{p}_{N}(\tau )=\displaystyle \frac{2}{{\tau }^{2}}\left(N+1-\displaystyle \frac{{\rm{\Gamma }}(N+2,\tau )}{N!}\right),\end{eqnarray} \tag{ 2 }$$

with Γ(a, b) is the upper incomplete gamma function, defined as ${\rm{\Gamma }}(a,b)={\int }_{b}^{\infty }{t}^{a-1}{e}^{-t}\,{dt}$. As found in Bramante et al. (2017), g_N (w) can be approximated with

$$\begin{eqnarray}&&{g}_{N}(w)={\rm{\Theta }}\left({v}_{\mathrm{esc}}{\left(1-\displaystyle \frac{{\beta }_{+}}{2}\right)}^{-N/2}-w\right),\end{eqnarray} \tag{ 3 }$$

where ${\beta }_{+}=\tfrac{4{m}_{\chi }{m}_{T}}{{\left({m}_{\chi }+{m}_{T}\right)}^{2}}$, with m_χ being the DM particle mass, m_T the mass of the target particle, and v_esc the escape velocity at the surface of the star. Additionally, w represents the velocity of a DM particle as it enters the star, and is related to its velocity infinitely far away (u) by ${w}^{2}={v}_{\mathrm{esc}}^{2}+{u}^{2}$. This last statement is just conservation of energy. In order to determine the total capture rate, we must sum the values of C_N for every value of N:

$$\begin{eqnarray}&&{C}_{\mathrm{tot}}=\sum _{N=1}^{\infty }{C}_{N}.\end{eqnarray} \tag{ 4 }$$

This is a complete analytical representation of the total capture rate. However, in order to perform a numerical calculation, it is impossible to sum to infinity. Therefore it is necessary to implement a cutoff condition. We continue summing the series up to N_cut, when we reach a desired level of accuracy, which we arbitrarily set to 0.1%; that is, until one additional iteration of C_N only changes C_tot by 0.1%. As shown in Ilie & Zhang (2019), convergence is attained when N_cut ∼ τ, i.e., whenever we sum up to the average number of collisions a DM particle experiences, per crossing, with targets inside the star.

Next we restrict our attention to a capturing object at rest with respect to the DM halo rest frame. In this situation, the DM velocity distribution, f(u), is simply the Maxwell–Boltzmann distribution, f₀ (Gould 1988):

$$\begin{eqnarray}&&{f}_{0}(u){du}={n}_{\chi }\displaystyle \frac{4}{\sqrt{\pi }}{x}^{2}\exp (-{x}^{2}){dx},\end{eqnarray} \tag{ 5 }$$

where n_χ is the number density of DM particles, x is a dimensionless quantity, defined as $x\equiv \sqrt{\tfrac{{m}_{\chi }}{2{T}_{\chi }}}u$, with T_χ representing the DM temperature, which can be related to the thermal average velocity of DM particles: $\bar{v}\equiv \sqrt{3{T}_{\chi }/{m}_{\chi }}$. For this case the integral representation C_N presented in Equation (1) has a closed-form analytic solution (Bramante et al. 2017; Ilie et al. 2020b):

$$\begin{eqnarray}&&{C}_{N}=\displaystyle \frac{1}{3}\pi {R}_{\star }^{2}{p}_{N}(\tau )\displaystyle \frac{\sqrt{6}{n}_{\chi }}{\sqrt{\pi }\bar{v}}\left(\left(2{\bar{v}}^{2}+3{v}_{{esc}}^{2}\right)-\left(2{\bar{v}}^{2}+3{v}_{N}^{2}\right)\exp \left(-\displaystyle \frac{3\left({v}_{N}^{2}-{v}_{{esc}}^{2}\right)}{2{\bar{v}}^{2}}\right)\right),\end{eqnarray} \tag{ 6 }$$

where ${v}_{N}={v}_{\mathrm{esc}}{\left(1-\langle z\rangle {\beta }_{+}\right)}^{-N/2}$, with 〈z〉 the average of the kinematic variable defined that accounts for the scatter angle, and for which a good approximation is 〈z〉 ≈ 1/2 (Bramante et al. 2017).

Several useful analytic approximations for the total capture rates based on summing the C_N 's of Equation (6) have been derived in Ilie et al. (2020b, 2021). We reproduce those results here, for convenience, and future reference:

$$\begin{eqnarray}{C}_{\mathrm{tot}}\approx \left\{\begin{array}{ll}{\left(\displaystyle \frac{2}{3\pi }\right)}^{1/2}\displaystyle \frac{\pi {R}_{\star }^{2}}{{\tau }^{2}}{n}_{\chi }\displaystyle \frac{3{v}_{\mathrm{esc}}^{2}+2{\bar{v}}^{2}}{\bar{v}}{N}_{\mathrm{cut}}({N}_{\mathrm{cut}}+3), & \mathrm{if}\ \,{R}_{v}\ll 1\\ {\left(\displaystyle \frac{6}{\pi }\right)}^{1/2}\displaystyle \frac{\pi {R}_{\star }^{2}}{{\tau }^{2}}{n}_{\chi }\displaystyle \frac{{v}_{\mathrm{esc}}^{4}}{{\bar{v}}^{3}}{\beta }_{+}\langle z\rangle {N}_{\mathrm{cut}}({N}_{\mathrm{cut}}+1)({N}_{\mathrm{cut}}+2)\left(1+\displaystyle \frac{{\beta }_{+}\langle z\rangle }{4}(1+3{N}_{\mathrm{cut}})\right), & \,\mathrm{if}\ \,{R}_{v}\gg 1,\end{array}\right.\end{eqnarray} \tag{ 7 }$$

where we defined ${R}_{v}\equiv \displaystyle \frac{3({v}_{N}^{2}-{v}_{\mathrm{esc}}^{2})}{2{\bar{v}}^{2}}$, and, as pointed out before, the series defining C_tot converges at N_cut ∼ τ.

For most astrophysical objects of interest (with Earth being an important exception), the escape velocity is much larger than the thermal velocity of DM (${v}_{\mathrm{esc}}\gg \bar{v}$). Assuming there is a definite hierarchy between m_χ and m_T , i.e., if m_χ ≫ m_T or m_χ ≪ m_T , Equation (6) could be simplified as

$$\begin{eqnarray}&&\begin{array}{l}{C}_{N}=\sqrt{24\pi }{n}_{\chi }{{GM}}_{\star }{R}_{\star }\displaystyle \frac{1}{\bar{v}}{p}_{N}(\tau )\left(1-\left(1+\displaystyle \frac{2{A}_{N}^{2}{\bar{v}}^{2}}{3{v}_{\mathrm{esc}}^{2}}\right){e}^{-{A}_{N}^{2}}\right),\\ \,\mathrm{where}\ {A}_{N}^{2}={Nk},\end{array}\end{eqnarray} \tag{ 8 }$$

where we defined the following dimensionless parameter:

$$\begin{eqnarray}&&k\equiv 3\displaystyle \frac{\min ({m}_{T};{m}_{\chi })}{\max ({m}_{T};{m}_{\chi })}\displaystyle \frac{{v}_{\mathrm{esc}}^{2}}{{\bar{v}}^{2}}.\end{eqnarray} \tag{ 9 }$$

Using the approximate form of C_N from Equation (8) in Ilie et al. (2020b, 2021) we found closed-form approximations of the total capture rate, which we reproduce below. ³ Those approximations are functionally different, depending on the region of the σ − m_χ parameter space. For the case of multiscattering capture (τ ≳ 1) we find two distinct regimes. First, in the region we called Region I (τ ≳ 1 and k τ ≲ 1):

$$\begin{eqnarray}&&\begin{array}{l}{C}_{\mathrm{tot}}^{I}\approx 5\times {10}^{54}{{\rm{s}}}^{-1}\left(\displaystyle \frac{{\rho }_{\chi }\times \sigma }{\,\mathrm{GeV}\,{\mathrm{cm}}^{-1}}\right){\left(\displaystyle \frac{{10}^{8}\,\,\mathrm{GeV}}{{m}_{\chi }}\right)}^{2}\\ \,\times {\left(\displaystyle \frac{10\,\mathrm{km}\,{{\rm{s}}}^{-1}}{\bar{v}}\right)}^{3}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{3}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-2}.\end{array}\end{eqnarray} \tag{ 10 }$$

In what we called Region II, defined as τ ≳ 1 (multiscatter) and k τ ≳ 1, we find that the capture rates are insensitive to the cross section σ. This essentially means that the cross section is so high ⁴ once we cross the boundary between regions I and II, that the entire DM flux crossing the object gets captured whenever σ is in Region II of the parameter space, and the capture rates saturates:

$$\begin{eqnarray}&&\begin{array}{l}{C}_{\mathrm{tot}}^{{II}}\approx 8\times {10}^{29}\,{{\rm{s}}}^{-1}\left(\displaystyle \frac{{\rho }_{\chi }}{\,\mathrm{GeV}\,{\mathrm{cm}}^{-3}}\right)\left(\displaystyle \frac{{10}^{2}\,\,\mathrm{GeV}}{{m}_{\chi }}\right)\\ \,\times \left(\displaystyle \frac{10\,\mathrm{km}\,{{\rm{s}}}^{-1}}{\bar{v}}\right)\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}.\end{array}\end{eqnarray} \tag{ 11 }$$

Moving on to the single-scattering regime (τ ≲ 1) we find two distinct functional forms of the capture rates, depending on the relative size of the parameter k when compared to unity. In what we called Region III, i.e., τ ≲ 1 and k ≳ 1, we find

$$\begin{eqnarray}&&\begin{array}{l}{C}_{\mathrm{tot}}^{{III}}\approx 4.3\times {10}^{64}\,{{\rm{s}}}^{-1}\left(\displaystyle \frac{{\rho }_{\chi }\times \sigma }{\,\mathrm{GeV}\,{\mathrm{cm}}^{-1}}\right)\left(\displaystyle \frac{{10}^{2}\,\,\mathrm{GeV}}{{m}_{\chi }}\right)\\ \,\times \left(\displaystyle \frac{10\,\mathrm{km}\,{{\rm{s}}}^{-1}}{\bar{v}}\right){\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{2}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 12 }$$

Finally, in Region IV, defined as τ ≲ 1 and k ≲ 1 we find, remarkably, that the capture rate has the exact same parametric form as that of Region I (τ ≳ 1 and k τ ≲ 1):

$$\begin{eqnarray}&&\begin{array}{l}{C}_{\mathrm{tot}}^{{IV}}\approx 5\times {10}^{54}\,{{\rm{s}}}^{-1}\left(\displaystyle \frac{{\rho }_{\chi }\times \sigma }{\,\mathrm{GeV}\,{\mathrm{cm}}^{-1}}\right){\left(\displaystyle \frac{{10}^{8}\,\,\mathrm{GeV}}{{m}_{\chi }}\right)}^{2}\\ \,\times {\left(\displaystyle \frac{10\,\mathrm{km}\,{{\rm{s}}}^{-1}}{\bar{v}}\right)}^{3}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{3}{\left(\displaystyle \frac{{R}_{\star }}{{R}_{\odot }}\right)}^{-2}.\end{array}\end{eqnarray} \tag{ 13 }$$

In summary, in this section we have briefly reviewed the multiscatter DM capture formalism of Bramante et al. (2017). Applying it to the case of zero stellar velocities, we reproduced useful closed-form analytic formulae for the total capture rates, previously obtained in the literature (Ilie et al. 2020b, 2021).

In the next section we move to the main aim of our paper, that of addressing the following question: Is it possible to obtain similar analytic, closed-form formulae for the total capture rates when considering the more general case of a capturing object moving with respect to the DM halo rest frame? As we will show shortly, the answer is yes. This could prove to be useful for future research, as full numeric calculations are computationally expensive, especially when coupling DM capture to stellar evolution codes in order to assess the effects of captured DM annihilations on the stellar structure and evolution.

In this section we present an analytical approximation of the suppression factor for the DM capture rates, in both the low (i.e., k ≫ 1) and high (i.e., k ≪ 1) DM mass regimes. This can be very useful when one needs to estimate the effects of the stellar velocity on DM capture rates, and implicitly on DM scattering cross-section bounds, since calculating numerically the capture rates including the full, boosted Maxwell-Boltzmann distribution can be quite computationally expensive. Our procedure allows one to calculate the simpler, and fully analytically solvable (Ilie et al. 2021), rates when the stellar velocity is neglected, and then apply the suppression factor we derive for any given η. Such a procedure is quite useful when considering capture of DM by astrophysical probes within the solar system neighborhood, where, based on DM profile (Lin & Li 2019) and dispersion velocity (Brown et al. 2009) estimates, we would expect η to be on the order of a few.

In principle, the capture rates of DM by astrophysical objects that have a nonzero velocity with respect to the DM halo rest frame are straightforward to calculate numerically. Essentially, C_tot is still a series obtained by summing the partial capture rates, C_N , as given by Equation (4). The only change now is that, when calculating each C_N numerically via the integral over the DM distribution given in Equation (1), we need to use the appropriate DM distribution. As shown in Gould (1987), this "boosted" distribution (f_η ) can be easily related to the Maxwell–Boltzmann distribution (f₀; see Equation (5)) that would be appropriate to use when the star is stationary:

$$\begin{eqnarray}&&{f}_{\eta }(u)={f}_{0}(u)\exp (-{\eta }^{2})\displaystyle \frac{\sinh \left(2x\eta \right)}{2x\eta },\end{eqnarray} \tag{ 14 }$$

with f₀(u) given in Equation (5). The parameter η represents the dimensionless stellar velocity, $\tilde{v}$, normalized to the dispersion velocity of DM particles in the halo:

$$\begin{eqnarray}&&\eta \equiv \sqrt{\tfrac{3}{2}}\displaystyle \frac{\tilde{v}}{\bar{v}}.\end{eqnarray} \tag{ 15 }$$

Rather than integrating numerically over the DM velocity distribution to calculate a "boosted" capture rate, we can find an equivalent analytical expression. In order to develop this, we followed a similar method to Bramante et al. (2017) and Ilie et al. (2020b). The main difference comes from the use of the boosted velocity distribution as outlined by Gould (1987) instead of the assumption of a Maxwell–Boltzmann distribution. Evaluating Equation (1) by substituting Equation (14) for f(u), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{C}_{N} & = & \displaystyle \frac{{n}_{\chi }\pi {p}_{N}(\tau ){R}^{2}}{2\sqrt{6}\bar{v}\eta }\left[\Space{0ex}{2.50em}{0ex}\displaystyle \frac{1}{\sqrt{\pi }}\exp \left(\displaystyle \frac{-3\left({v}_{N}^{2}-{v}_{\mathrm{esc}}^{2}\right)}{{\bar{v}}^{2}}-2{\eta }^{2}\right)\bar{v}\left\{\Space{0ex}{2.6em}{0ex}4\exp \left(\displaystyle \frac{3({v}_{N}^{2}-{v}_{\mathrm{esc}}^{2})}{{\bar{v}}^{2}}+{\eta }^{2}\right)\bar{v}\eta \right.\right.\\ & & +\ \exp \left(\displaystyle \frac{-6{v}_{\mathrm{esc}}^{2}+6{v}_{N}^{2}-4\sqrt{6}{v}_{\mathrm{esc}}\bar{v}\eta \sqrt{-1+\tfrac{{v}_{N}^{2}}{{v}_{\mathrm{esc}}^{2}}}+4{\bar{v}}^{2}{\eta }^{2}}{4{\bar{v}}^{2}}\right)\left(\sqrt{6}{v}_{\mathrm{esc}}\sqrt{-1+\tfrac{{v}_{N}^{2}}{{v}_{\mathrm{esc}}^{2}}}-2\bar{v}\eta \right)\\ & & -\ \left.\exp {\left(\displaystyle \frac{\sqrt{\tfrac{3}{2}}{v}_{\mathrm{esc}}\sqrt{-1+\tfrac{{v}_{N}^{2}}{{v}_{\mathrm{esc}}^{2}}}}{\bar{v}}+\eta \right)}^{2}\left(\sqrt{6}{v}_{\mathrm{esc}}\sqrt{-1+\tfrac{{v}_{N}^{2}}{{v}_{\mathrm{esc}}^{2}}}+2\bar{v}\eta \right)\right\}\\ & & +\ \left(3{v}_{\mathrm{esc}}^{2}+{\bar{v}}^{2}\left(1+2{\eta }^{2}\right)\right)\mathrm{erf}\left(\displaystyle \frac{\sqrt{\tfrac{3}{2}}{v}_{\mathrm{esc}}\sqrt{-1+\tfrac{{v}_{N}^{2}}{{v}_{\mathrm{esc}}^{2}}}}{\bar{v}}-\eta \right)+2\left(3{v}_{\mathrm{esc}}^{2}+{\bar{v}}^{2}\left(1+2{\eta }^{2}\right)\right)\mathrm{erf}\left(\eta \right)\\ & & -\ \left.\left(3{v}_{\mathrm{esc}}^{2}+{\bar{v}}^{2}\left(1+2{\eta }^{2}\right)\right)\mathrm{erf}\left(\displaystyle \frac{\sqrt{\tfrac{3}{2}}{v}_{\mathrm{esc}}\sqrt{-1+\tfrac{{v}_{N}^{2}}{{v}_{\mathrm{esc}}^{2}}}}{\bar{v}}+\eta \right)\right],\end{array}\end{eqnarray} \tag{ 16 }$$

where the quantity v_N is defined as follows (Bramante et al. 2017):

$$\begin{eqnarray}&&{v}_{N}={v}_{\mathrm{esc}}{\left(1-\displaystyle \frac{{\beta }_{+}}{2}\right)}^{-\tfrac{N}{2}}={v}_{\mathrm{esc}}{\left|\displaystyle \frac{{m}_{\chi }^{2}+{m}_{T}^{2}}{{\left({m}_{\chi }+{m}_{T}\right)}^{2}}\right|}^{-N/2}.\end{eqnarray} \tag{ 17 }$$

As a sanity check, we verify that in the limit of η = 0, the expression reduces to that found in Equation (6):

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\mathrm{lim}}\limits_{\eta \to 0}{C}_{N}=\sqrt{\tfrac{2\pi }{3}}\displaystyle \frac{{n}_{\chi }{p}_{N}(\tau ){R}^{2}}{\bar{v}}\left[(2{\bar{v}}^{2}+3{v}_{\mathrm{esc}}^{2})\right.\\ \,\left.-\exp \left(\displaystyle \frac{-3({v}_{N}^{2}-{v}_{\mathrm{esc}}^{2})}{2{\bar{v}}^{2}}\right)(3{v}_{N}^{2}+2{\bar{v}}^{2})\right],\end{array}\end{eqnarray} \tag{ 18 }$$

which corresponds to the case of zero stellar velocity, i.e., η = 0. The expression of C_N presented in Equation (16) is not particularly illuminating. However, it could be used to evaluate the total capture rate as a series, by adding each C_N from N = 1 until the series has reached the desired level of convergence. This would avoid a full numeric calculation, where each C_N is obtained via Equation (1). Below we present an alternative to this procedure, based on a predictable ratio between the total capture rates when stellar velocities are nonzero and the total capture rate when the stellar velocity is zero, with all other parameters kept the same. A velocity-dependent suppression of the capture rates, for the case of single scattering, was a generic result found in Gould (1987), and simple analytic estimates were provided for two limiting regimes, high and low DM mass (see Equation (2.30) of Gould 1987). In this paper we provide a full analytic, exact, closed-form solution for the suppression factor, valid in the single-scattering regime. Additionally, we generalize this to the case of multiscatter capture of DM. Keeping the same notation as Gould (1987), we define the "suppression factor," labeled as ξ_η , as the ratio between the total capture rate when η is nonzero to the total capture rate for η = 0:

$$\begin{eqnarray}&&{\xi }_{\eta }\equiv \displaystyle \frac{{C}_{\mathrm{tot}}(\eta )}{{C}_{\mathrm{tot}}(\eta =0)}.\end{eqnarray} \tag{ 19 }$$

We next estimate an upper bound and a lower bound on the suppression factor defined above, and identify the conditions under which those two are equal, and as such equal to ξ_η itself. We start by noting that the role of the function g_N (w) in the integrals defining C_N , and correspondingly in the definition of the suppression factor, is to impose a cutoff on the integral. Namely, this amounts to accounting for the DM particles in the tail of the DM distribution that are too fast to be slowed down and captured after N collisions. From Equation (3) one can show that

$$\begin{eqnarray}&&{u}_{\max }(N)={v}_{\mathrm{esc}}{\left({\left(1-\displaystyle \frac{{\beta }_{+}}{2}\right)}^{-N}-1\right)}^{\tfrac{1}{2}}.\end{eqnarray} \tag{ 20 }$$

Defining $a\equiv \min \{{u}_{\max }(N)\}$ and $b\equiv \max \{{u}_{\max }(N)\}$, for future convenience, and by combining Equations (1) and (4), one can show that ξ_η lies between a lower and and upper limit given by

$$\begin{eqnarray}&&\displaystyle \frac{{\int }_{0}^{a}{f}_{\eta }(u)\tfrac{{du}}{u}{w}^{2}}{{\int }_{0}^{b}{f}_{0}(u)\tfrac{{du}}{u}{w}^{2}}\leqslant {\xi }_{\eta }\leqslant \displaystyle \frac{{\int }_{0}^{b}{f}_{\eta }(u)\tfrac{{du}}{u}{w}^{2}}{{\int }_{0}^{a}{f}_{0}(u)\tfrac{{du}}{u}{w}^{2}}.\end{eqnarray} \tag{ 21 }$$

Whenever a = b then the upper bound and the lower bound on the suppression factor (ξ_η ) coincide, and, moreover, they are equal to the suppression factor itself. For the single-scattering case, this happens naturally, since N = 1 so there is only one term in the series $\{{u}_{\max }(N)\}$. Below we provide analytic formulae for the bounds limiting the suppression factor and investigate in detail the conditions under which this can be approximated, not only constrained, in the multiscatter-capture case. Changing variables to the dimensionless ${x}^{2}\equiv \tfrac{3{u}^{2}}{2{\bar{v}}^{2}}$, and introducing the following convenient notations:

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\eta }(t)\equiv \exp \left(-{\eta }^{2}\right){\int }_{0}^{t}\exp \left(-{x}^{2}\right)\left(\displaystyle \frac{3{{xv}}_{\mathrm{esc}}^{2}}{2{v}^{2}}+{x}^{3}\right)\\ \,\times \displaystyle \frac{\sinh (2\eta x)}{2\eta x}\,{dx},\end{array}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{I}_{0}(t)\equiv {\int }_{0}^{t}\exp \left(-{x}^{2}\right)\left(\displaystyle \frac{3{{xv}}_{\mathrm{esc}}^{2}}{2{v}^{2}}+{x}^{3}\right)\,{dx},\end{eqnarray} \tag{ 23 }$$

one can show that the lower and upper bounds of the suppression factor can be recast as

$$\begin{eqnarray}&&{\xi }_{\eta }^{L.B.}=\displaystyle \frac{{I}_{\eta }(x(a))}{{I}_{0}(x(b))},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\xi }_{\eta }^{U.B.}=\displaystyle \frac{{I}_{\eta }(x(b))}{{I}_{0}(x(a))},\end{eqnarray} \tag{ 25 }$$

with $x(a)=\sqrt{\tfrac{3{a}^{2}}{2{\bar{v}}^{2}}}$, and $x(b)=\sqrt{\tfrac{3{b}^{2}}{2{\bar{v}}^{2}}}$, where a and b are the minimum and maximum, respectively, of the sequence $\{{u}_{\max }(N)\}$, defined in Equation (20). Note that I₀(t) is just I_η (t), in the limit of η = 0. For the generic integral I_η (t) we find the following closed-form solution:

$$\begin{eqnarray}&&{I}_{\eta }(t)=\displaystyle \frac{{e}^{-{\eta }^{2}}}{16\eta }\left(4\eta +2\sqrt{\pi }{e}^{{\eta }^{2}}\left(2{\eta }^{2}+1\right)\mathrm{erf}(\eta )+A(\eta ;t)+B(\eta ;t)\right),\end{eqnarray} \tag{ 26 }$$

with

$$\begin{eqnarray}&&A(\eta ;t)\equiv {e}^{-t(2\eta +t)}\left[\sqrt{\pi }\left(2{\eta }^{2}+1\right){e}^{{\left(\eta +t\right)}^{2}}\left(\mathrm{erf}(t-\eta )-\mathrm{erf}(\eta +t)\right)-2\eta -2{e}^{4\eta t}(\eta +t)+2t\right],\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&B(\eta ;t)\equiv \displaystyle \frac{3\sqrt{\pi }{e}^{{\eta }^{2}}{v}_{\mathrm{esc}}^{2}(2\mathrm{erf}(\eta )+\mathrm{erf}(t-\eta )-\mathrm{erf}(\eta +t))}{{\bar{v}}^{2}}.\end{eqnarray} \tag{ 28 }$$

Equations (24)–(28) can be used to compute exactly the lower and upper bounds for the suppression factor. Moreover, for single-scattering capture (i.e., τ ≪ 1) the same set of equations predict exactly the value of the suppression factor, since $x(a)=x(b)=\sqrt{\tfrac{3{u}_{\max }{\left(1\right)}^{2}}{2{\bar{v}}^{2}}}$. To gain further insight it is instructive to take limiting cases. First, we explicitly write the analytic form we find for I₀(t):

$$\begin{eqnarray}&&{I}_{0}(t)=\displaystyle \frac{1}{2}\left(\displaystyle \frac{3{e}^{-{t}^{2}}\left({e}^{{t}^{2}}-1\right){v}_{\mathrm{esc}}^{2}}{2{\bar{v}}^{2}}-{e}^{-{t}^{2}}\left({t}^{2}+1\right)+1\right).\end{eqnarray} \tag{ 29 }$$

We will next restrict our attention to the case when there is a definite hierarchy between m_χ and m_T , i.e., when one of those mass scales is larger than the other. In this case the parameter β₊ is much less than unity, and can be approximated as

$$\begin{eqnarray}&&{\beta }_{+}\approx 4\displaystyle \frac{\min ({m}_{T};{m}_{\chi })}{\max ({m}_{T};{m}_{\chi })}.\end{eqnarray} \tag{ 30 }$$

We can now approximate the terms in the sequence $\{{u}_{\max }(N)\}$, defined by Equation (20):

$$\begin{eqnarray}&&{u}_{\max }(N)\approx {v}_{\mathrm{esc}}{\left(N{\beta }_{+}/2\right)}^{1/2}.\end{eqnarray} \tag{ 31 }$$

The last two equations combined with the definition of k from Equation (9) can be used to show that $x(a)\approx \sqrt{k}$ and $x(b)\approx \sqrt{k\max (\tau ;1)}$. In the last step we used the fact that for multiscatter capture ${N}_{\max }\approx \tau$. In order to keep the treatment of single scatter and multiscatter unified we used $\max (\tau ;1)$, since ${N}_{\max }=N=1$ for single scattering. However, as discussed before, for single scatter the upper and lower bounds of Equations (24)–(25) coincide, since in that case $x(a)=x(b)=\sqrt{k}$. We are now in a position to derive the limiting cases of the lower and upper bounds of the suppression factor of Equations (24)–(25). We start with the case of single-scattering capture, for which there are two natural regimes: the k ≫ 1 regime and the k ≪ 1 regime. From the definition of k in Equation (9) one can find that k ≫ 1 is valid whenever ${\left(3{v}_{\mathrm{esc}}^{2}/{\bar{v}}^{2}\right)}^{-1}{m}_{T}\ll {m}_{\chi }\ll (3{v}_{\mathrm{esc}}^{2}/{\bar{v}}^{2}){m}_{T}$, whereas k ≪ 1 otherwise. In the k ≫ 1 limit (corresponding to low m_χ ) we perform an asymptotic expansion of I_η and I₀ from Equations (26)–(29) around $x(a)=x(b)=\sqrt{k}\to \infty$. Keeping only leading-order terms, we get

$$\begin{eqnarray}&&{\xi }_{\eta }^{k\gg 1}\approx \displaystyle \frac{{e}^{-{\eta }^{2}}\left(4\eta +\tfrac{6{e}^{{\eta }^{2}}\sqrt{\pi }{v}_{\mathrm{esc}}^{2}\mathrm{erf}\left(\eta \right)}{{\bar{v}}^{2}}+2{e}^{{\eta }^{2}}\sqrt{\pi }\left(1+2{\eta }^{2}\right)\mathrm{erf}\left(\eta \right)\right)}{4\left(2+\tfrac{3{v}_{\mathrm{esc}}^{2}}{{\bar{v}}^{2}}\right)\eta }.\end{eqnarray} \tag{ 32 }$$

Whenever ${v}_{\mathrm{esc}}^{2}\gg {\bar{v}}^{2}$, and η is not much larger than unity we can further simplify the previous result to ${\xi }_{\eta }^{k\gg 1}\approx \tfrac{\sqrt{\pi }\mathrm{erf}(\eta )}{2\eta }$, which matches the result of Gould (see Equation (2.30) of Gould 1987).

While still in the single-scatter regime, at either very high or very low m_χ , the parameter k becomes much less than unity. Therefore we simply Taylor-expand I_η and I₀ from Equations (26)–(23) around $x(a)=x(b)=\sqrt{k}\approx 0$. Neglecting terms of ${ \mathcal O }({k}^{3})$, we find that the suppression factor, ξ_η , can be approximated with

$$\begin{eqnarray}&&\begin{array}{l}{\xi }_{\eta }^{k\ll 1}\approx {e}^{-{\eta }^{2}}+\displaystyle \frac{1}{3}{e}^{-{\eta }^{2}}{\eta }^{2}k\\ \,+\displaystyle \frac{{e}^{-{\eta }^{2}}{\eta }^{2}{k}^{2}\left(10{\bar{v}}^{2}+3\left(4{\eta }^{2}-5\right){v}_{\mathrm{esc}}^{2}\right)}{270{v}_{\mathrm{esc}}^{2}}.\end{array}\end{eqnarray} \tag{ 33 }$$

We emphasize once more that, for the single-scattering regime, the full, nonapproximated, functional form of the suppression factor that can be obtained from ξ_η = I_η (x(a))/I₀(x(a)), with $x(a)=\sqrt{3{u}_{\max }^{2}(1)/2{\bar{v}}^{2}}$, and I_η given in Equation (26) and I₀ from Equation (29). However, the approximations derived above allow one to gain some additional insight. In the low-m_χ regime, defined by the k ≫ 1 condition, the suppression factor has a roughly constant value, given by Equation (32). Once m_χ becomes either extremely low or extremely large, such that k crosses unity, and now becomes less than one, the suppression factor starts to change significantly, in an approximately polynomial fashion, according to Equation (33). Whenever k becomes much less than unity, the suppression factor asymptotes to ${e}^{-{\eta }^{2}}$, matching the result found by Gould (see Equation (2.30) of Gould 1987). We note here that k ≪ 1 is equivalent to $\bar{v}\gg \sqrt{3\tfrac{\min ({m}_{T};{m}_{\chi })}{\max ({m}_{T};{m}_{\chi })}}{v}_{\mathrm{esc}}$. Therefore, whenever the DM dispersion velocity is high compared to the escape velocity, only a small fraction of the DM particles will be captured, as most of them will have speeds larger than the escape velocity.

We next move our focus to the multiscatter-capture case. For all objects we considered, it turns out that k τ ≪ 1, given the present bounds on σ from direct-detection experiments. In turn, that means that for the multiscatter case $x(b)\approx \sqrt{k\tau }\ll 1$. Moreover, the same bounds on σ imply that $x(a)\approx \sqrt{k}\ll 1$, whenever τ ≫ 1, i.e., in the multiscatter regime. Expanding around $x(a)\approx \sqrt{k}\approx 0$ and $x(b)\approx \sqrt{k\tau }\approx 0$ we get the following approximations for the lower and upper bounds of ξ_η , defined in Equations (24)–(25):

$$\begin{eqnarray}&&{\xi }_{\eta }^{L.B.}\approx \displaystyle \frac{{e}^{-{\eta }^{2}}}{\tau },\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{\xi }_{\eta }^{U.B.}\approx \min ({e}^{-{\eta }^{2}}\tau ;1).\end{eqnarray} \tag{ 35 }$$

While those bounds can be useful, it turns out that in most cases of interest the suppression factor itself can be well approximated with ${\xi }_{\eta }\approx {e}^{-{\eta }^{2}}$, as shown below. Each of the C_N in the definition of the total capture rate (C_tot) is defined as per Equation (1). Under the conditions we explore here, i.e., when there is a mass hierarchy between m_χ and m_T , and using the integrals I_η defined in Equations (26)–(23), one can show that

$$\begin{eqnarray}&&{C}_{N}\sim {p}_{N}(\tau ){I}_{\eta }({kN}),\end{eqnarray} \tag{ 36 }$$

up to constants independent of N. For the case of η = 0 we have C_N (η = 0) ∼ p_N (τ)I₀(kN). As explained above, bounds on σ from direct-detection experiments imply that, in most cases of interest, k ≪ 1 and kN_cut ≪ 1, once τ ≫ 1, i.e., for multiscatter capture. Therefore we can expand both I_η and I₀ around zero. Keeping leading-order terms we have, up to constants independent of N, the following scaling relations: ${C}_{N}\sim {e}^{-{\eta }^{2}}{Nk}$ and C_N (η = 0) ∼ Nk. It is important to note that in both terms the same constants were "ignored." As such, the suppression factor in the multiscatter regime becomes simply ${\xi }_{\eta }^{M.S.}\approx {e}^{-{\eta }^{2}}$. Note that this is a smooth continuation of the asymptotic behavior found in the single-scattering regime (see discussion in the paragraph following Equation (33)).

In Figure 1 we validate our analytic results for the suppression factor, against the full numeric result, using the Sun as a sample capturing object. In order to explicitly show the dependence of the suppression factor on the stellar velocity, in Figure 2 we consider a Sun-like star for which we arbitrarily set η = 5, with all other parameters being fixed. Note the significant suppression in this case, when contrasted to an object in the solar system neighborhood, where η ≈ 1. We point out that for the Milky Way, as demonstrated by observed rotation curves, the value of the stellar velocities, and in turn the value of η, is roughly constant, for stars farther than a few kiloparsecs from the Galactic center. As such, our choice of η = 5 should be viewed as a hypothetical example, only for the purpose of illustrating how rapidly the exponential suppression factor can reduce capture rates, even for order-unity values of η.

Zoom In Zoom Out Reset image size

In summary, in this section we have derived and validated simple analytical formulae for the suppression factors of the capture rates in terms of the dimensionless stellar velocity $\eta \equiv \sqrt{\tfrac{3}{2}}\tfrac{\tilde{v}}{\bar{v}}$ for both single and multiscatter capture of DM. In the next section we explore the effects of the suppression of the capture rates by stellar velocities in the context of Population III stars as DM probes.

In the previous section, we found an analytic closed form for the suppression in capture rate due to the relative velocity between a star and the DM halo. In this section we address the following question: If a Population III star does not form precisely at the center of the DM halo and, therefore, has some orbital velocity, how will this affect the constraining power of Population III stars on DM parameters such as the DM–nucleon scattering cross section?

In most cases, a Maxwell–Boltzmann velocity distribution has been applied to calculate DM capture rates of Population III stars. This is because it is typically assumed that they form at the center of DM halos. As a result, they do not have a velocity relative to the halo. Simulations demonstrate that Population III stars would form near the center of DM halos in low multiplicity (Barkana & Loeb 2001; Abel et al. 2002; Bromm & Larson 2004; Yoshida et al. 2006, 2008; Loeb 2010; Bromm 2013; Machida & Doi 2013; Klessen 2010). These stars would orbit around the center of the DM halo. They therefore have a relative velocity directly related to the star's distance from the halo's center. We point out that the formalism developed here, and the analytical approximations, are valid for any DM-capturing object, such as stars, neutron stars, and brown dwarfs, and we use Population III stars just as an example of how to apply it.

In order to isolate the effects of the stellar velocity on the capture rate, we first consider the extreme case, where the star forms at the scale radius of the halo. This scenario is highly unlikely, as Population III stars form much closer to the center of the DM halo; as such, the suppression due to Population III stars' stellar velocities is expected to be always less than whatever suppression we will find for this benchmark, overly conservative case. At first pass we assume that the halo follows a Navarro–Frenk–White (NFW) profile (Navarro et al. 1997):

$$\begin{eqnarray}&&{\rho }_{\mathrm{halo}}=\displaystyle \frac{{\rho }_{0}}{\tfrac{r}{{r}_{s}}{\left(1+\tfrac{r}{{r}_{s}}\right)}^{2}},\end{eqnarray} \tag{ 37 }$$

where r is the distance from the center and r_s is the scale radius, and for DM mini-halos in which Population III stars form, it has a value that ranges between 3 and 300 parsecs. ρ₀ is the central density, defined as

$$\begin{eqnarray}&&{\rho }_{0}=\displaystyle \frac{200}{3}\displaystyle \frac{{c}_{\mathrm{vir}}^{3}}{\mathrm{ln}(1+{c}_{\mathrm{vir}})-\tfrac{{c}_{\mathrm{vir}}}{{c}_{\mathrm{vir}}+1}}{\rho }_{c},\end{eqnarray} \tag{ 38 }$$

where c_vir represents a concentration parameter ${c}_{\mathrm{vir}}=\tfrac{{r}_{\mathrm{vir}}}{{r}_{s}}$ and ranges in value from 1 to 10 (Freese et al. 2009). ρ_c is the critical density and depends on the redshift, z, in accordance with the Friedmann equation.

Knowing the density distribution of DM in the halo, we calculate the mass enclosed in the stellar orbit, and thus can easily find the speed at which a Population III star located at this point would orbit around its center. The stellar velocity will be encoded in a parameter called η, a dimensionless quantity which is defined as in Equation (15), which we reproduce here for convenience:

$$\begin{eqnarray*}&&\eta \equiv \sqrt{\tfrac{3}{2}}\displaystyle \frac{\tilde{v}}{\bar{v}},\end{eqnarray*}$$

with $\tilde{v}$ representing the stellar velocity and $\bar{v}$ the dispersion velocity of DM. Adopting the parameters described in the caption to Figure 3, we expect an object located at the scale radius of the halo to have an orbital velocity of $\tilde{v}\approx 5.22\times {10}^{5}\,\mathrm{cm}\,{{\rm{s}}}^{-1}$ when placed in a standard NFW profile. Of course, changing the redshift or concentration parameter, for instance, would yield slightly different values; we provide an analysis adopting z = 15 and c = 10 in order to illustrate one example in depth. Including the effects of the adiabatic compression (Young 1980; Blumenthal et al. 1986; Freese et al. 2009; Gnedin et al. 2011) on the DM-density profile would not much affect this value, since the mass enclosed within the scale radius will stay roughly constant, as the adiabatic compression operates at smaller, subparsec scales. The relation between the value of η and the distance from the halo center is shown in Figure 3 for a standard NFW profile as well as two adiabatically contracted (AC) profiles. As the baryonic molecular cloud collapses to form a protostar, the DM orbits respond to this enhancement of the gravitational potential by becoming more tightly packed, a consequence of conservation of adiabatic invariants, such as angular momentum or radial action. This is, in essence, what in the literature is called "adiabatic contraction." As commonly done in the literature (see Freese et al. 2009, for example) we use the standard Blumenthal et al. (1986) formalism to estimate the DM densities. This formalism assumes circular DM orbits, and, as such, the only relevant adiabatic invariant being angular momentum.

Zoom In Zoom Out Reset image size

We elaborate below some of the details of the calculation of the dimensionless stellar velocity η. The mass profile of the halo is found by integrating over the density profile considered:

$$\begin{eqnarray}&&M(r)={\int }_{0}^{{r}_{\mathrm{vir}}}{\rho }_{\mathrm{halo}}(r)\times 4\pi {r}^{2}{dr},\end{eqnarray} \tag{ 39 }$$

where, for an NFW profile, we obtain

$$\begin{eqnarray}&&M(r)=4\pi {r}_{s}^{3}{\rho }_{0}\left[-\displaystyle \frac{r}{(1+\tfrac{r}{{r}_{s}}){r}_{s}}+\mathrm{ln}\left(1+\displaystyle \frac{r}{{r}_{s}}\right)\right].\end{eqnarray} \tag{ 40 }$$

After substituting the mass profile into η, we obtain the following expression for a standard NFW profile:

$$\begin{eqnarray}&&\eta (r)=\displaystyle \frac{\sqrt{\pi }}{{\rm{50,000}}}\sqrt{\tfrac{{c}_{\mathrm{vir}}^{3}G{\rho }_{c}{{\rm{r}}}_{s}^{3}\left(\mathrm{ln}\left(\tfrac{{\rm{r}}}{{{\rm{r}}}_{s}}+1\right)-\tfrac{{\rm{r}}}{{{\rm{r}}}_{s}\left(\tfrac{{\rm{r}}}{{{\rm{r}}}_{s}}+1\right)}\right)}{r\left(\mathrm{ln}({c}_{\mathrm{vir}}+1)-\tfrac{{c}_{\mathrm{vir}}}{{c}_{\mathrm{vir}}+1}\right)}}.\end{eqnarray} \tag{ 41 }$$

Knowing the velocity, and the value of η, for a Population III star at a given distance, we can now apply Equation (14). We choose to select the scale radius of the DM halo as a reasonable maximum bound to use when considering boosted capture of Population III stars, because, in practice, these stars are expected to form well inside the scale radius of DM halos. In turn, this will lead to the highest possible suppression on the previously calculated capture rates in Ilie & Zhang (2019) and Ilie et al. (2020a, 2021). Note that in our capture-rate calculation, since we are mainly focusing on the effect of stellar velocity, we take the assumption that the DM density is fixed, i.e., is the same value at the scale radius as at the halo center.

In order to numerically calculate the capture rate of DM, we need to adopt parameters of Population III stars from numerical simulations. Although Population III stars are still theoretical objects and have not been observed, simulations have been done, such as for example in Iocco et al. (2008) and Ohkubo et al. (2009). In Ilie & Zhang (2019), it has been shown that Population III stars have two different homology scaling relations (in two different mass regimes), where stars with a mass M_⋆ < 20 M_⊙ follow ${R}_{\star }\propto {M}_{\star }^{0.21}$, and larger mass stars follow ${R}_{\star }\propto {M}_{\star }^{0.56}$.

Since our aim in this paper is to understand and quantify the effects of the stellar velocity on DM capture, we assume, for now, the same ambient density at the location of the star, in order to disambiguate between the suppression due to an increase in the stellar velocity, and the decrease in the DM density. Both of those lead to a suppression in the capture rate. For the latter, the effect is trivial, since the total capture rate scales linearly with the DM density: C_tot ∼ ρ_χ . Our aim is to obtain a simple, analytic procedure that would estimate the suppression rate on capture rates by any astrophysical object, if the parameter η is known. Previous works in the literature that use compact astrophysical objects as DM probes typically neglect the effects of the stellar velocity. For example, neutron stars were considered by Bramante et al. (2017), and exoplanets by Leane & Smirnov (2021), and both works neglect the possible role of the relative velocity between the capturing object and the DM halo. The formalism we will develop in Section 3 can be easily applied to any such scenario, if the location (and therefore velocity) of the object in question is known.

In Figure 4 we contrast the total capture rates of DM by an arbitrary Population III star, first placed at the center of the DM halo (as previously assumed) and then placed at the scale radius of the DM halo. We note that, to a good approximation, the capture rates remain unaffected by the inclusion of the stellar velocity, for the case of Population III stars. When the boosted distribution is applied, the DM capture rate (see Equation (1)) is suppressed, as one may expect. We next take the ratio between the capture rates calculated using a boosted (η ≠ 0) and a regular (η = 0) Maxwell–Boltzmann distribution to illustrate the amount by which capture is suppressed. As shown in Figure 5, the ratio plateaus for low and high DM masses. Notice that the drastic change in this ratio occurs when the DM mass reaches 10⁵ GeV, which corresponds exactly to the m_χ for which the quantity $k=3\tfrac{{m}_{T}}{{m}_{\chi }}\tfrac{{v}_{\mathrm{esc}}^{2}}{{\bar{v}}^{2}}$, defined in Equation (9), reaches a value of 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Prior works such as Ilie et al. (2020a, 2021) constrain the bounds on the cross section of interaction between DM and baryonic particles due to the impact DM has on the luminosity of Population III stars. Any object that is gravitationally bound, such as a star, will have an upper bound on how bright it can shine, at a given mass, i.e., the Eddington limit:

$$\begin{eqnarray}&&{L}_{\mathrm{Edd}}\leqslant {L}_{\mathrm{nuc}}+{L}_{\mathrm{DM}},\end{eqnarray} \tag{ 42 }$$

where L_nuc is the luminosity due to nuclear fusion, and L_DM is the additional luminosity provided by DM annihilations, which is directly related to the amount of DM captured:

$$\begin{eqnarray}&&{L}_{\mathrm{DM}}=f\,{C}_{\mathrm{tot}}\,{m}_{\chi },\end{eqnarray} \tag{ 43 }$$

where f is the efficiency with which DM annihilation contributes to the luminosity of the star, i.e., the amount of energy thermalized with the star. The remainder 1 − f is lost to products of annihilation that escape, such as neutrinos. Because C_tot is dependent on σ, we can numerically calculate the maximum expected value of the cross section by finding the maximum value of L_DM. Recall that the Eddington luminosity is given by

$$\begin{eqnarray}&&{L}_{\mathrm{Edd}}=\displaystyle \frac{4\pi \,c\,G\,{M}_{\star }}{{\kappa }_{\rho }},\end{eqnarray} \tag{ 44 }$$

where c is the speed of light, G is the gravitational constant, M_⋆ is the stellar mass, and κ_ρ is the opacity of the stellar atmosphere.

The value of L_nuc is dependent on the mass of the star (x ≡ M_⋆/M_⊙), and here we use the fitting form found in Ilie et al. (2021), which we reproduce here for convenience:

$$\begin{eqnarray}&&{L}_{{nuc}}\simeq {10}^{\tfrac{\mathrm{log}\left(3.71\times {10}^{4}{L}_{\odot }{\rm{s}}\,{\mathrm{erg}}^{-1}\right)}{1+\exp (-0.85\mathrm{log}(x)-1.95)}}\cdot {x}^{\tfrac{2.01}{{x}^{0.48}+1}}\ \mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 45 }$$

We note here that the above equation does not take into account the effect DM heating has on the internal structure of the star, specifically on the core temperature that directly affects the nuclear luminosity. Moreover, we ignored the DM heating effects on the stellar radius, which in turn affects the capture rate. We have used these models in order to facilitate comparison with earlier work. An accurate calculation requires incorporating DM heating into a stellar structure code, which is beyond the scope of this paper. ⁵

Since the capture rate is suppressed when including the effects of the stellar velocity, the bounds shift upwards by exactly a factor of ${\xi }_{\eta }^{-1}$ and become less stringent, as illustrated in Figure 6. However, note that the values of η = 0.69 and η = 0.75 considered here are, for Population III stars, unrealistically high. That is because they correspond to the star at the scale radius of the DM halo, which is many orders of magnitude above the typically expected maximum tens of au from the center where Population III stars form. Even with these exaggerated values of η we note that the suppression in the capture rates, and correspondingly the weakening of the cross-section bounds are, at most, approximately 55%. This suppression has a negligible effect on the constraints placed on the DM–nucleon cross section, as demonstrated in Figure 6. There is no significant difference in the bounds on the cross section for all of the values of η tested. Note that the bounds on σ shown here are calculated both numerically assuming a boosted distribution throughout the calculation, and by rescaling bounds found under the assumption η = 0 by a factor of the inverse of the suppression factor as in ξ_η . We point out that both of these methods produce an exact match (in the figure, the green and pink overlap exactly).

Zoom In Zoom Out Reset image size

In this paper we derived and validated an analytic closed form of the suppression factor for the capture rates of DM by astrophysical objects that have a nonzero velocity with respect to the DM halo: Equations (26) through (29). One of the most useful applications of these formulae is that they allow the immediate rescaling of any bounds previously obtained, for any object, under the assumption of zero stellar velocity. Namely, if the stellar velocity is determined, all one needs to do is to rescale the previously obtained bounds σ(η = 0) with the inverse of the suppression factor ξ_η . For the case of Population III stars as DM probes, we find that the role of the stellar velocity can be safely neglected, and all of our previous results, where Population III were considered to be at rest with respect to the DM halo, remain largely unchanged. This is because the DM capture rate is suppressed by a factor of 57% at the most. This happens for high-mass DM particles (m_χ ≳ 10⁷ GeV), and when the star is considered to have formed—or migrated—all the way to the scale radius of the DM halo, which is a highly unrealistic scenario. In most cases Population III stars will live much closer to the center of the DM halo, within the inner 10 au or so, leading to much higher suppression rates. Of course, the instance of DM halo mergers would change these results, as the location of stars could change significantly in the process. Our formalism is even more relevant for astrophysical objects within the Milky Way that act as DM probes, such as neutron stars, brown dwarfs, and exoplanets. In this case, the most promising location in terms of high DM density, the center of the Milky Way, is the site of a supermassive black hole, which would lead to large orbital velocities, when compared to the centers of high-redshift DM microhalos, and therefore larger suppression factors for the DM capture. Moreover, for very dim probes, such as neutron stars, the most optimal location would be in the vicinity of the solar system, where the suppression factor would be even more significant, and therefore important to take into account and estimate.

J.P. thanks the financial support from Colgate University, via the Research Council student wage grant, and the Justus '43 and Jayne Schlichting Student Research Funds.