Dark Matter and General Relativistic Instability in Supermassive Stars

Abstract

We calculate the extent to which collisionless dark matter impacts the stability of supermassive stars $(M\gtrsim 10^{4}\,M_{\odot})$. We find that, depending on the star’s mass, a dark matter content in excess of ${\sim}1\%$ by mass throughout the entire star can raise the critical central density for the onset general relativistic instability, in some cases by orders of magnitude. We consider implications of this effect for the onset of nuclear burning and significant neutrino energy losses.

[lines=5]In this paper, we examine how dark matter could impact the onset of general relativistic instability in massive self-gravitating configurations, i.e., supermassive stars (SMSs), as well as their nuclear and neutrino evolution. SMSs are supported principally by radiation pressure, making them vulnerable to general relativistic instability. If these extreme objects had existed in the early universe, and then collapsed to black holes, they could provide massive seeds that facilitate early production of supermassive black holes (SMBHs).

Indeed, SMBHs seem to be extant at an embarrassingly early epoch in the history of the universe. For example, recent high redshift surveys from the James Webb Space Telescope (JWST) [1] and even the Chandra X-Ray Observatory [2] suggest the existence of compact objects with masses ${\gtrsim}10^{7}\,M_{\odot}$ at redshift $z\sim 10$, when the universe is a scant ${\sim}400\,{\rm Myr}$ in age. These observations are consistent with lower redshift observations (e.g. [3, 4, 5, 6, 7, 8]) that also suggest a significant inventory of SMBHs early on. Moreover, these black holes may be growing too quickly to be consistent with building up the obsevationally-inferred SMBH masses given a ${\sim}1\text{--}10\,M_{\odot}$ progenitor black hole, as expected in the standard pictures of stellar evolution and mass accretion.

By contrast, at the end of their lives, SMSs would likely collapse directly to black holes, providing large SMBH “seeds” that would be able to grow to the SMBH sizes observed in the early universe. Such scenarios are favored by some modelers since they do not invoke super-Eddington accretion rates to explain the high masses observed [9, 10]. Bogdán et al. [2] reports having observational evidence that supports this argument. The collapse of the SMS itself is potentially observable with future gravitational wave observatories like DECIGO [11]. This gravitational wave signal stems, in part, from an assumed (slightly) non-spherically symmetric neutrino burst generated from the collapsing homologous core [12].

For the purposes of this paper, an “SMS” is any hydrostatic and self-gravitating gas cloud with a mass in excess of ${\sim}10^{4}$ $M_{\odot}$. These are fully convective, high entropy systems whose pressure support is almost entirely derived from radiation (i.e. relativistic particles). They will eventually suffer general relativistic instability [13, 14, 15, 16]. Such configurations, those that in general rely on pressure from relativistic particles, are “trembling on the verge of instability,” and are most widely realized in compact scenarios like white dwarfs and neutron stars. In Newtonian gravitation, radiation dominated systems are nearly neutrally buoyant under radial perturbations. If the pressure is entirely derived from radiation, the outward pressure forces always perfectly match the inward gravitational forces as the object is expanded and contracted, meaning such perturbations cost no energy to perform, i.e., a metastable configuration.

The situation changes, however, when first order general relativistic (GR) corrections are included in the analysis. Being non-linear, GR predicts that both stress-energy (primarily baryon rest mass, in this case) and spacetime curvature itself cause local spacetime curvature. This makes the metastable system under Newtonian gravity actually marginally unstable under GR. This is the Feynman-Chandrasekhar post-Newtonian instability [17, 18].

Specifically, upon the star’s core exceeding a critical central density, it will become unstable to collapse and unable to establish hydrostatic equilibrium. It will most likely collapse directly into a black hole, but may explode completely if energy derived from nuclear reactions (or some other source) can be added quickly enough to make up for the in-fall kinetic energy of collapse and neutrino energy losses [19]. If the energy source in this scenario indeed stems from nuclear reactions, then an explosion is only possible if the star has sufficient metallicity to burn Hydrogen with the more efficient Carbon-Nitrogen-Oxygen (CNO) cycle instead of through the Weak Interaction limited proton-proton (PPI) chain [16].

While the formation of such SMSs seems dubious at first, especially when considering cloud fragmentation, Ilie et al. [20] argues that tantalizing evidence of such objects may have been observed by JWST. This, however, begs the question as to how they would have been observed in the first place since they are so close to instability. Freese et al. [21] puts forward the idea that these SMSs could be supported through particle dark matter self-annihilation into standard model particles which would keep them “puffy” (i.e. lower the central density) and prevent GR instability, at least until they exhaust their dark matter reserves.

Alternatively, a SMS also could be stabilized by dark matter without assuming such a self-annihilation cross section [22, 23]. If it is assumed that the Dark Matter particles are collisionless and only interact with the Standard Model component of the star through gravitation, then it is possible to derive how much its presence alone raises the critical central density beyond which the star suffers the Feynman-Chandrasekhar instability.

Non-intuitively, adding more matter to the meta-stable star, specifically a non-interacting non-relativistic fluid, actually tends to stabilize the configuration as opposed to tipping it closer to collapse. The physical mechanism of this increased stability arises from the dark matter acting as a source of non-relativistic stress.

Absent a significant dark matter content, heavier (${\gtrsim}10^{5}\,M_{\odot}$) SMSs will undergo the Feynman-Chandrasekhar instability prior to, or just after Hydrogen burning ignition, meaning they effectively have no main sequence phase (see e.g. Fuller et al. [16]). An interesting question is then whether a sufficient Dark Matter content in these objects could delay the onset of instability, allowing a hydrostatic burning phase, and thereby changing the downstream nuclear and entropy evolution of these stars.

In Section I, we establish the overall framework to analyze stability of SMSs by considering the case without any dark matter. We determine the critical central density from extremizing the total energy of the configuration with a 1 Order GR correction term to derive the standard result.

In Section II, we examine two kinematic limits of the dark matter based on its specific kinetic energy relative to the gravitational potential at the SMS’s surface, determine the total energy added by this dark matter, and derive how its presence impacts the critical central density as a function of the dark matter mass fraction. We then present our results for several SMS masses to examine how the effect also depends on the total mass.

In Section III, we approximate how the total fraction of dark matter changes with radial perturbations of the SMS, providing insight into the physical mechanism for the enhanced stability discussed in the previous section.

Finally, in Section IV, we investigate nuclear burning rates, neutrino emission rates, and their respective luminosities as a fraction of the stars’ Eddington luminosities. This allows us to ascertain how “important” nuclear burning is to an SMS’s evolution and when neutrino emission overpowers it. We also determine by how much these luminosities change with the addition of dark matter, and what that may mean for the downstream stability and entropy evolution.

I Stability of Supermassive Stars

The analysis in this section will closely follow Butler et al. [24]’s analytical treatment, specifically their “Approach II” in Section 2.5, to derive the onset of stability for SMSs. The SMSs of concern in this paper are taken to be spherically symmetric, with primordial metallicity, fully ionized and fully convective. They are high entropy gas clouds supported against gravity almost entirely by radiation pressure $P_{\text{rad.}}$, with a marginal but important contribution from classical Boltzmann gas pressure $P_{\text{gas}}$. We define the parameter $\beta$ as the ratio of gas pressure to total pressure (the final approximation following from Fowler [25])

$$\beta=\frac{P_{\text{gas}}}{P_{\text{gas}}+P_{\text{rad.}}}\approx\frac{4.3}{\mu}\left(\frac{M}{M_{\odot}}\right)^{-1/2},$$ (1)

where $\mu=0.59$ is the mean molecular weight (amu) per particle for a primordial composition gas.

Given our approximations, the equation of state (EOS) for this configuration can be modeled accurately as a polytrope:

$$P(s,\rho)=K(s)\rho^{\Gamma}=K(s)\rho^{1+1/n},$$ (2)

where $P$ is the total pressure, $s$ is the entropy per baryon, $\rho$ is the mass density, and $\Gamma$ is the so-called polytropic exponent, with the associated quantity $n$ being the polytropic index. For a radiation dominated system, we can safely assume that the radiation field dominates the entropy per baryon $s\approx s_{\text{rad.}}=\tfrac{4}{3}aT^{3}/n_{\text{b}}$, where $T$ is the plasma temperature, $n_{\text{b}}=\rho/(\mu m_{\text{b}})$ is the baryon number density and $m_{\text{b}}\approx 931.5$ MeV is the mass of an atomic mass unit. Here, in natural units, $a=\pi^{2}/15$ is the radiation constant. Since $s$ is constant throughout the fully convective SMS, we can use Eq. (1) and $P_{\text{gas}}=n_{\text{b}}T$ and $P_{\text{rad.}}=aT^{4}/3$ to derive

$$s\approx\frac{4\mu}{4.3}\left(\frac{M}{M_{\odot}}\right)^{1/2}-4,$$ (3)

where here we express $s$ in units of Boltzmann’s constant $k_{\text{B}}$ per baryon. Note that for the SMSs considered here the entropy per baryon is large, e.g., $s\sim 1000$ for $M\sim 10^{6}\,M_{\odot}$. This large entropy implies a large number of photons per baryon, resulting in implications for the $e^{\pm}$ content of the star during its eventual collapse which, in turn, has implications for neutrino emissivity [26, 27, 28].

Using $P(s,\rho)=P_{\text{gas}}+P_{\text{rad.}}$ and $n=3$, $K(s)$ is found to be

$$K(s) =\frac{1}{1-\beta}\frac{a}{3}\left(\frac{3}{4a}\frac{s}{\mu m_{\text{b}}}\right)^{4/3},$$

$$\approx\frac{0.5801}{1-\beta}\left(\frac{0.59}{\mu}\right)^{4/3}\left(\frac{s}{1000}\right)^{4/3}\,\rm{MeV}^{-4/3}.$$ (4)

For the stars discussed in this paper, their EOS and structure are well approximated by an $n=3$ ($\Gamma=4/3$) polytropic model.

Since these stars are fully convective, we can take their entropy per baryon to be constant throughout. With constant entropy, the polytropic exponent $\Gamma=1+1/n$ is equivalent to Chandrasekhar’s adiabatic index $\Gamma_{1}$, defined as [29]

$$\Gamma_{1}\equiv\left(\frac{\partial\log P}{\partial\log\rho}\right)\Big{|}_{s}.$$ (5)

In the purely Newtonian gravitation case, the frequency of oscillations in response to radial perturbations is proportional to $(\left<\Gamma_{1}\right>_{P}-4/3)^{1/2}$, where the pressure-averaged value of $\Gamma_{1}$ throughout the star is $\left<\Gamma_{1}\right>_{P}$. Therefore, in order for oscillations to have a real frequency and thereby stable against perturbations, $\left<\Gamma_{1}\right>_{P}$ must be strictly greater than $4/3$. When including 1 Order GR corrections, however, this condition of stability changes to

$$\left<\Gamma_{1}\right>_{P}>\frac{4}{3}+\mathcal{O}\left(\frac{R_{\text{S}}}{R}\right),$$ (6)

where $R$ is the radius of the star and $R_{\text{S}}=2GM$ is the star’s Schwarzschild radius. For the types of stars considered in this paper, this correction is typically no more than one part in ten-thousand.

As the star quasistatically contracts, however, $\left<\Gamma_{1}\right>_{P}$ asymptotically approaches $4/3$ from above, meaning that the adiabatic index will eventually fall below this threshold and make the whole configuration unstable. This point highlights a striking fact about the nature of these stars: while the internal structure is accurately determined by just Newtonian gravitation, their stability is entirely determined by GR.

To quantify this instability point, we can extremize the total energy of the star. As a function of the central density $\rho_{\text{c}}$, the total energy $E$ of this star is [30]

$$E(\rho_{\text{c}})=k_{1}MK\rho_{\text{c}}^{1/n}-k_{2}GM^{5/3}\rho_{\text{c}}^{1/3}-k_{4}G^{2}M^{7/3}\rho_{\text{c}}^{2/3},$$ (7)

where $k_{1}\approx 1.7558$, $k_{2}\approx 0.6390$, and $k_{4}\approx 0.9183$ are structure coefficients calculated for an $n=3$ polytrope (see Shapiro and Teukolsky [30] for derivation). The first two terms of Eq. (7) are the standard internal gas energy and Newtonian gravitational binding energy, respectively, and the final term represents the non-linear gravitational self coupling term arising from 1 Order GR corrections. A convenient dimensionless re-scaling follows from performing the transformations

$$\tilde{E} =\frac{G^{3/2}}{K^{\tilde{n}/2}}E,$$ (8a)

$$\tilde{M} =\frac{G^{3/2}}{K^{\tilde{n}/2}}M,$$ (8b)

$$x =GM^{2/3}\rho_{\text{c}}^{1/3},$$ (8c)

where here $\tilde{n}=1/(\Gamma_{1}-1)$ and we take the standard adiabatic index for a radiation pressure dominated plasma to be [25]:

$$\Gamma_{1}\approx\frac{4}{3}+\frac{\beta}{6}.$$ (9)

We assume $\Gamma_{1}$ is constant throughout the SMS since it is fully convective. Substituting the above into Eq. (7) gives

$$\tilde{E}=k_{1}\tilde{M}^{1/3-\beta/3}x^{1+\beta/2}-k_{2}\tilde{M}x-k_{4}\tilde{M}x^{2}.$$ (10)

The extremization follows on finding the first and second derivatives of Eq. (10) with respect to $x$, setting them both to zero, and solving simultaneously. The first derivative vanishing determines the point of hydrostatic equilibrium. The second derivative also vanishing determines where this hydrostatic configuration changes from stable to unstable. If we ignore the correction from GR and take the $\beta\rightarrow 0$ limit for the first derivative, we find the standard Newtonian rescaled mass for an $n=3$ polytrope:

$$\tilde{M}=\left(\frac{k_{1}}{k_{2}}\right)^{3/2}\approx 4.5547,$$ (11)

which is independent of $x$ (and hence the central density) as expected for an $n=3$ polytrope. In other words, a star with Newtonian gravitation and supported entirely by radiation (particles with relativistic kinematics) has zero total energy, and hence there is no energy cost for expanding or contracting this configuration. This self gravitating configuration is therefore very close to instability.

Such metastability means that even very small effects can have a dramatic effect on stability to push a configuration over the edge. As we have argued above, GR gives a negligible effect on the structure of the star, i.e., the run of pressure with density. However, the inclusion of a very simple 1-Order correction from GR, namely one that captures the non-linearity of gravitation through self-coupling, can make stable configurations in Newtonian gravitation unstable.

Including general relativistic corrections in this way and simultaneously solving the two extremization equations, leads to a critical value of $x$ where a hydrostatic configuration changes from stable to unstable:

$$x_{\text{crit.}}=\frac{k_{2}}{4k_{4}}\beta+\mathcal{O}(\beta^{2}).$$ (12)

This corresponds to a critical central density of

$$\rho_{\rm{c,\,crit.}}^{0} \approx\left(\frac{k_{2}}{4k_{4}}\right)^{3}\frac{\beta^{3}}{G^{3}M^{2}}$$

$$\approx 3.98\left(\frac{0.59}{\mu}\right)^{3}\left(\frac{{10}^{5}\,M_{\odot}}{M}\right)^{7/2}\,{\rm g}\,{\rm cm}^{-3}$$ (13)

where we have used Eq. (1) to relate $\beta$ to the star’s composition and mass. The physical interpretation of Eq. (13) has two parts. First, note that the rest mass of an SMS stems almost entirely from baryons, whereas the pressure is mostly from radiation. Increasing the mean molecular weight $\mu$ by fusing baryons in composite nuclei has the effect of lowering the gas pressure contribution to the total pressure. That, in turn, decreases the critical density for the onset of instability. Secondly, increasing the mass $M$ of the SMS requires a higher pressure for hydrostatic equilibrium, implying a higher temperature and a higher pressure contribution from radiation. Hence, increasing mass also serves to push the critical central density lower.

II Collisionless Dark Matter

We will now consider the effects of including a cloud of collisionless dark matter that permeates the entire star. Assuming standard spherical accretion, we expect the density profile of this dark matter cloud to be given by [30]

$$\rho_{\text{DM}}(r)=\rho_{\text{DM}}^{\infty}\left(1+\frac{2GM(r)}{v_{\infty}^{2}r}\right)^{1/2}$$ (14)

where $\rho_{\infty}$ and $v_{\infty}$ are the typical dark matter density and dark matter particle velocities far from any gravitational sources, respectively, and $M(r)$ is the total mass enclosed in a sphere of radius $r$.

The final term in the above equation represents the ratio of the gravitational potential of a dark matter particle to its specific kinetic energy. It permits two natural limits: $\tfrac{1}{2}v_{\infty}^{2}\gg GM(r)/r$ and $\tfrac{1}{2}v_{\infty}^{2}\ll GM(r)/r$ which, following Ref. [24]’s convention, we will call “hot” and “cold” dark matter, respectively. It should be noted that this nomenclature is not to be confused with the cosmological notion of “Hot Dark Matter” and “Cold Dark Matter,” which refers to the kinematics of dark matter at its decoupling epoch and how these kinematics affect early structure formation. Here, we will assume that both our “hot” and “cold” limits refer to dark matter particles with non-relativistic kinematics.

II.1 “Hot” Dark Matter

In the “hot” dark matter limit, Eq. (14) reduces to the simple relation

$$\rho_{\text{DM}}(r)=\rho_{\text{DM}}^{\infty}\equiv\rho_{\text{DM}},$$ (15)

i.e. a constant density fluid to first order. The energy added to the system by the dark matter is calculated by first determining the gravitational potential created by the fluid. From Poisson’s Equation, this dark matter has the harmonic potential

$$\Phi_{\text{HDM}}(r)=\frac{2\pi}{3}G\rho_{\text{DM}}r^{2}+C,$$ (16)

where we have assumed regularity at the origin and C is a constant of integration. To calculate the total gravitational energy contributed by this fluid, we use

$$E_{\text{HDM}} =\int\mathrm{d}m\,\Phi_{\text{HDM}}$$

$$=\frac{8\pi^{2}}{3}G\rho_{\text{DM}}\int_{0}^{R}\rho(r)r^{2}\,\mathrm{d}r+CM.$$ (17)

We will assume that $M$, $\mathrm{d}m$, and $\rho(r)$ are dominated by the baryon rest mass. Note that $\rho(r)$ very closely follows the profile of an $n=3$ polytrope. For such a polytrope in Newtonian gravitation, $M$ is independent of the central density, so we will ignore it for the purposes of extremization with respect to the central density $\rho_{\text{c}}$. The integral can be evaluated by employing the well known solutions to the Lane-Emden Equation of order 3, yielding

$$E_{\text{HDM}}=k_{\text{HDM}}MK\rho_{\text{DM}}\rho_{\text{c}}^{-2/3},$$ (18)

where $k_{\text{HDM}}=3.5845$ is the structure coefficient for the included dark matter and $K$ is the same from Eq. (2). Following the same transformation in Eqs. (8a)-(8c), the total energy now reads

$$\tilde{E}=k_{1}\tilde{M}^{1/3-\beta/3}x^{1+\beta/2}-k_{2}\tilde{M}x-k_{4}\tilde{M}x^{2}\\ +k_{\text{HDM}}\tilde{M}^{1/3}x_{\text{DM}}^{3}x^{-2}.$$ (19)

Extremization of this yields the following set of equations to be solved simultaneously after expanding in powers of the ratio of gas pressure to total pressure $\beta$:

$$0 =\frac{k_{1}}{\tilde{M}^{2/3}}\left(1-2\frac{k_{\text{HDM}}}{k_{1}}\frac{x_{\text{DM}}^{3}}{x_{\text{crit.}}^{3}}\right)-k_{2}-2k_{4}x_{\text{crit.}}+\mathcal{O}(\beta),$$ (20a)

$$0 =-2k_{4}+\frac{k_{1}}{2\tilde{M}^{2/3}}\frac{\beta}{x_{\text{crit.}}}\left(1+12\frac{k_{\text{HDM}}}{k_{1}\beta}\frac{x_{\text{DM}}^{3}}{x_{\text{crit.}}^{3}}\right)+\mathcal{O}(\beta^{2}).$$ (20b)

After identifying that the ratio $x_{\text{DM}}^{3}/x^{3}$ is related to the dark matter mass fraction $f$ via

$$\frac{x_{\text{DM}}^{3}}{x_{\text{crit.}}^{3}}=\frac{\rho_{\text{DM}}}{\rho_{\text{c,\,crit.}}}=\frac{\rho_{\text{DM}}}{C_{3}\rho_{\text{avg.}}}=\frac{M_{\text{DM}}}{C_{3}M}=\frac{f}{C_{3}},$$ (21)

where $C_{3}\approx 54.18$ is the core-to-average-density ratio (a.k.a “condensation parameter”) for a $n=3$ polytrope, we find that

$$\rho_{\text{c,\,crit.}}\approx\frac{1}{G^{3}M^{2}}\left(\frac{k_{2}}{4k_{4}}\beta+\frac{3k_{1}k_{\text{HDM}}}{k_{2}k_{4}}\frac{f}{C_{3}}\right)^{3},$$ (22)

which for typical values of $\mu$, $M$, and $f$ is

$$\rho_{\text{c,\,crit.}}^{\text{HDM}}\approx\rho_{\text{c,\,crit.}}^{0}\times\\ \left[1+0.197\left(\frac{\mu}{0.59}\right)\left(\frac{M}{10^{5}\,M_{\odot}}\right)^{1/2}\left(\frac{f}{0.01}\right)\right]^{3},$$ (23)

where $\rho_{\text{c,\,crit.}}^{0}$ is the critical central density without dark matter given in Eq. (13). Here, we have ignored $\mathcal{O}(\beta^{4})$ and $\mathcal{O}(f^{4})$ terms. It should be noted that all higher order terms in $f$ have positive coefficients, so Eq. (23) give a slight underestimate to the true increase in $\rho_{\text{c, crit.}}$ for larger $f$. Eq. (23) recovers the standard case when $f\rightarrow 0$, as expected.

In Figure 1, we show the critical central density $\rho_{\text{c,\,crit.}}$ at the onset of GR instability as a function of dark matter mass fraction $f$ in the hot dark matter case (dashed lines) for various SMS masses between $10^{4}$ and $10^{8}\,M_{\odot}$. While the relative increase is not too significant for smaller mass SMSs, the effect is quite strong for larger masses (e.g. ${\gtrsim}~{}10^{5}\,M_{\odot}$), increasing $\rho_{\text{c,\,crit.}}$ by as much as 4 orders of magnitude at $f\sim 20\%$.

Figure 1: Critical central density for the onset of general relativistic instability, $\rho_{\rm c,crit}$, for SMSs vs. total dark matter mass fraction $f$ for different SMS masses. Dashed lines are the “hot” dark matter limit and solid lines are the “cold” dark matter limit (Eqs. (23) and (32), respectively).

II.2 “Cold” Dark Matter

The alternative limit of Eq. (14) reduces to the profile

$$\rho_{\text{DM}}(r)=\frac{\rho_{\text{DM}}^{\infty}}{v_{\infty}}\sqrt{\frac{GM}{r}},$$ (24)

where we have also made the same approximation as Butler et al. [24] that since an $n=3$ polytrope is so centrally condensed, it is reasonable to treat the SMS as a point source of mass $M$. Without this assumption, the dark matter density profile would feature a flattened core as opposed to the above divergent cusp at the origin. Such a flat shape would be closer to the limit described in Section II.1, so the cuspy profile offers sufficiently disparate extreme from the “hot” dark matter case to allow us to explore the full range of dark matter effects we would expect to see.

The mass fraction of this dark matter is found by integrating the density throughout the entire star and dividing by the mass:

$$f=\frac{M_{\text{DM}}}{M} =4\pi\sqrt{\frac{2G}{M}}\frac{\rho_{\text{DM}}}{v_{\infty}}\int_{0}^{R}\mathrm{d}r\,r^{3/2},$$

$$=\frac{8\pi\sqrt{2}}{5}\sqrt{\frac{GR^{5}}{M}}\frac{\rho_{\text{DM}}}{v_{\infty}},$$

$$f =C_{\text{CDM}}\frac{x_{\text{DM}}^{3}}{x^{5/2}},$$ (25)

where $C_{\text{CDM}}\approx 60.0151$ is the “cold” dark matter equivalent to the condensation parameter in the previous section and is calculated in part from solutions to the Lane-Emden equation of order 3, $x$ is the reparameterized central baryonic density given by Eq. (8c), and $x_{\text{DM}}$ is defined as

$$x_{\text{DM}}^{3}=G^{3}M^{2}\frac{\rho_{\text{DM}}^{\infty}}{v_{\infty}}.$$ (26)

This definition is similar to $x_{\text{DM}}$ of the previous section with the inclusion of $v_{\infty}$, the typical DM particle velocity far from the SMS. Since the only the ratio $\rho_{\text{DM}}^{\infty}/v_{\infty}$ appears in our calculations, both have been collapsed into the single quantity $x_{\text{DM}}$.

Solving Poisson’s Equation again for the gravitational potential caused by this “cold” dark matter, we find that

$$\Phi_{\text{CDM}}(r)=\frac{16\pi\sqrt{2}}{15}G^{3/2}M^{1/2}\frac{\rho_{\text{DM}}^{\infty}}{v_{\infty}}r^{3/2},$$ (27)

where we again assume regularity at the origin and have ignored the constant of integration for the same reason as the previous section. The total gravitational energy is then calculated to be

$$E_{\text{CDM}} =\int\mathrm{d}m\,\Phi_{\text{CDM}}$$

$$=k_{\text{CDM}}MK^{3/2}\frac{\rho_{\text{DM}}^{\infty}}{v_{\infty}}\rho_{\text{c}}^{-1/2}$$ (28)

where $k_{\text{CDM}}\approx 14.066$ is calculated following a similar process to the previous section. The rescaled energy per Eqs. (8a)-(8c) and Eq. (26) is then

$$\tilde{E}_{\text{CDM}}=k_{\text{CDM}}x_{\text{DM}}^{3}x^{-3/2},$$ (29)

so the total rescaled energy is

$$\tilde{E}=k_{1}\tilde{M}^{1/3-\beta/3}x^{1+\beta/2}-k_{2}\tilde{M}x-k_{4}\tilde{M}x^{2}\\ +k_{\text{CDM}}x_{\text{DM}}^{3}x^{-3/2}.$$ (30)

Following the same process as the “hot” dark matter case with the mass fraction in Eq. (25), we find that

$$\rho_{\text{c,\,crit.}}\approx\frac{1}{G^{3}M^{2}}\left(\frac{k_{2}}{4k_{4}}\beta+\frac{15k_{2}k_{\text{CDM}}}{8\tilde{M}^{1/3}k_{1}k_{4}}\frac{f}{C_{\text{CDM}}}\right)^{3},$$ (31)

where $\tilde{M}=4.5547$ (per Eq. (11)) and we have again ignored 4 order terms in $\beta$ and $f$. Again, for typical values of $\mu$, $M$, and $f$, this reduces to

$$\rho_{\text{c,\,crit.}}^{\text{CDM}}\approx\rho_{\text{c,\,crit.}}^{0}\times\\ \left[1+0.262\left(\frac{\mu}{0.59}\right)\left(\frac{M}{10^{5}\,M_{\odot}}\right)^{1/2}\left(\frac{f}{0.01}\right)\right]^{3},$$ (32)

which has a slightly stronger dependence on the dark matter mass fraction than the “hot” dark matter case, as evidenced in Figure 1.

III DM Response to SMS Contraction

As the star quasistatically contracts, we would like to quantify how much the mass fraction of dark matter within the star changes. If the fraction of dark matter decreases for a hypothetical decrease in radius (enclosing the baryons), the SMS would be enclosing less total mass. That may help stabilize the configuration. Alternatively, or perhaps in addition to this effect, if SMS contraction serves to increase the velocity dispersion of the dark matter, the effective pressure provided by the dispersion would increase and also stabilize the configuration.

From Eqs. (21) and (25), the fractional change in dark matter fraction $f$ is

$$\frac{\delta f}{f}=\frac{\delta\rho_{\text{DM}}}{\rho_{\text{DM}}}-\kappa\frac{\delta\rho_{\text{c}}}{\rho_{\text{c}}},$$ (33)

where $\kappa=1$ or $5/6$ for the “hot” and “cold” limits, respectively.

The largest potential deviation in dark matter density from its background value would occur where $M(r)/r$ achieves its maximum value, which is at $r=R$, so we will consider that region and small changes in the SMS’s radius $\delta R<0$ where $\delta R/R\ll 1$. Since $\rho_{\text{c}}\sim M/R^{3}$,

$$\frac{\delta\rho_{\text{c}}}{\rho_{\text{c}}}=-3\frac{\delta R}{R}$$ (34)

for a fixed mass, and from Eq. (14) we find that

$$\frac{\delta\rho_{\text{DM}}}{\rho_{\text{DM}}}=-\frac{GM}{v_{\infty}^{2}R}\left(1+\frac{2GM}{v_{\infty}^{2}R}\right)^{-1}\frac{\delta R}{R}.$$ (35)

We will now consider both the “hot” and “cold” dark matter limits to calculate this fractional change.

III.1 “Hot” Dark Matter Case

To first order in $GM/(v_{\infty}^{2}R)$, Eq. (35) reduces to

$$\frac{\delta\rho_{\text{DM}}}{\rho_{\text{DM}}}\approx-\frac{GM}{v_{\infty}^{2}R}\frac{\delta R}{R},$$ (36)

so

$$\frac{\delta f}{f}\approx 3\left(1-\frac{GM}{3v_{\infty}^{2}R}\right)\frac{\delta R}{R}.$$ (37)

Therefore, the total mass fraction of dark matter would tend to decrease since $\delta R<0$, unless $GM/(6R)>\tfrac{1}{2}v_{\infty}^{2}$ which is the opposite of the limit we are considering, so this can be further reduced to

$$\left(\frac{\delta f}{f}\right)_{\text{HDM}}\approx 3\frac{\delta R}{R}.$$ (38)

III.2 “Cold” Dark Matter Case

If $GM/(v_{\infty}^{2}R)\gg 1$, then Eq. (35) becomes simply

$$\frac{\delta\rho_{DM}}{\rho_{DM}}\approx-\frac{1}{2}\frac{\delta R}{R},$$ (39)

so

$$\left(\frac{\delta f}{f}\right)_{\text{CDM}}\approx 2\frac{\delta R}{R},$$ (40)

which is again negative for $\delta R<0$, indicating that the total dark matter mass fraction also decreases in this limit but by a less severe amount compared to the “hot” dark matter case.

III.3 Physical Origin of Enhanced Stability

The stability of a purely standard model star (i.e., composed of baryons, electrons, photons, etc.) is determined by its pressure-averaged adiabatic index as given in Eq. (5). A heuristic physical picture for this effective “spring constant” for the star’s material is as follows: Squeeze the star radially and rapidly (so rapidly that heat does not flow and so entropy remains constant) to higher density and ask how much the pressure increases – i.e., how much does the material “push back” when quickly perturbed to higher density? Of course, when the pressure support for a self-gravitating configuration stems from particles with relativistic kinematics (e.g., photons in the SMS case), that configuration is close to instability. For example, this case would correspond to a zero total energy configuration in purely Newtonian gravitation, therefore pressure forces will grow in lock-step with the increased gravitational force from compressing the star. This means the configuration is in a meta-stable state as it costs no energy to expand or contract the star, making it ripe for instability once the non-linearity of GR comes into play.

When collisionless dark matter is added to this picture there are several effects to consider. We envision that the standard model portion of the star, contained inside a radius $R$, is embedded in an extended halo of dark matter. When the star is radially and adiabatically squeezed to higher baryon density there are two questions that arise: (1) What is the effect of reducing $R$ in the limit where the dark matter distribution in space is fixed so that less dark matter is enclosed inside the SMS radius?; and (2) If the standard model star is squeezed to higher density, what is the response of the dark matter spatial distribution and how does this affect stability?

To the first issue, as discussed previously, enclosing less total mass as the SMS is squeezed is a contributing factor in the enhanced stability afforded by dark matter. Indeed, as evidenced in Eqs. (37) and (40), the total dark matter mass fraction decreases in both dark matter kinematic limits. Curiously, this total mass decrement effect is less pronounced for the “cold” case despite it producing a stronger stabilizing effect (c.f. the coefficients in the correction terms in Eqs. (32) and (23) and Figure 1). This indicates that dark matter-induced stabilization is not only due to enclosing less total mass.

As for the second issue, the gravitational response of the dark matter spatial distribution to a perturbation in the standard model star is second order, so it is not accounted for in our calculations. The “cold” dark matter case possesses a more centrally condensed dark matter distribution, which changes the outcome of the extremization calculation as compared to the “hot” case. In the end, adding non-relativistic stress from dark matter to the supermassive star in general, and doing so in a way that provides a nearly fixed background, is inherently stabilizing.

IV Nuclear Burning and Neutrino Losses

It should be expected that, at least for some masses in the range ${\sim}10^{4}\text{--}10^{8}\,M_{\odot}$, nuclear (hydrogen) burning on the PPI chain is taking place at a high enough rate to significantly contribute to the total energy budget of the SMS. This rate depends on both the density and temperature throughout the star, both of which decrease as the mass increases. Specifically, at the onset of collapse where the density and temperature are highest, $\rho_{\text{c}}\sim M^{-7/2}$ and $T_{c}\sim M^{-1}$, the latter following from the definition of an $n=3$ polytrope. Since including dark matter can raise the critical central density (and hence the temperature) required before collapse, we can also investigate how the dark matter mass fraction impacts the energy released by the PPI chain throughout the star. Chen et al. [31] investigates Helium burning in SMSs. This could be relevant for SMSs light enough to have a stable Hydrogen burning phase.

While significant nuclear burning would potentially lead the star through a stable “main-sequence” phase of evolution, neutrino emission from the high entropy plasma could carry that extra energy generation away from the core, providing no direct pressure support and hence adding no more stability to the configuration. At such high entropies and temperatures, the neutrino luminosity may outshine that of nuclear burning completely, causing the star to contract more quickly and leading it closer to GR instability.

IV.1 PPI Chain Nuclear Burning

For zero metallicity, protons $p$ can only fuse to Helium nuclei $\alpha$ via the PPI chain $4p\rightarrow\alpha+2e^{+}+2\nu_{e}$, where $e^{+}$ indicates positrons and $\nu_{e}$ indicates electron neutrinos. In total, each such reaction releases about $26$ MeV of energy to the plasma while the neutrinos carry away about $0.5$ MeV away from the star as the likelihood they will transfer energy to the plasma is vanishingly small [32]. The energy generation rate for this reaction is calculated to be (in $\mathrm{erg\,cm^{-3}\,s^{-1}}$)

$$Q_{\text{PPI}}=c_{1}\rho^{2}X_{\text{H}}^{2}T_{6}^{-2/3}\exp\left(-c_{2}T_{6}^{-1/3}\right)\\ \times\left(1+c_{3}T_{6}^{1/3}+c_{4}T_{6}^{2/3}+c_{5}T_{6}\right),$$ (41)

where $\rho$ is the density in $\mathrm{g\,cm^{-3}}$, $X_{\text{H}}\approx 0.75$ is the total hydrogen mass fraction, $T_{6}=T/(10^{6}\,\mathrm{K})$ is the temperature in units of $10^{6}$ Kelvins, and ${c_{1}=2.319\times 10^{6}}$, $c_{2}=33.81$, $c_{3}=0.0123$, $c_{4}=0.0109$, $c_{5}=0.00095$ are numerical constants (see Clayton [32] for derivation).

To determine if the PPI chain provides a significant contribution to the energy budget of the star, we will compare the integrated luminosity at the surface, i.e.

$$L_{\text{PPI}}=\int_{0}^{R}Q_{\text{PPI}}4\pi r^{2}\mathrm{d}r$$ (42)

to the Eddington luminosity, which is approximately [30]

$$L_{\text{Edd.}}\approx 1.3\times 10^{38}\left(\frac{M}{M_{\odot}}\right)\,\mathrm{erg\,s^{-1}},$$ (43)

which comes from considering at what photon luminosity momentum transfer via photon-electron Thomson scattering balances the gravitational force.

IV.2 Neutrino Losses

The high entropy per baryon ($s\gtrsim 1000$) required for hydrostatic equilibrium in SMSs implies a significant $e^{\pm}$ pair content in electromagnetic equilibrium $e^{+}~{}+~{}e^{-}~{}\leftrightarrow~{}2\gamma$. These electrons and positrons can both annihilate and interact with the local radiation field to produce neutrinos via [26, 27, 33]:

$$e^{-}+e^{+} \rightarrow\nu+\bar{\nu}\,$$ (pair production)

$$\gamma+e^{\pm} \rightarrow e^{\pm}+\nu+\bar{\nu}\,$$ (photo-emission)

where $\gamma$ represents an incident photon and $\nu$ ($\bar{\nu}$) refers to (anti-)neutrinos of any flavor. These processes are expected to be the dominant source of neutrino emissivity in the low density and low-to-moderate temperature regimes in SMSs. There is also a process where plasmons ($\tilde{\gamma}$) decay to neutrino pairs ($\tilde{\gamma}\rightarrow\nu+\bar{\nu}$), but this only dominates in high density and temperature regimes.

If the electrons and positrons are non-relativistic and non-degenerate, the energy density loss rates for the two dominant neutrino production processes are calculated to be [32, 26] (in $\mathrm{erg\,cm^{-3}\,s^{-1}}$):

$$Q_{\text{pair}} \approx\left(4.9\times 10^{18}\right)T_{9}^{3}\exp\left(-\frac{11.86}{T_{9}}\right),$$ (44)

$$Q_{\text{photo.}} \approx\left(0.98\times 10^{8}\right)\frac{\rho}{\mu}\,T_{9}^{8},$$ (45)

where $\rho$ is the mass density in $\mathrm{g\,cm^{-3}}$, $\mu$ is the mean molecular weight, and $T_{9}=T/(10^{9}\,\mathrm{K)}$ is the temperature in units of $10^{9}$ Kelvins. While the pair production mechanism has a much higher pre-factor than the photo-emission process, the former suffers an exponential suppression for low temperatures that the latter channel does not, actually making that channel dominant in higher mass SMSs.

Using a similar process as in the previous section, we calculate the combined luminosity from these neutrino emission processes to be

$$L_{\mathrm{\nu,\,Tot.}}=\int_{0}^{R}\left(Q_{\mathrm{pair}}+Q_{\mathrm{photo.}}\right)4\pi r^{2}\,\mathrm{d}r.$$ (46)

While the PPI chain does emit two $\nu_{e}$ per reaction, the energy they carry away is far smaller than the energy released from fusing four protons into one $\alpha$ particle, so we will ignore that contribution from the total neutrino luminosity.

Figure 2: Symmetric log-log triptych plot showing total luminosity (burning minus neutrino losses) normalized to each SMSs’ Eddington luminosity versus the SMS’s central density as a fraction of its critical central density for different DM mass fractions. Solid (dashed) curves correspond to the “cold” (“hot”) kinematic limit. The top panel is for a star with no DM, the middle for a star with 1% DM by mass, and the bottom for a star with 10% DM by mass, all for various SMS masses between $10^{4}\text{--}10^{8}\,M_{\odot}$. Indicated in each panel is a dashed gray line at 1% Eddington luminosity to indicate when burning is “important” for the evolution of the star. Note that these should not be taken to be evolutionary plots since the mean molecular weight $\mu$ is fixed.

To examine how the PPI chain luminosity compares to the neutrino luminosity as the star contracts, we can simply scale down the central density and let it approach the critical central density from below. Assuming that the run of density and temperature closely follows an $n=3$ polytrope, this is a valid analysis as such a polytrope is entirely self similar (i.e. the mass is independent of the central density and only depends on $s$ and $\mu$ via the polytropic constant of proportionality $K$).

In Figure 2, we present this result in three symmetric log-log plots showing the total luminosity from both sources ($L_{\text{PPI}}~{}-~{}L_{\mathrm{\nu,\,Tot.}}$) as a fraction of the Eddington luminosity for various SMS masses and DM mass fractions. We subtract the neutrino luminosity since it carries away energy from the star, as opposed to creating energy like the PPI process. This plot should not be treated as an evolutionary plot since the curves shown in the figure are calculated with a fixed mean molecular weight, $\mu=0.59$, i.e., the value corresponding to a fully ionized primordial gas. Nuclear evolution, in the sense of burning via incorporating baryons into nuclei, should serve to increase the mean molecular weight $\mu$ as higher mass nuclei are created.

We would also expect that when nuclear burning is “important” (for the purposes of this paper chosen to be ${\sim}1\%$ of a star’s Eddington luminosity) the central density is held almost constant until its fuel is exhausted. This would prevent an unabated contraction of the star to the critical central density for GR instability, if only for a relatively short duration.

Since adding dark matter serves to increase the critical central density, more dark matter can push massive stars to higher central densities without causing them to become unstable. Both the PPI process and the relevant neutrino processes depend sensitively on the central temperature (and hence the central density), so a delicate balancing act ensues between energy generation from proton fusion and energy loss from neutrino emission.

As evidenced from the top plot (i.e. zero dark matter) in Figure 2, specifically for SMSs of masses ${\lesssim}10^{5}\,M_{\odot}$, increasing the central density increases the fraction of PPI energy generation, and hence net energy deposition, but only up to a point. Eventually, the neutrino losses from the increasing temperature win out, and the star is net losing energy instead of generating energy.

Should the quasi-statically contracting SMS encounter this situation during its evolution, the neutrino emission would cause the star to contract more quickly. This is because the neutrinos carry away the plasma’s entropy faster than Hydrogen burning can generate it throughout the configuration. Per the solutions to the Lane-Emden equations, an $n=3$ polytrope’s radius scales as $R\sim K^{1/2}\sim s^{2/3}$, and since $\rho_{\text{c}}\sim R^{-3}$, the central density will increase as the radius decreases as neutrinos remove entropy from the SMS. This serves to push the SMS closer to instability.

The situation changes slightly for the ${\sim}1\%$ dark matter mass fraction case (middle panel) and quite dramatically for the ${\sim}10\%$ case (bottom panel). At a given central density fraction $\rho_{\text{c}}/\rho_{\text{c,\,crit}}$, increasing DM mass fraction corresponds to a higher absolute density on each mass curve. If we examine the $10^{5}\,M_{\odot}$ curve, for example, this mechanism allows for the energy generation fraction from PPI to quite closely approach $1\%$ of the Eddington luminosity as the DM mass fraction increases. This potentially permits such stars to burn stably on the main sequence instead of a direct collapse as expected in the standard case. This would drastically change their evolution as opposed to the case without dark matter.

V Conclusion

These SMSs are quite interesting case studies, primarily because their existence provides a potential physical explanation for several observational mysteries that have only deepened as more advanced telescopes have come online. Not only has JWST, as an example, seen plenty of high mass quasars at high redshift which puts the standard picture of solar mass remnant accretion into doubt, it may also have observed SMSs supported in some way by dark matter prior to collapse in the JADES survey, per Freese et al. [34].

Considering that these stars are “trembling on the verge of instability” due to the Feynman-Chandrasekhar Instability, it begs the question as to how they would be observed at all prior to collapse. One possible avenue to support these stars is through some form of non-interacting and non-relativistic dark matter. The enhanced stability would certainly also have implications for the nuclear and neutrino evolution of the stars.

In Section I, we derived the standard result for the critical central density beyond which a hydrostatic SMS is unstable to collapse. This rested upon extremizing the total energy of the star with respect to the central density. We found that, as expected, higher mean molecular weight plasmas and higher mass SMSs are more unstable to collapse. Both of these factors serve to decrease the critical central density.

In Sections II.1 and II.2, we derived analytically how the critical central density changes as a function of the dark matter mass fraction within the star. We considered two limiting cases of the dark matter kinematics, the “hot” case and the “cold” case, which compared the typical DM particle’s specific kinetic energy far away from the SMS to the gravitational potential felt at the SMS surface, and found that in both cases the critical central density for the onset of general relativistic instability increased as the dark matter mass fraction increases. The “cold” case shows a slightly stronger effect, and the relative increase is greater for more massive stars in both kinematic limits.

In Section III we investigated a possible mechanism for the increased stability. We showed that if the star contracts, the total mass fraction of dark matter inside the star decreases, so the total mass enclosed by the star also decreases. This added stability is therefore completely agnostic to any particle physics specific to the dark matter, so long as any consideration of such physics keeps the DM non-relativistic and non-interacting with the baryonic matter. In fact, if the dark matter possesses some self-annihilation cross-section, like in the model discussed in Freese et al. [34], this may further bolster the stability by providing an extra source of energy and entropy.

Finally, in Section IV, we considered by what degree Hydrogen burning could impact the evolution of these SMSs and how dark matter influences that burning. We also analyzed how neutrino emission from the high-entropy plasma would take away energy, potentially faster than what is generated from fusion, as the star approaches the critical central density from below. We found that, considering both energy generation from fusion and energy losses from neutrino emission, dark matter could allow heavier stars than previously predicted to stably burn on the main sequence at some point in their lives, increasing the rough upper limit of ${\sim}30,000\,M_{\odot}$ to about $50,000\,M_{\odot}$, depending on when one considers nuclear burning to be “important” to the evolution of the star.

More investigation is needed to understand these extreme objects. If they are to explain the high mass SMBHs in the early universe by providing a supermassive seed black hole, there should be enough of them formed to account for all or most of them. This is complicated by a suspect formation history, since if the collapsing gas is cool enough to contain molecular Hydrogen instead of atomic Hydrogen, the enhanced cooling efficiency via rotational and vibrational modes could favor a fragmented formation into many stars instead of one SMS [35]. However, enhanced cooling from early molecular Hydrogen formation via Beyond the Standard Model (BSM) channels (e.g. in Biermann and Kusenko [36]) can lead to very early star formation. Whether SMS formation may accompany this early star formation is an open question.

We would expect that nuclear burning, i.e. fusing bare nucleons into larger composite nuclei, would lower the critical central density for the onset of collapse since the mean molecular weight increases. The gas pressure is provided primarily by non-relativistic baryons, which is proportional to the total number of particles, so nuclear fusion reduces this number and hence lowers the gas pressure if sufficient burning takes place. While sufficient burning would serve to slow the SMS’s contraction by injecting energy in the region surrounding the core, the critical central density is being reduced nonetheless, hastening the approach to instability if the chemical composition of the star is significantly changed.

Including dark matter further complicates this picture. We saw that sufficient dark matter will change the evolution of some SMS masses that would normally collapse via GR instability, making them also vulnerable to burning induced instability. As dark matter raises the critical central density, it also permits these stars to host greater central temperatures, making neutrino emission more efficient at carrying away entropy. This would cause the SMS to collapse on a different adiabat, specifically one at lower entropy than expected, potentially altering the prospects for post-collapse detection.

If such stars existed and underwent sufficient burning, they could be detected from their collapse directly if they then explode (instead of directly onto a black hole), and indirectly from enriching their surrounding medium with burning products both pre and post-collapse [37, 38]. The enhanced neutrino emission from dark matter would also alter the gravitational wave signature from an anisotropic neutrino pulse as the homologous core’s entropy would be reduced for a given mass (see e.g. Li et al. [12] for details on this neutrino driven gravitational radiation for standard SMSs).

The authors would also caution those who may be unfazed by the supposedly ‘modest’ amount of dark matter invoked in this study. Consider a $10^{6}\,M_{\odot}$ SMS with ${\sim}1\%$ dark matter by mass. If the dark matter is at a constant density within the star, such a DM mass fraction corresponds to a density of ${\sim}10^{20}\,\mathrm{GeV\,cm^{-3}}$, which would be an over-dense region of dark matter about 20 orders of magnitude more dense than the expected background at redshift $z\sim 10$. If these DM supported SMSs had existed, this huge accumulation and concentration of dark matter greatly constrains models for dark matter particles. Specifically, models where these particles cannot alter their phase space density and clump together through self and/or standard model interactions would clearly be disfavored in this regard. An example of a model where particles can alter their phase space density appropriately is discussed in Feng et al. [39], which in turn also provides a different mechanism for seeding SMBHs.

Moreover, the huge concentration of dark matter that we discuss here likely also constrains the manner and environments of SMS formation. Whether such environments can be found, and whether a component of the dark matter can be compatible with the concentrations discussed here, remain intriguing and open questions. In the end, such speculation may be warranted given our relative ignorance of dark matter and dark sector physics, and given the existence of SMBHs at very high redshift.

Here we have used a very simplistic model, a non-rotating and isolated SMS, to explore the physical effects of a large dark matter content. Nonetheless, our findings pose several open questions. First, with the inclusion of dark matter, what is the interplay of nuclear burning, hydrodynamics, energy transport, neutrino emission, etc., in the evolution of the SMS up to instability and the subsequent collapse to a black hole? Second, are there specific observational “handles” on the collapse of a SMS with or without dark matter? For example, is the gravitational radiation signature [12] different in these two cases, and to what extent would it depend on the amount of dark matter? Third, in this paper we have not considered the effects of dark matter on scenarios where a supermassive core undergoes rapid accretion, as in the hylotropic configurations proposed in Ref. [40]. The stabilizing effect of dark matter on such hylotropic configurations has been recently considered [41]. These scenarios may provide the most plausible route for the production of an SMBH from an SMS-like object.

Finally, perhaps the most significant issue is how an SMS would form with the huge quantities of dark matter considered in this work and in others that invoke dark matter effects in the production of SMBHs. However, there is urgency for consideration of these issues because of the existence of SMBHs at very early epochs in the history of the Universe, and because both the sources and properties of dark matter remain mysterious.