Multiscatter capture of superheavy dark matter by Pop. III stars

TL;DR

This paper investigates how superheavy dark matter particles can be captured by Population III stars through multiscatter interactions, potentially affecting star evolution and providing a novel way to constrain dark matter properties.

Contribution

It introduces an analytical formalism for multiscatter capture of superheavy dark matter by Pop. III stars and explores its implications for star evolution and dark matter constraints.

Findings

Upper bounds on dark matter capture rates by Pop. III stars.

Potential influence of dark matter annihilation heat on star evolution.

Constraints on Pop. III star masses based on dark matter interactions.

Abstract

If captured by the gravitational field of stars or other compact objects, dark matter can self-annihilate and produce a potentially detectable particle flux. In the case of superheavy dark matter ($ m_{X} \gtrsim 10^{8} GeV $), a large number of scattering events with nuclei inside stars are necessary to slow down the dark matter particles below the escape velocity of the stars, at which point the Dark Matter (DM) particle becomes trapped, or captured. Using the recently developed analytical formalism for multiscatter capture, combined with the latest results on the constraints of dark-matter-baryon scattering cross-section, we calculate upper bounds on the capture rates for superheavy dark matter particles by the first (Pop. III) stars. Assuming that a non-zero fraction of the products of captured superheavy dark matter (SHDM) annihilations can be trapped and thermalized inside the star we find that this additional heat source could influence the evolutionary phase of Pop. III stars. Moreover, requiring that Pop. III stars shine with sub-Eddington luminosity, we find upper bounds on the masses of the Pop. III stars. This implies a DM dependent cutoff on the initial mass function (IMF) of Pop. III stars, thus opening up the intriguing possibility of constraining DM properties using the IMF of extremely metal-poor stars.

Tables

Table 1. Table 1 : Stellar mass ( M ⋆ subscript 𝑀 ⋆ M_{\star} ), radius ( R ⋆ subscript 𝑅 ⋆ R_{\star} ), surface escape velocity ( V e ​ s ​ c subscript 𝑉 𝑒 𝑠 𝑐 V_{esc} ), and luminosity ( L ⋆ subscript 𝐿 ⋆ L_{\star} ) in solar units, for the Pop. III models of  [ 10 ] and  [ 50 ] we consider in this paper.

$M_{\star }\left(M_{\odot }\right)$	$R_{\star }\left(R_{\odot }\right)$	$V_{esc}\left(V_{\odot }\right)$	$L_{\star }\left(L_{\odot }\right)$
5	1.203	2.06	$8.41\times {10}^2$
7	1.286	2.36	$2.39\times {10}^3$
9	1.347	2.61	$5.07\times {10}^3$
12	1.43	2.93	$1.18\times {10}^4$
15	1.48	3.22	$2.19\times {10}^4$
20	1.65	3.52	$5.00\times {10}^4$
40	2.57	3.98	$2.86\times {10}^5$
100	4.25	4.90	$1.45\times {10}^6$
200	6.14	5.76	$3.97\times {10}^6$
400	9.03	6.72	$9.89\times {10}^6$
600	11.24	7.38	$1.61\times {10}^7$
1000	14.83	8.29	$3.50\times {10}^7$

Equations

$$N\approx n_T{\sigma }_{n}2R_{\star }$$

$$0\le \frac{\Delta E_i}{E_i}\le {\beta }_{+}.$$

$$\frac{dC}{dV{d}^{3}u}=dF(n_X,u,v_{star},v_{esc}^{halo})\Omega (n_T(r),\omega (r),{\sigma }_n,m_n,m_X),$$

$$C_{tot}(m_X)=\sum _{N=1}^{\infty }{C}_N,$$

$$p_N(\tau )=\frac{2}{{\tau }^2}\left(N+1-\frac{\Gamma (N+1,\tau )}{N!}\right).$$

$$E_N=\prod _{i=1}^N(1-z_i{\beta }_{+})E_0,v_N=\prod _{i=1}^N(1-z_i{\beta }_{+}{)}^{1/2}w.$$

$$g_N(w)={\int }_0^{1}d{z}_1{\int }_0^{1}d{z}_2\cdot \cdot \cdot {\int }_0^{1}d{z}_N\Theta \left(v_{esc}\prod _{i=1}^N(1-z_i{\beta }_{+}{)}^{-1/2}-w\right).$$

$$g_N(w)=\Theta \left(v_{esc}\prod _{i=1}^N(1-z_i{\beta }_{+}{)}^{-1/2}-w\right).$$

$$C_N=\pi R^2{p}_N(\tau ){\int }_{v_{esc}}^{\infty }f(u)\frac{dw}{u^2}{w}^3{g}_N(w),$$

$$C_N$$

$$A_N^2$$

$${\rho }_{halo}=\frac{{\rho }_0}{\frac{r}{r_s}{\left(1+\frac{r}{r_s}\right)}^2},$$

$${\rho }_0={\rho }_{crit}(Z)\frac{200}{3}\frac{c^3}{\ln (1+c)-c/(c+1)},$$

$$⟨{\overline{v}}^2⟩=\frac{\overline{W}}{M_{halo}},$$

$$W=-4\pi G\int {\rho }_{halo}{M}_{halo}(r)rdr$$

$$\overline{v}=10km/s.$$

$$\begin{gathered} C_N\simeq 1.35\times 10^{43}{s}^{-1}\left(N+1-\frac{\Gamma (N+2,\tau )}{N!}\right)\left[1-\left(1-2\times 10^{-8}N\left(\frac{10^{8}GeV}{m_X}\right)\right)e^{-A_N^2}\right] \\ \left(\frac{10^{8}GeV}{m_X}\right){\left(\frac{1.26\times 10^{-40}c{m}^2}{{\sigma }_n}\right)}^2\left(\frac{10km/s}{\overline{v}}\right)\left(\frac{{\rho }_X}{10^{9}GeV/c{m}^3}\right)\left(\frac{M_{\odot }}{M_{\star }}\right){\left(\frac{R_{\star }}{R_{\odot }}\right)}^5, \end{gathered}$$

$$A_N^2$$

$$\tau$$

$$N+1-\frac{\Gamma (N+2,\tau )}{N!}\approx \frac{{\tau }^{N+2}}{N!(N+2)}+\mathcal{O}({\tau }^{N+3}).$$

$${\sigma }_n(m_X)\lesssim 1.26\times 10^{-40}\left(\frac{m_X}{10^8\mathrm{GeV}}\right).$$

$$|(C_{tot,N_{cutoff}+1}/C_{tot,N_{cutoff}})-1|\le \mathrm{0.001.}$$

$$\left[1-\left(1-\frac{2A_N^2{\overline{v}}^2}{3v_{esc}^2}\right)e^{-A_N^2}\right]\approx A_N^2\approx 1.10\times 10^{-4}N\left(\frac{10^8\mathrm{GeV}}{m_X}\right){\left(\frac{10km/s}{\overline{v}}\right)}^2\left(\frac{M_{\star }}{M_{\odot }}\right)\left(\frac{R_{\odot }}{R_{\star }}\right).$$

$$C_N$$

$$\frac{R_{\star }}{R_{\odot }}$$

$$N_{cutoff}\propto N_{average}\propto n_T{R}_{\star }\propto \frac{M_{\star }}{R_{\star }^2}\propto M_{\star }^{1-2p}.$$

$$C_{tot}$$

$$\sum _{N=1}^{N=\infty }{p}_N(\tau )N=N_{avg}\simeq \tau$$

$$C_{tot}\simeq 10^{24}{s}^{-1}\left(\frac{{\sigma }_n}{1.26\times 10^{-40}c{m}^2}\right)\left(\frac{{\rho }_X}{10^{9}GeV/c{m}^3}\right){\left(\frac{10^{8}GeV}{m_X}\right)}^2{\left(\frac{10km/s}{\overline{v}}\right)}^3{\left(\frac{M_{\star }}{M_{\odot }}\right)}^3{\left(\frac{R_{\odot }}{R_{\star }}\right)}^2$$

$$\dot{N}=C-2{\Gamma }_A,$$

$${\Gamma }_A=\frac{1}{2}C.$$

Full text

1 Introduction

The first stars formed via the gravitational collapse of zero metallicity primordial baryonic gas clouds that contain H and He from big bang nucleosynthesis. The gas collapses on the gravity well provided by dark matter, at the center of mini-halos ($M_{halo}\sim 10^{6}M_{\odot}$). Those DM halos form in the early universe via the gravitational growth of the initial density perturbations provided by cosmic inflation, and subsequent mergers, i.e. hierarchical structure formation. The emergence of the first stars marks the end of the “dark ages” of the universe. Using hydrodynamical numerical simulations of gas cloud collapse in the early universe, the following picture emerges: the first generation of stars (Pop. III stars) have formed at redshifts $z=10-50$, when the universe was roughly 200 million years old. At that stage, the main cooling mechanism that allows the baryonic gas cloud to collapse towards the center of the halo and form a protostar is molecular hydrogen cooling. Collapse itself leads to heating of the protostellar gas cloud, and when a balance between heating and cooling is achieved, a protostar is formed (i.e. an object in hydrostatic equilibrium supported against further collapse by radiation and thermal gas pressure). This object continues to accrete material and become hotter and hotter, until eventually hydrogen burning has begun, and a zero age main sequence (ZAMS) zero metallicity star (a Pop. III star) is formed. For reviews of the standard picture of the formation Pop. III stars, see [1, 2, 3, 4, 5, 6, 7, 8].

One important question is the following: can Dark Matter heating alter the formation or evolution of the first stars? This was first addressed by [9], for the case of self-annihilating weakly interacting particles (WIMPs) dark matter, the most promising -on theoretical grounds- DM scenario at the time. The authors found that given three generic assumptions -see details below- the collapse of the protostellar cloud can be halted by dark matter heating. This leads to a new phase of stellar evolution, a “Dark Star (DS).” The three conditions identified by [9] that could lead to the formation of a DS are: (1) high DM density at the location of collapsing gas cloud, (2) some of the annihilation products can become trapped and thermalized inside the gas cloud/star, and (3) DM heating dominates over other heating or cooling mechanisms. Once the protostellar collapse has stopped, due to DM heating, a Dark Star is formed. This object will continue to accrete mass from the baryon cloud, and since it is typically cooler than a corresponding Pop. III star, it can become much more massive, as feedback effects that stop the accretion process are directly proportional to the surface temperature. The dark star phase has been subsequently confirmed by [10].

There are be two distinct sources for dark matter that can pile up and annihilate inside a star: (i) Adiabatic Contraction (AC) and (ii) capture. Adiabatic contraction [11] is a mechanism that is responsible for driving up the density of dark matter at the center of a collapsing gas cloud. Simply put, DM orbits shrink in response to the increase in the gravitational potential, and an increase in the baryon number density leads to an increase in the DM density. Capture is a mechanism via which dark matter particles can be trapped inside compact baryonic objects. For the case of WIMP dark matter, this phenomenon was studied in the literature when the baryonic object is: the Sun [12], the Earth [13, 14], a dark star [15, 16, 17], a neutron star or a white dwarf [18], etc. The case of strongly interacting superheavy DM capture by the sun was studied by [19]. For all of the above scenarios the physics is similar: if the orbit of a DM particle crosses the said object, there is a chance that collisions with the baryons will lead to the DM particle losing enough energy to become gravitationally trapped (i.e. whenever its speed becomes less than the escape velocity at the surface of the object).

Within the WIMP paradigm, the evolution and stellar structure of DSs has been extensively studied over the past decade. We summarize here some of the main results. Dark stars are born with masses $\sim 1M_{\odot},$ for a large variety of WIMP masses, and grow to much larger masses via accretion [20]. As alluded before, they are much puffier, and hence somewhat cooler than Pop. III stars, with radii of about $10~{}\,\mathrm{AU}$ and surface temperatures of $\sim 10^{4}~{}\,\mathrm{K}$. A dark star is essentially an object in hydrostatic equilibrium that is supported against gravitational collapse by the pressure due mostly to the heat deposited inside the star via the thermalization of products of DM-DM annihilations happening inside the star. At a certain phase during the evolution of a DS, the adiabatically contracted DM reservoir will be depleted by annihilations. If one assumes spherical DM halos, this will happen when the DS has annihilated about $1~{}M_{\odot}$ of DM. [20, 21] found that this is sufficient to allow the DS to reach, via accretion, a mass of $\sim 10^{3}M_{\odot},$ after about $300,000$ yr. At this stage a short contraction phase ensues, as now the DM annihilations are no longer efficient enough to prevent gravitational collapse. As the core temperature and density increase, nuclear fusion and/or DM capture can become efficient and prevent further collapse. Once a new equilibrium phase is reached, the DS enters the Zero Age Main Sequence (ZAMS).

N-body simulations indicate that DM halos are prolate-triaxial [22, 23, 24, 25], which leads to a large part of the orbits passing arbitrarily close to the center, i.e. centrophilic orbits [26, 27, 28]. As such, the AC DM reservoir can be refilled efficiently, and hence the DS phase prolonged significantly [29]. Within this “extended AC” phase, a DS could become supermassive, i.e. $M_{DS}\gtrsim 10^{6}M_{\odot}$. A separate mechanism that could lead to the formation of a Supermassive Dark Star (SMDS) is capture of Dark Matter, if the ambient DM density is high enough. SMDSs are very bright and could be observed with the James Webb Space Telescope (JWST), as shown by [29]. In view of their relatively cooler temperatures, compared to Pop. III stars, one could also distinguish them from early galaxies formed by Pop. III stars [30]. After the DM fuel runs out, a SMDS will quickly undergo a nuclear burning phase and then collapse to a Supermassive Black Hole (SMBH). We mention that in the standard picture, where the first stars are Pop. III stars, the formation of the SMDS required to power the very bright observed high redshift quasars222For example J1342+0928 [31] is a quasar at $z=7.45$ that is powered by a behemoth BH with a mass of $8\times 10^{8}M_{\odot}$ poses significant theoretical challenges. It requires either super-Eddington accretion or non-standard direct collapse to black holes (DCBH), two non-orthodox scenarios. In contrast, SMDS, and the SMBH they will form, are a natural solution for this problem [29, 32].

Most of the work on Dark Stars has been done under the assumption of a thermal WIMP DM particles, i.e. thermal relics with masses in the $\,\mathrm{GeV}-\,\mathrm{TeV}$ range, generated in matter-antimatter collisions in the early universe. Considering a variety of thermal DM models [33] finds that in general this kind of models lead to the formation of DSs. However, at the LHC no particle has been identified as dark matter, while at the same time, direct detection experiments such as LUX and XENON have ruled out a large swath of the preferred parameter space for thermal WIMP DM. This places very tight constraints on the realization of the standard thermal relic WIMP paradigm. Therefore, non-thermal models started to become more and more relevant. Of those we mention only two: at the lower end of the mass spectrum one has the axions [34], whereas WIMPZILLAs [35] are one example of superheavy dark matter (SHDM) particles. We note that there are many more DM particle models that are still within experimental bounds. For reviews on the state of particle dark matter, see [36, 37]. For a more general review on the state of dark matter in the universe, in view of all experimental constraints, see [38].

It is therefore important to revisit the DS picture when considering non-thermal DM models. In this work we analyze if the evolution of Pop. III stars could be altered by capture of SHDM particles with masses in the $10^{8}-10^{15}~{}\,\mathrm{GeV}$ range 333This range of masses for the Dark Matter candidate is representative for non-thermal relics produced during inflation, such as WIMPZILLAs [35]. Additionally, we assume that DM annihilations is not preventing a Pop. III star from forming and therefore we only consider its effects after the Pop. III star has reached the zero metallicity zero age main sequence and DM capture can become important. In a future study we plan to investigate under what conditions, if any, heating due to annihilations of adiabatically contracted superheavy dark matter could halt the collapse of protostellar gas cloud prior to ignition of nuclear fusion, and therefore lead to the formation of a Dark Star instead of a Pop. III star. Using the recently developed formalism for multiscatter capture of dark matter [39], we calculate the capture rates and associated luminosity due to captured SHDM annihilations inside Pop. III stars. Assuming a fiducial value of $\rho_{X}=10^{9}\,\mathrm{GeV/cm^{3}}$ for the ambient DM density, we find that the luminosity due to dark matter annihilations (DMAs) is typically orders of magnitude below the one due to nuclear fusion. For the most massive Pop. III stars, that shine at or near the Eddington limit, the additional heat source would quickly destabilize them and lead to an unbound object. Therefore one possible effect of DMA is to impose an upper-bound on the masses of the Pop. III stars. Since the DM heating is proportional to the ambient DM density, this upper bound is lower for higher values of $\rho_{X}$. As mentioned before the ambient dark matter density can be enhanced by many orders of magnitude near the center of a micro dark matter halo during the collapse of a protostellar gas cloud via a process called adiabatic contraction. For example, we find that, for an extreme case of $\rho_{X}=10^{18}\,\mathrm{GeV/cm^{3}}$, the Eddington limit places a bound on the masses of Pop. III stars of a few $M_{\odot}$. All of our results are upper bounds, obtained by using the assumption that the DM-nucleon scattering cross section ($\sigma_{n}$) has the highest possible value allowed by the exclusion limits from the latest direct detection experiments for strongly interacting massive particle (SIMP) dark matter obtained in [40].

This paper is organized as follows: in section 2 we review the physics of dark matter capture, focusing on the main ingredients in the analytical formalism of multiscatter capture of dark matter introduced by [39]. In section 3 we calculate the capture rates and associated luminosities due to SHDM annihilations inside Pop. III stars. Section 4 is dedicated to a discussion of one important effect of the additional heating from DM annihilations on Pop. III stars: a dark matter dependent upper bound on Pop. III stellar masses. We end with a discussion and conclusions in section 5.

2 Capture of Dark Matter: brief review

As it transits a gravitationally bound object (e.g. star, earth, dark star, neutron star, etc) of radius $R_{\star}$, a dark matter particle would have collided with its constituents (e.g. nuclei) an average number of times given by:

[TABLE]

Here $n_{T}$ is the average number density of collision targets (e.g. nuclei) inside the star, and $\sigma_{n}$ is the cross section of DM-nucleon scattering cross section. For the case of capture of weakly interacting dark matter particles with mass not much larger than $\sim 100~{}\,\mathrm{GeV}$ by the Sun [12] or the Earth [13, 14], the average number of such collisions is always less than unity, in view of the small scattering cross sections. For instance the XENON1T experiment places the most recent upper bounds for a $100~{}\,\mathrm{GeV}$ WIMP on the spin independent nucleon DM interaction cross section at the order $\sigma_{SI}\lesssim 10^{-46}~{}\,\mathrm{cm}^{2}$ [41], while for spin dependent cross section one has: $\sigma_{SD}\lesssim 10^{-41}~{}\,\mathrm{cm}^{2}$ [42]. Note that for a fixed capturing object mass ($M_{\star}$), the average number of collisions scales as $N\sim\sigma_{n}/R_{\star}^{2}$. Even for the case of the most compact objects, neutron stars, the average number of collisions is much less than unity, given the above bounds on the scattering cross sections. Therefore for WIMPs, one needs not consider the multiscatter capture formalism. Instead, in order to calculate the capture rates, one can use the single scattering capture formalism of [12, 13, 14].

Note that the experimental bounds via underground direct detection experiments are limited at large DM particle masses by the lower DM flux. Namely, the flux scales as: $\Phi\sim n_{X}\sim\rho_{X}/m_{X}\sim m_{X}^{-1}$. This leads to sensitivities that drop as $1/m_{X}$. In turn this amounts to constraints on the nucleon-DM scattering cross sections that scale as $\sigma_{bound}\sim m_{X}$, and therefore weaker bounds at high masses. [19] considers the regime of strongly interacting ($\sigma_{n}\sim 10^{-24}\,\mathrm{cm}^{2}$) superheavy ($m_{X}>10^{10}\,\mathrm{GeV}$) dark matter particles by the sun. In view of the large cross section, the mean free path of the DM particles inside the sun is very short ($\sim 10^{-2}\,\mathrm{cm}$), which amounts to a very large number of collisions each DM particle has, on average, as it traverses the Sun. After each collision, depending on the scattering angles, the factional energy loss by a DM particle is uniformly distributed in the following interval:

[TABLE]

In the above equation, $E_{i}$ represents the energy the dark matter particle of mass $m_{X}$ has before the $i$-th collision, $\Delta E_{i}$ is the energy it has lost due to the collision with a target nucleus of mass $m_{n}$, and $\beta_{\pm}\equiv 4m_{X}m_{n}/(m_{X}\pm m_{n})^{2}$. Therefore the average fractional energy loss per collision is, assuming equal probability for all values in the kinematically allowed range of equation 2.2, just half of the upperbound (i.e. $\beta_{+}/2$). The most efficient energy transfer happens when $m_{X}=m_{n}$, as expected. In the case of SHDM ($m_{X}\gg m_{n}$) one can see that the average fractional energy loss becomes $\Delta E/E\approx(1/2)(m_{n}/m_{X})\ll 1$. Again, this is to be expected. By analogy, a bowling ball would lose a very small fraction of its energy after a collision with a ping pong ball. Therefore, only DM particles at the low end of the velocity distribution can be captured in one collision. For most SHDM particles, it would take much more than one collision to slow down below the escape velocity. This is the case considered by [19], where analytic formulae for capture rates of a class of strongly interacting SHDM known as SIMPZILLAs can be found. However, experimental constraints on SHDM-nucleon scattering cross sections have closed the window on such DM models [43, 40], and as such the formalism developed in [19] has little practical use for DM capture physics. Capture of DM in the intermediary regime, where the average number of collisions it takes for a DM particle to be trapped is of order unity, can be treated using the formalism recently developed by [39]. This regime will be the one relevant for this work. In the reminder of this section we give a brief review of the main ingredients from [39] that we will subsequently use in this work.

There are two factors that control the rate of capture: the flux ($F$) of DM particles through the star and the probability ($\Omega$) that each DM particle crossing will lose enough energy via collisions such that it becomes bound by the gravitational field of the star. For all the regimes described above, the schematic differential capture rate is given as [12, 13, 14, 39]:

[TABLE]

where $u$ represents the DM velocity far from the gravitational field of the star, $v_{star}$ is the relative velocity of the star with respect to the halo, $v_{esc}^{halo}$ denotes the escape speed of the DM halo, $n_{T}(r)$ is the number density of scattering targets inside the star at a radius $r$ from the center, $w(r)^{2}=u^{2}+v_{esc}(r)^{2}$ is the velocity of the dark matter particles inside the star, with $v_{esc}(r)^{2}\equiv 2GM_{\star}(r)/r$. Since Pop III stars form usually at the center of DM halos, we can assume no relative motion of the star in the dark matter halo($v_{star}\rightarrow 0$). Additionally, the following assumptions are made: a very large escape speed for dark matter halo($v_{esc}^{halo}\rightarrow\infty$), a uniform density for the star, and a fixed escape velocity ($v_{esc}(r)=v_{esc}(R_{\star})$) where $R_{\star}$ is the radius of the star. The total capture rate corresponding to each dark matter particle mass($m_{X}$) is

[TABLE]

where $C_{N}$ is the capture rate after exactly N scattering events. The optical depth is conveniently defined as: $\tau\equiv n_{T}\sigma_{n}(2R_{\star})$, with $n_{T}$ representing the average number density of targets the DM particle can collide with as it traverses the star. In view of equation 2.1, this also represents the average number of collisions the DM particle will experience as it traverses the star. In addition to the DM flux, there are two separate factors at play when one calculates $C_{N}$: the probability that a DM particle with optical depth $\tau$ will actually participate in $N$ scatterings ($p_{N}(\tau)$), and the probability that the velocity of the DM particle drops below the escape speed at the surface of the star after exactly N collisions ($g_{N}(w)$). Taking into account the geometry of the different possible incidence angles, [39] defines $p_{N}(\tau)=(2/N!)\int_{0}^{1}\,dy\,ye^{-y\tau}(y\tau)^{N}$. We note that the integral can be performed analytically, with the following result which we will use throughout our calculations:

[TABLE]

We next consider $g_{N}(w)$. The initial kinetic energy of dark matter particle as it enters the stars is: $E_{0}=m_{X}w^{2}/2$. The energy loss after each elastic collision follows simple kinematics rules as $\bigtriangleup E=z\beta_{+}E_{0}$, where $z\in[0,1]$ is a kinematic variable related to the scattering angle, and $\beta_{\pm}\equiv 4m_{X}m_{n}/(m_{X}\pm m_{n})^{2}$. Therefore, the kinetic energy after one collision reduces to $E_{i}=(1-z_{i}\beta_{+})E_{i-1}$ and the velocity of the DM particle becomes $v_{i}=(1-z_{i}\beta_{+})^{1/2}v_{i-1}$. In our case, since $m_{X}\gg m_{n}$ we can approximate $\beta_{+}\approx 4m_{n}/m_{X}$. After $N$ scatters, the velocity and energy for a dark matter particle becomes[19, 39]:

[TABLE]

Considering the capture condition ($v_{N}<v_{esc}$) and taking into account the different possible paths a DM particle can trace inside the star while undergoing $N$ scatterings (i.e. different scattering angles for each collision), one can express the probability $g_{N}(w)$ as [39]:

[TABLE]

The $\Theta\left(v_{esc}\prod_{i=1}^{N}(1-z_{i}\beta_{+})^{-1/2}-w\right)$ factor amounts to the probability of capture after a path through the star described by $N$ collisions, each with scattering angle determined by $z_{i}$. For the case of Spin Independent cross sections, considered in this work, one can further assume no preferred scattering direction (i.e. cross section for nucleon-DM scattering is independent of scattering angle). Therefore the average value of $z_{i}$ can be assumed to be equal to 1/2, leading to a simplified form for $g_{N}(w)$ [39]:

[TABLE]

In order to calculate the capture rate after exactly $N$, scattering events one needs to multiply the particle rate through the surface of the star with the probability that a DM particle will be captured after $N$ collisions (i.e. the product $p_{N}(\tau)*g_{N}(w)$). The capture rate after exactly $N$ collisions ($C_{N}$) can be expressed as the following phase space integral [39]:

[TABLE]

where $f(u)$ is the DM velocity distribution. The assumption that the star is not moving with respect the dark matter halo is implicit. The integral can be analytically evaluated assuming a Maxwellian distribution with an average speed $\overline{v}$ and using the probability $g_{N}(w)$ from equation 2.8. For the full result see [39]. Under the following two assumptions: (i) $v_{esc}\gg\bar{v}$ and (ii) $m_{X}\gg m_{n}$ (which are always valid for the cases we will consider in this work) the following simplified expression for $C_{N}$ can be obtained [39]:

[TABLE]

where $\rho_{X}$ is the DM density, and $M_{\star}$ and $R_{\star}$ are the mass and radius of the star, $m_{n}$ mass of nucleons in the star, $\overline{v}$ the dispersion velocity of DM, $v_{esc}$ the escape velocity at the surface of the star.

In the next section we will implement a numerical procedure to calculate the infinite series in equation 2.4 and use it for the case of SHDM being captured by Pop. III stars. It is worth emphasizing that a constant density star is also assumed. Relaxing this assumption, and considering a variable density star, leads to an increase in the capture rates, as pointed out in [39]. With this in mind, all the bounds we compute in this paper should be viewed as conservative. In particular, this applies to the upper bounds on the total capture rates from section 3 and the associated Pop. III stellar mass limits obtained in section 4.

3 Superheavy dark matter capture by Pop. III stars

In this section, we obtain upper bounds on capture rates and associated DM-DM annihilation luminosity in the case of superheavy dark matter capture by Pop. III stars with mass ranging from a few to a thousand solar masses. As mentioned in section 1, the first stars formed at the center of dark matter mini-halos ($M_{halo}=10^{5}-10^{6}M_{\odot}$) at redshifts of $z=10-50$. We assume that the density profile of the dark matter halo follows the Navarro-Frenk-White[44] profile:

[TABLE]

where $\rho_{halo}$ is the density of dark matter halo, $\rho_{0}$ the density at the center of the DM halo,$r$ the distance away from the center of the DM halo, and $r_{s}$ is the scale radius. And the central density $\rho_{0}$ is related to the critical density $\rho_{crit}$ at given redshift $z$, via

[TABLE]

where $c\equiv r_{vir}/r_{s}$ is the concentration parameter and $r_{vir}$ is the scale radius. From the virial theorem, the dispersion velocity $\overline{v}$ of DM inside the halo is

[TABLE]

where

[TABLE]

is the gravitational potential of the dark matter halo.

We adopt here the same parameters for the NFW profile describing the DM halo where the first stars form as those used in [15], who studied the capture of WIMPs by Pop. III stars. Namely, the following ranges: concentration parameter $c=1-10$, redshift $z=10-50$, and $r_{s}=15-100~{}\,\mathrm{pc}$. For the DM dispersion velocity one gets values ranging from $\overline{v}=1-15~{}\,\mathrm{km/s}$. Therefore we adopt the following fiducial value:

[TABLE]

For the DM density at the location of the star (i.e. the center of the DM halo), we assume $\rho_{X}=\rho_{0}=10^{9}~{}\,\mathrm{GeV/cm^{3}}$. From equation 2.10 using numerical values and substituting $p_{N}(\tau)$ from equation 2.5 we obtain the following scaling relation:

[TABLE]

with the following approximate numerical scaling relationships for the average collision number ($\tau$) and the exponential factor ($A_{N}^{2}$):

[TABLE]

Note that one could use mass-radius “homology” relationships to eliminate the $R_{*}$ dependence in equations 3.6 and 3.8. For example, for massive ($M*\gtrsim 100M_{\odot}$) Zero Age Main Sequence (ZAMS) stars of very low metallicity, $R_{\star}\propto M_{\star}^{5/11}$ [45]. Also, see figure 2 and equations 3.12 and 3.13, where we obtain the homology mass-radius relations for our set of Pop. III stars.

At first glance it seems that the capture rate in equation 3.6 is proportional to $\sigma_{n}^{-2}$, which is contrary to expectation. A larger scattering cross section should, in principle, lead to a larger capture probability, and hence a larger rate of particles captured. This is certainly the case for the single scatter capture, where $C\propto\sigma_{n}$. The $\sigma_{n}^{-2}$ factor comes from $p_{N}(\tau)\propto 1/\tau^{2}\propto 1/\sigma_{n}^{2}$. However, $p_{N}(\tau)$ contains the additional factor: $\left(N+1-\dfrac{\Gamma(N+2,\tau)}{N!}\right)$, which also has a $\tau$ (and hence $\sigma_{n}$) dependence. It is straightforward to show that in the case of $\tau\ll 1$:

[TABLE]

As expected, for $\tau\ll 1$, which is to say when the dark matter particle has a very large (compared to $R_{\star}$) mean free path through the star, $C_{N}$ will quickly vanish for $N>1$. So, only $C_{1}$ is relevant, or, in other words the multiscatter formalism, in the appropriate limit, recovers naturally the single scatter formalism of [12, 13, 14]. Note that for low mass (WIMP-like) dark matter capture, the exponential factor $A_{N}^{2}\gg 1$ and the entire term in the square brackets of equation 3.6 reduces to one [15]. Therefore, isolating the scaling of $C_{1}$ with DM parameters we get: $C_{1}\propto\sigma_{n}\rho_{X}/m_{X},$ as expected and predicted by the single scatter capture analytic formulae of [14].

For the remainder of this work, we will adopt the XENON1T 2018 bounds [41] on spin independent dark mater nucleon scattering cross section. Since the first stars are mostly ionized hydrogen (nucleons), this bound applies to the cross section of the scattering of DM with the targets inside the first stars. Therefore, all of our results on the capture rates and associate quantities such as DM-DM luminosity are upper bounds. Specifically, by fitting the XENON1T one year exposure upper limits in the $10^{2}\,\mathrm{GeV}\leq m_{X}\leq 10^{15}\,\mathrm{GeV}$, we get the following linear dependence:

[TABLE]

We emphasize that this dependence is not due to any particle physics interaction; it is just a reflection of the fact that at large mass the sensitivity of direct detection experiments are limited by the DM flux. For theoretical expectations of the $\sigma_{n}$ dependence with $m_{X}$ for a large variety of DM models see [46, 47, 48].

In the calculation of the capture rates we will need parameters -radii and masses- for the stars of interest (i.e. Pop. III stars). As those are not yet observed, and it is unlikely even JWST will be able to observe them in isolation, so one has to resort to numerical simulations. We adopt models from [10], for stars with mass in the $5~{}M_{\odot}-600~{}M_{\odot}$ range. In addition we include a $1000M_{\odot}$ Pop. III star [49, 50]. We have checked that all models we considered are consistent with more recent simulations, such as those given in [51], for example. See table 1 for a summary of the parameters of interest for the Pop. III models we use.

Our first aim is to calculate $C_{tot}(m_{X})$, as given by equations 2.4 and 2.10. In principle, the series defining the total capture rate is infinite, but in view of $p_{N}$ decreasing for large $N$, the series converges after a number of steps $N_{cut}$ that is $\tau$ dependent. Defining the partial sum $C_{tot,N}\equiv\sum_{i=1}^{N}C_{i}$, our convergence criterion is the following:

[TABLE]

We code this criterion when estimating numerically the total capture rate $C_{tot}$ for each case we consider. It is important to note that for all the Pop. III models we considered $A_{N}^{2}\ll 1$, even for the largest values of $N$. In view of equation 3.7, this is due to the large values of $m_{X}$, and the fact that $p_{N}(\tau)$ drops to zero quickly after $N>\tau$, so values of $N\gtrsim 1000$ are never needed. In fact, the largest value for $N_{cutoff}$ is $818$, corresponding to the capture of $10^{15}\,\mathrm{GeV}$ dark matter particles by a $20M_{\odot}$ Pop. III star. Note that in the $A_{N}^{2}\ll 1$ and $\overline{v}\ll v_{esc}$ limits the term in the square brackets of the equations defining $C_{N}$ can be Taylor expanded. Keeping only the leading order terms we find:

[TABLE]

When taking this approximation into account, $C_{N}$ becomes:

[TABLE]

Note here one important distinction, when compared to the low mass case (WIMPs), where $C_{1}\propto\sigma_{n}\rho_{X}/m_{X}$. For SHDM, in the case of $\tau\ll 1$, using the expansion from equation 3.9 we get $C_{1}(SHDM)\approx\sigma_{n}\rho_{X}/m_{X}^{2}$. The extra factor of $1/m_{X}$ suppression comes from the different way the term in the square brackets of equation 2.9 scales with mass: $\propto m_{X}^{0}$ for WIMPs vs. $\propto 1/m_{X}$ for SHDM, as previously explained.

In figure 1 we plot the coefficients $C_{N}$ (capture rate after exactly $N$ scatterings) calculated for the case of $m_{X}=10^{14}\,\mathrm{GeV}$ as a function of $N$ for several Pop. III models. As expected, as large $N$ the rate $C_{N}$ drops significantly, approaching zero very fast after a certain threshold ($N_{cutoff}$). This is due to the decrease in the probability of capture if the number of scatterings ($N$) exceeds by a significant margin the average number of scatterings ($\tau$). In addition we note that the maximum value for $C_{N}$ increases with an increase in the mass of the star. This can be explained in the following way: $C_{N}\propto R_{\star}^{4}$ (see equation 3.11). Assuming a homology scaling relation $R_{\star}\propto M_{\star}^{p}$ we get $C_{N}\propto M_{\star}^{4p}$. The exponent is positive for any values of $p>0$. Numerically we find that the mass radius dependence can be well fit in the following way (see figure 2):

[TABLE]

Going back to figure 1, we point out that for the higher mass stars ($M_{\star}\gtrsim 20M_{\odot}$) the value of $N_{cutoff}$ has a mild inverse dependence with the mass of the star, i.e. higher mass stars have lower $N_{cutoff}$. This is to be expected, as

[TABLE]

For the higher mass stars we expect $N_{cutoff}\propto M_{\star}^{-0.12}$ in view of the numerically obtained value of $p=0.56$. We also confirmed this result using numerical fitting. For the lower mass stars, based on the same type of analysis, we have $N_{cutoff}\propto M_{\star}^{0.6}$.

Next we investigate the dependence of $N_{cutoff}$ with the mass of the dark matter particle. We expect a transition to happen whenever $\tau(m_{X})\sim 1$ is reached. The exact value of $m_{X}$ where this transition takes place depends on the mass and radius of the star (see equation 3.8), via the following combination $M_{\star}/R_{\star}^{2}$. As pointed out before, for Pop. III stars with $M_{\star}\gtrsim 20M_{\odot}$, this converts into a very mild dependence with mass of the star: $M_{\star}^{-0.12}$. For all cases considered, this transition happens when $m_{X}\sim 10^{12}\,\mathrm{GeV}$. For lower $m_{X}$, when $\tau\leq 1$, we recover the single scattering formalism, and therefore $N_{cutoff}\approx 1$. As $\tau\propto\sigma_{n}$, and the upper bound on the scattering cross section scales linearly with $m_{X}$, we expect a linear dependence of $N_{cutoff}$ with $m_{X}$. All those trends can be confirmed in figure 3.

We next do two consistency checks, to reproduce results previously published in the literature. First, in the case of WIMP (i.e. $m_{X}\sim 100\,\mathrm{GeV}$), and for constant scattering cross section, the capture can be calculated using the single scattering formalism, and rate scales as: $C\propto 1/m_{X}$ [15]. We have numerically checked that the multi scattering formalism reproduces similar values for the capture rates for the same Pop. III stars as those considered in [15] for capture of WIMPs. We additionally checked numerically that for the case of strongly interacting SHDM of [19] our implementation of the multiscattering capture formalism reproduces the results calculated analytically under the $\tau\gg 1$ assumption in [19]. After those consistency checks we proceed to calculate the upper bounds on the total capture rates for SHDM of mass in the $10^{8}\,\mathrm{GeV}-10^{15}\,\mathrm{GeV}$ by Pop. III stars (see figure 4).

It is remarkable that for the entire range of dark matter mass considered, the upper bound on the capture rate scales as $C_{tot}\propto m_{X}^{-1}$. As one can see from figure 3, for all stars considered, at masses below $\sim 10^{12}\,\mathrm{GeV}$ the single scattering formalism holds (i.e. $N_{cutoff}=1$). Therefore, as explained in the discussion below equation 3.11, $C_{1}\propto\sigma_{n}\rho_{X}/m_{X}^{2}$. The upper bound on the scattering cross section scales linearly with $m_{X}$(see equation 3.10). Therefore $C_{1}\propto\rho_{X}/m_{X}$, as evidenced by trend at the lower mass range in figure 4. At the higher mass range, when $\tau>1$ we now have to consider adding multiple terms that each scale, in the large N limit, approximately like: $C_{N}\propto N^{2}\rho_{X}/(\sigma_{n}m_{X})^{2}$. Note that we have dropped the $\Gamma(N+2,\tau)/N!$ term, since this term is approximately zero whenever $N<\tau$ (i.e. for $N<N_{cutoff}$). Therefore:

[TABLE]

For the masses considered ($m_{X}\gtrsim 10^{12}\,\mathrm{GeV}$) we expect $N_{cutoff}\propto\tau\propto\sigma_{n}$. This conclusion is reinforced from fitting the data in figure 3, as alluded to before. Therefore, we find $C_{tot}\propto\sigma_{n}\rho_{X}/m_{X}^{2}$ in this regime as well. The fact that the two regimes, single scatter ($\tau\lesssim 1$) and multiple scatter ( $\tau\gtrsim 1$), have the same scaling for the total capture rate can be traced back to the fact that $A_{N}^{2}\ll 1$ for both cases, and therefore the approximation for the term in the square brackets used when obtaining equation 3.11 holds for both cases. From equation 3.7 and table 1 we estimate that the transition to $A_{N}^{2}\gtrsim 1$ happens somewhere in the $10^{4}-10^{6}\,\mathrm{GeV}$ range for all the Pop. III stars considered, with the lower bound corresponding to the $5M_{\odot}$ star, and the upper bound to the $1000M_{\odot}$ case. 444In estimating the $m_{X}$ values quoted we kept $\overline{v}=10\,\mathrm{km/s}$. Furthermore, for this mass range $\tau\ll 1$, so effectively $N=1$ (i.e. we only have to consider single scattering). When $A_{1}\sim 1$ the scaling of $C_{1}$ and, hence of $C_{tot}$ with the physical parameters of interest is more complicated. At even smaller $m_{X}$, when $A_{N}^{2}\ll 1$, one recovers the expected scaling for WIMPs: $C_{tot}=C_{1}\propto\sigma_{n}\rho_{X}/m_{X}$.

We point out that for the case of $A_{N}^{2}\ll 1$, relevant to this work, the sum defining the total capture rate can be estimated analytically by substituting $(N+1-\Gamma(N+2,\tau)/N!)$ with $\tau^{2}p_{N}(\tau)/2$ in equation 3.11 and then noting that

[TABLE]

With this, and using the approximation for $C_{N}$ from equation 3.11 we get the following estimates for the total capture rate:

[TABLE]

For the scaling of the upper bound on the total capture rate, in view of $\sigma_{n}\propto m_{X}$ we find, as expected from figure 4, $C_{tot}\propto m_{X}^{-1}$. It is worth reminding the reader that here when calculating the upper bounds on $C_{N}$ we assume the scattering cross section at the limit allowed by the XENON1T one year exposure limits for spin independent (SI) DM-nucleon scattering. The upper limits on spin dependent (SD) scattering cross section are weaker by a few orders of magnitude, so if one were to use the SD scattering, our upper bounds on the total capture rate would be enhanced by the same factor.

We summarize here our most important results so far. For the case of SHDM with masses $m_{X}\gtrsim 10^{12}\,\mathrm{GeV}$, when multiscatter capture becomes important, the interplay between the $\tau$ dependence of $N_{cutoff}$ and the scaling of each term in the series with $N$ leads to a total capture rate that scales as: $C_{tot}\propto(\sigma_{n}\rho_{X}/m_{X}^{2})$. This is one of the most important, and perhaps counter intuitive results, of our work. Additionally, we find that for SHDM with $m_{X}\lesssim 10^{12}\,\mathrm{GeV}$ the multiscatter formalism is redundant, when considering Pop.III stars, and therefore $C_{tot}=C_{1}\propto\sigma_{n}\rho_{X}/m_{X}^{2}$. It is remarkable that the total capture rates has the same scaling in both regimes. Regarding upper bounds of the total capture rate, we note that for the entire SHDM mass range considered ($10^{8}-10^{15}\,\mathrm{GeV}$) we obtain: $C_{tot}\propto\rho_{X}/m_{X}$. We will investigate the possible consequences of this finding on the maximum mass of Pop. III stars in the next section.

4 Pop. III stellar mass limits

After DM particles get captured by any object, they can self-annihilate, assuming they are their own anti-partners. This could lead to a new source of energy for the star, and potentially could even disrupt its evolution or impose a cutoff on the maximum mass for a Pop. III star. For the case of WIMPs this was studied by [15]. We proceed here to do a similar analysis for the case of SHDM capture.

The number of dark matter particles in the star can be modeled by the competition between capture and annihilation, described mathematically by the differential equation below: 555We neglect evaporation effects, as those are only relevant for sub-$\,\mathrm{GeV}$ DM particle mass.

[TABLE]

where $C(s^{-1})$ is the capture rate, $\Gamma_{A}(s^{-1})$ the annihilation rate.

For WIMPs, this equilibrium is quickly reached within the lifetime of the star by a time scale $\tau$ [15]. We have checked numerically that, assuming a distribution where most of the DM captured is located near the core of the star, this conclusion holds in our case as well. At the equilibrium where the rate of change of the numbers of the dark matter particles in the star equals to zero, one has the following relationship between capture and annihilation:

[TABLE]

At this time the DM in the core can provide an additional stabilizing energy source, with a total luminosity :

[TABLE]

since an energy of $2m_{X}$ is released after each annihilation event. Henceforth $f$ represents the efficiency of energy conversion from dark matter particle to power up the star, for which we will assume, following [9], a value of $2/3$. This convention amounts to one third of the annihilation energy being lost to neutrinos, and the other two thirds being deposited in the star. In figure 5 we plot the luminosity due to dark matter annihilations, assuming equilibrium between capture an annihilations. We mention in passing that if this assumption is dropped, one would need to develop a Monte Carlo simulation, tracking each DM particle, as it gets captured and annihilated.

We next proceed to estimate an upper bound on the mass of Pop. III stars, by using the Eddington limit. As a star becomes purely radiation pressure dominated, its luminosity will become linearly proportional with mass, in what is known as the Eddington limit. For a given stellar mass $M_{\star}$, no star can shine brighter than the Eddington luminosity, as any further accretion of mass is disrupted by the radiation pressure. The Eddington luminosity is defined as:

[TABLE]

where $G$ is the gravitation constant, $c$ is the speed of light, $M_{\star}$ is the mass of the star and $\kappa_{\rho}$ is the opacity of the stellar atmosphere. Since the first stars are metal free with hot atmosphere, the opacity is mainly due to Thompson electron scattering. For ease of comparison with previous results from the literature for capture of WIMPs [15], we adopt the same opacity, and hence the same value for the Eddington limit:

[TABLE]

In reality, for Thompson electron scattering $\kappa$ depends only on the hydrogen fraction (X) and takes a particularly simple form: $\kappa=\kappa_{es}=0.2(1+X)\,\mathrm{cm^{2}s^{-1}}$. For any stars the value of X will change, as the star burns H into He, and thus even while on the main sequence $X$ is not a constant. However, for $X$ between $1-0.5$ we get the value of the prefactor in equation 4.5: ranging between $3.25-4.33$, therefore for ZAMS Pop. III stars, to a good approximation the Eddington limit is represented by the mass dependent value in equation 4.5.

We use two criteria for finding an approximation for the maximum mass of Pop. III stars, when including the effects of heating from captured DM-DM annihilations:

[TABLE]

For $L_{nuc}$ we take the values of the stellar luminosities for Pop. III stars, where no DM heating is included, as tabulated in table 1. The first criterion will lead to somewhat more stringent bounds; however, it assumes that the rate of hydrogen burning will not be affected by the additional heat source from DM-DM annihilations. The second criterion is very strict, in the sense that it completely disregards any energy from nuclear fusion. In reality, the actual value for $M_{max}$ will be somewhere in between the two bounds, closer to the one obtained using the first criterion. For a more precise estimation one would need to run a stellar evolution code with the effect of DM heating annihilation included. However, for our purposes, where an order of magnitude estimate of the upper limit of the stellar mass is what we are looking for, this approach would be overkill.

In figure 6 we plot the relevant luminosities as a function of stellar mass. We note that $L_{DM}$ can be very well approximated with a broken power law. In fact, combining the estimate for the total capture rate from equation 3.16 with the luminosity due to dark matter, when the equilibrium between capture and annihilation has been reached: $L_{DM}=fC_{tot}m_{X}$, we get:

[TABLE]

Using the homology scaling relation for the higher mass Pop. III stars found in equation 3.13 one gets $L_{DM}\propto M_{\star}^{1.88}$. This is a remarkable result and is in contrast to the scaling obtained for WIMP capture by [15], where $L_{WIMP}\propto M_{\star}^{1.55}$. For SHDM capture, the increase in the DM luminosity with the mass of the star is significantly faster than the one for WIMPs. As explained before, when discussing the scaling of the total capture rate, the differences between the scaling relationships for WIMPs vs SHDM can be attributed to the $A_{N}^{2}\ll 1$ (WIMPs) vs. $A_{N}^{2}\gg 1$ (SHDM), and therefore the different approximations valid for the term in the square brackets of equation 2.10.

For the nuclear fusion luminosity of ZAMS zero metallicity stars (Pop. III stars) we find the following approximate fitting formula:

[TABLE]

where $x\equiv\frac{M_{\star}}{M_{\odot}}$ and $L_{\odot}\equiv 3.846\times 10^{33}\,\mathrm{erg/s}$. This interpolates between the expected scaling regimes $L_{nuc}\propto M_{\star}^{3}$ for lower mass stars, and $L_{nuc}\propto M_{\star}$, for high mass stars, which are radiation pressure dominated and shine at the Eddington limit.

Since $L_{DM}\propto\rho_{X}$, an increase in the DM ambient density will lead to a lower value of $M_{max}$. One mechanism identified in the literature that can lead to an increase in $\rho_{X}$ at the center of a DM halo during the formation of a star is adiabatic contraction [11]. For instance, [9] shows that during the formation of the first stars adiabatic contraction can lead to an enhancement of $\rho_{X}$ by many orders of magnitude. If adiabatic contraction operates until the protostellar core gas has a number density $n\sim 10^{22}\,\mathrm{cm^{-3}}$, the DM density would reach $\rho_{X}\sim 10^{18}\,\mathrm{GeV/cm^{3}}$. x

In figure 7 we plot the maximum Pop. III stellar mass versus $\rho_{X}$. The kink at $\rho_{X}\approx 10^{16}\,\mathrm{GeV/cm^{3}}$ is due to the different radius-mass homology relations for low vs high mass Pop. III stars. Note that if the DM ambient density reaches $\rho_{X}\sim 10^{14}\,\mathrm{GeV/cm^{3}}$ the effects of captured SHDM annihilations are able to prevent stars with masses larger than roughly $1000M_{\odot}$. At lower densities the maximum stellar mass is likely to be controlled by other effects, such as fragmentation of the protostellar gas cloud. However, when $\rho_{X}\gtrsim 10^{14}\,\mathrm{GeV/cm^{3}}$, SHDM capture can place very stringent bounds on the maximum mass a Pop. III star can have. For the extreme case of $\rho_{X}\approx 10^{18}\,\mathrm{GeV/cm^{3}}$ the maximum Pop. III mass becomes roughly one solar mass!

In the case of the less stringent criterion for $M_{max}$ we find, via equations 4.5 and 4.8, the following approximation:

[TABLE]

If the DM ambient density is not higher than roughly $10^{16}\,\mathrm{GeV/cm^{2}}$ we can use the mass-radius scaling relations of equation 3.13, valid for $M_{\star}\gtrsim 20M_{\odot}$. For higher DM densities $M_{max}\lesssim 20M_{\odot}$, so one will need to use the mass-radius scaling relations of equation 3.12. Assuming the cross section given by the upper bound of the XENON1T one year exposure (see equation 3.10) we find:

[TABLE]

We confirmed this result using a numerical fit of the relevant points from figure 7. Including the effect of the nuclear fusion luminosity when calculating $M_{max}$, leads to a milder dependence of $M_{max}$ with $\rho_{X}$. By fitting with a power law we find $M_{max}\propto\rho_{X}^{-0.81}$ (when $\rho_{X}\lesssim 10^{16}\,\mathrm{GeV}$) and $M_{max}\propto\rho_{X}^{-0.63}$ (when $\rho_{X}\gtrsim 10^{16}\,\mathrm{GeV}$). From figure 7 note that the two criteria are essentially indistinguishable at $\rho_{X}\gtrsim 10^{16}\,\mathrm{GeV/cm^{3}}$. This is due to the fact that at those high DM densities the luminosity due to DM annihilations dominates over the one due to fusion for the stars with masses close to the respective upper bounds.

In conclusion, annihilation of SHDM captured by Pop. III stars can lead to a maximum value of the stellar mass that is inversely proportional to the ambient DM density $\rho_{X}$. This confirms results previously obtained in [15] for WIMPs and extends them to the capture of SHDM. If the DM ambient density is as high as $10^{18}GeV/cm^{3}$, i.e. if we assume adiabatic contraction throughout the entire history of the protostellar gas cloud collapse, we find that heating due to annihilations of captured SHDM prevents Pop. III stars from having masses larger than roughly a solar mass! We mention in passing that this argument can also be turned around and used to constrain the following combination of dark matter parameters: $(\rho_{X}\sigma_{n}/m_{X})$, since the dark matter heating is directly proportional to this combination of parameters. Namely, for any future observed Pop. III star of $M_{obs}$, one can rule out the combination $(\rho_{X}\sigma_{n}/m_{X})$ that would lead to $L_{DM}>L_{Edd}(M_{obs})$. This opens up the possibility to place constraints on properties of dark matter from the observed masses of the first stars, which should be possible with JWST.

5 Conclusion

In this paper we have investigated the capture of SHDM by Pop. III stars. Using the multi scattering formalism, we find, for SHDM with masses in the $10^{8}-10^{16}\,\mathrm{GeV}$ range, that the total capture rate has the following scaling with the relevant parameters: $(\rho_{X}\sigma_{n}/m_{X}^{2})(M_{\star}^{3}/R_{\star}^{2})$. This leads to an upper bound on the capture rate that is proportional with $1/m_{X}$, which, in turn implies a $m_{X}$ independent upper bound on the luminosity due to captured SHDM annihilations. We also show that if a non-zero fraction of the annihilation products can thermalize and deposit energy inside the star, the masses of Pop. III stars have an upper bound ($M_{max}$) that is inversely proportional to the ambient dark matter density. We estimate values for $M_{max}$ for ambient DM densities ranging from $10^{9}-10^{18}\,\mathrm{GeV/cm^{3}}$. At $\rho_{X}\sim 10^{16}\,\mathrm{GeV/cm^{3}}$, and assuming the scattering cross section given by the XENON1T one year detection bounds for one year exposure, we obtain a maximum value of the mass of a Pop. III star of roughly a few tens of solar masses. In the extreme case of $\rho_{X}\sim 10^{18}\,\mathrm{GeV/cm^{3}}$ Pop. III stars could not grow much past a few solar masses, in view of the Eddington limit. The existence of a dark matter dependent upper limit on Pop. III stellar masses raises the possibility to constrain properties of dark matter by using upcoming observations of first stars with JWST, and their inferred masses.

Appendix A Captured SHDM Annihilations

Dark Matter captured by any star will settle to the core on a time scale that is determined by their drift velocity. Specifically, $\tau_{drift}\sim r/v_{drift}$, with the drift velocity being found by imposing the equilibrium between the gravitational force and the viscous drag force due to scatterings with nucleons [52, 19]:

[TABLE]

with $\bar{\rho}$ and $\bar{T}$ being the typical mass density and temperature of the star. Since we are only interested in an order of magnitude estimate we can set $\bar{\rho}\sim n_{N}m_{N}$, leading to $\tau\propto\sigma_{N}m_{X}^{-1}(\bar{T})^{1/2}$. For all the DM masses considered in our paper $\tau_{drift}$ is very short, much less than the lifetime of the star. Therefore, quickly after being captured SHDM particles sink to the core of the Pop. III star, where they could annihilate.

We next detail the steps taken to estimate the time scale after which an equilibrium between capture of supermassive dark matter by a compact object and its subsequent annihilations is attained. We are following the same procedure as [15] did for the case of WIMPs being captured. The main difference comes from the annihilation cross section, which for thermal relics, such as WIMPs, is set by the dark matter relic abundance. For the case of non-thermal relics, such as the SHDM considered in this paper, the annihilation cross section is not tied to the relic abundance. However, the unitarity bound restricts the annihilation cross section:

[TABLE]

The number of dark matter particles in the star can be modeled by the competition between capture and annihilation, described mathematically by the differential equation below:

[TABLE]

where $C(s^{-1})$ is the capture rate, and $\Gamma_{A}(s^{-1})$ is the annihilation rate,

[TABLE]

is defined as an $N-$independent annihilation coefficient. Solving this equation, we find

[TABLE]

At the equilibrium between annihilations and capture, when the rate of change of the numbers of the dark matter particles in the star equals to zero, one has the following relationship between capture and annihilation:

[TABLE]

Note that for $t\gtrsim\tau$, the annihilation rate of equation A.5 becomes time independent, as the hyperbolic tangent approaches unity and $\Gamma_{A}\to 1/2C$. Therefore $\tau$ corresponds to the time scale for annihilation/capture equilibrium.

In general, the annihilation rate is given by the following integral over the volume of the star: $\Gamma_{A}=\int\,dVn^{2}_{X}(r)\langle\sigma v\rangle_{\mathrm{ann}}$. In view of our discussion following equation A.1 we can assume that most of the DM captured is located near the core of the star. Therefore, the DM density profile can be approximated using isothermal distribution

[TABLE]

where $n_{c}$ the central number density of DM and $T_{c}$ is the central temperature of the star,

[TABLE]

is the gravitational potential at radius $r$ with respect to the center, and $M(r)$ is the mass interior to $r$. For order of magnitude estimates we assume a constant density profile and therefore $M(r)\sim\rho_{c}4\pi/3r^{3}$, which leads to the following form for the distribution of DM particles [53]: $n_{X}(r)\propto\exp(-r^{2}/2r_{X}^{2})$. Note that $r_{X}$ roughly determines the size of the core within which the capture DM particle are concentrated and it can be expressed as:

[TABLE]

Inserting numerical values we can estimate the following scaling relation:

[TABLE]

which shows the size of the DM core within which annihilations could happen efficiently is very small, compared to the full radius of the star. This holds for the entire range of parameters in both DM mass and stellar temperature and densities considered here.

We next proceed to calculate the time scale $\tau$ that sets the equilibrium between capture and annihilations. One can define effective volumes

[TABLE]

where $m_{pl}$ is the Planck mass, and $\rho_{c}$ is the core mass density of the star. The name “effective volume” is suggestive since we have $N=n_{0}V_{1}$. The total annihilation rate is given as

[TABLE]

where the annihilation cross section

[TABLE]

is assumed to be at the unitarity bound. The $N-$independent annihilation coefficient becomes:

[TABLE]

With this expression for $C_{A}$ we evaluate the equilibrium timescale $\tau$, defined in equation A.5. For the capture rate $C$ we used our approximation from equation 3.16. For the stellar and dark matter parameters considered in our work this time scale ranges from $10^{-3}\,\mathrm{yr}$ to $10\,\mathrm{yrs}$. This is well within the lifetime of the star, and therefore we can assume an equilibrium between capture and annihilation rates. Therefore, the luminosity due to annihilations of captured dark matter can be related to the capture rate: $L_{DM}=f\Gamma_{A}2m_{X}=fCm_{X}$, as we did in section 4.