Dark Matter Annihilation and Pair-Instability Supernovae

The advent of gravitational wave astronomy has revolutionised the study of compact objects. With the publication of catalogues of transient (merger) events, the population statistics of astrophysical black holes can be studied for the first time. This information enables us to probe the processes governing the massive star progenitors of these black holes, allowing us to test stellar structure theory. For example, it is possible to measure the rates of nuclear reactions, most importantly the ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction Farmer et al. (2019, 2020); Mehta et al. (2022). In addition, new physics beyond the Standard Model (BSM) can alter the structure, evolution, and formation of these objects either via their effects on the evolution of the stars through the pulsational pair-instability supernovae (PPISN) and pair-instability supernovae (PISN) phases Croon et al. (2020a, b); Straight et al. (2020); Sakstein et al. (2020); Ziegler and Freese (2021); Sakstein et al. (2022); Fernandez et al. (2022); Ziegler and Freese (2022), or via their effects on the formation of black hole binaries Croon and Sakstein (2023). Examples of the new physics studied include light ($m_{\rm DM}\lesssim 10$ keV) weakly interacting particles, intermediate mass particles ($m_{\rm DM}\lesssim 1$ GeV) strongly coupled to the Standard Model (SM), and dark energy. Continuing this program, in this work we study the effects of dark matter (DM) with $m_{\rm DM}\sim\mathcal{O}(1)$ GeV. In this mass range, the DM is too heavy to be thermally produced by the stellar plasma, but can be gravitationally captured if the ambient dark matter density is sufficiently large. The subsequent annihilation results in an injection of energy into the stellar material.

The results of our study can bu summarized as follows:

1. 1.

   When DM is heavier than $\sim 1$ GeV it is highly concentrated in the star’s core. The energy injection from annihilation augments the lifetime of core helium burning because less energy is required to maintain hydrostatic equilibrium. This gives more time for the ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction to operate, resulting in a larger ¹⁶O abundance at core helium depletion. This drives a more violent explosion, which causes the star to lose more mass.  This results in a reduction in the location of the lower edge of the upper black hole mass gap.

2. 2.

   When DM is lighter than $\sim 1$ GeV it is more diffuse, injecting energy thought the entire star. Energy injected into non-nunclear burning regions post-helium depletion drives a phase of evolution where the star is partly supported by DM annihilation. During this phase, the star’s central density and temperature decrease, allowing it to evade the pair-instability in the $\rho$–$T$ plane. Stars that would have otherwise undergone PPISN instead directly core collapse, forming heavier black holes than predicted in the absence of DM annihilation. This effect vanishes in heavier stars which require more energy for support than is be provided by DM annihilation, and these objects experience PPISN and PISN, resulting in a mass gap.

We exemplify the scenarios above by studying two DM masses: $m_{\rm DM}=1$ GeV and $m_{\rm DM}=0.2$ GeV. Our results differ from those obtained in earlier work on dark heat injections and PPISN Ziegler and Freese (2021, 2022). We discuss this discrepancy below.

This paper is organized as follows. In section II we describe the DM model we adopt. In section III we describe our implementation of DM annihilation into our stellar structure code. Our results are presented in section IV. We discuss our results and conclude in section V.

To model the energy injections due to DM annihilation, we must first consider the appropriate dark matter number density profile in the star. This is a complicated procedure that requires a careful consideration of the effects of DM capture, annihilation, and evaporation — all of which depend on the size and structure of the object under consideration, the timescale for this object’s evolution, and the DM mass. Since our study is the first of its kind, we will make a number of simplifying assumptions that enable us to investigate the effects of DM annihilation across a wide range of masses and over the entire evolution of the star. One goal of this work is to determine whether a more comprehensive study is required, and over which range of DM masses. These approximations facilitate this.

First, we assume the dark matter is thermalised in the core. The thermalisation time can be estimated using the number of scatterings needed (given by the DM-SM mass ratio) times the DM mean free path ($\lambda=(\sigma_{0}n_{\rm SM})^{-1}$, where $\sigma_{0}$ is the cross section for DM-SM scattering) divided by the DM dispersion velocity in the star Iocco et al. (2008). Taking the latter to be the escape velocity at the surface (or the surface of the core), one obtains

where the last equality approximates the star as having uniform density. We note that the latter approximation is for illustrative purposes only. We expect that the result — that thermalisation happens on much smaller timescales than the typical timescales of stellar evolution — continues to hold true for more realistic density profiles.

Our second approximation is that annihilation equilibrium has been established on timescales larger than $\tau_{\rm eq}\sim\left(\Gamma_{\rm ann}\Gamma_{\rm cap}\right)^{-1/2}$. The capture rate we can estimate as $\Gamma_{\rm cap}=\Phi\pi R^{2}$ where $\Phi$ is the flux density of DM which is captured, given by Leane and Smirnov (2022)

where $f_{\rm cap}$ denotes the fraction of DM which is captured as it streams through the star, and where $\rho_{\rm DM}$ is the DM halo energy density. We can estimate the DM annihilation rate per effective volume as $\Gamma_{\rm ann}=\langle\sigma_{\rm ann}v\rangle/V_{\rm eff}$, with $1/V_{\rm eff}=\int_{V}n_{\chi}^{2}/\int_{V}n_{\chi}$ as in Croon and Smirnov (2023). This gives

where we have assumed an ambient DM density of $10^{9}\,\rm GeV\,cm^{-3}$, motivated by the results we present below that show deviations from the SM emerging for these values. Using Asteria Leane and Smirnov (2023), we find that the reference capture fraction $f_{\rm cap}=10^{-5}$ is relevant for DM-nucleon cross sections of $\sigma_{0}\lesssim 10^{-42}\,\rm cm^{2}$. For $\sigma_{0}=10^{-35}\,\rm cm^{2}$, $f_{\rm cap}$ approaches unity. To maximize the potential signal, we take $f_{\rm cap}=1$ below, implying that our results hold for $\sigma_{0}\gtrsim 10^{-35}\,\rm cm^{2}$.

Our third approximation, implicit in the approximation that the DM is in capture-annihilation equilibrium, is that evaporation is negligible. This is approximately valid for DM masses for which the evaporation rate is below the capture rate. We estimate the evaporation rate as $\Gamma_{\rm evap}=4\pi\Phi_{\rm J}(T)R_{\star}^{2}$, where $\Phi_{\rm J}(T)$ is the Jeans flux

with $n$ the DM number density in the star, $v=\sqrt{2T/m_{\rm DM}}$, and $v_{\rm esc}=\sqrt{2GM/R}$ being the surface escape velocity. This leads to evaporation masses on the order of $m_{\rm evap}\sim\mathcal{O}(0.1)$ GeV for $65{\rm M}_{\odot}$ stars, with a mild dependence on the cross section in the local thermal equilibrium regime (see below). Note that the evaporation mass can be lower in the presence of long-range forces in the dark sector Acevedo et al. (2023). The inclusion of such forces would lead to a slight modification of the profiles below.

Our fourth assumption is that the DM diffusion through the star is fast compared to both the timescales on which the total DM particle number significantly changes, and the stellar evolution. The diffusive timescale can be estimated as $t_{\rm diff}\sim(\Delta x)^{2}n_{\rm SM}\sigma_{0}/v_{T}\sim\rm\mathcal{O}(few)\,min$  Leane and Smirnov (2022). The advective timescale for these stars is typically much larger. Turning to the stellar evolution timescales, the relevant timescales in our simulation are small compared to (1) and (3) with the exception of the pulsations which can occur on timescales $\mathcal{O}(\rm days)$. If annihilation equilibrium would be assumed at this stage the energy injection would grow with the increase in radius of the star, resulting in a PISN whenever a PPISN is triggered. To conservatively model the behaviour of DM for smaller cross sections, we freeze its profile at the onset of the first pulsation.

The spatial profile that the DM assumes in the star depends on the Knudsen number. In the limit of small Knudsen numbers — that is, dark matter scatterings are frequent; the characteristic mean free path $\lambda$ is short compared to the core radius $\sim 0.5R_{\odot}$ as well as compared to the scale of temperature variations $|\nabla\ln T|^{-1}$ — the Local Thermal Equilibrium (LTE) assumption holds. We find that a $72{\rm M}_{\odot}$ star has $|\nabla\ln T|^{-1}>\mathcal{O}(10^{-1})\rm R_{\odot}$ in the core on the ZAHB and at the onset of the first pulsation. This means that $\sigma_{0}n_{\rm SM}\geq 10^{-11}\rm cm^{-1}$ defines the short mean free path regime. For a number density of $10^{25}\,\rm cm^{-3}$, this corresponds to $\sigma_{0}=10^{-36}\,\rm cm^{2}$. In the LTE regime, the DM profile $n_{\rm DM}$ in the star is given by the Maxwell-Boltzmann distribution Gould and Raffelt (1990):

Here $g(r)$ is the gravitational acceleration and $\alpha(r)$ is a “separation constant”, related to the diffusion coefficient in Leane and Smirnov (2022) by $\kappa=\alpha-5/2$. Equation (5) can be simplified in the limit that $\alpha$ is independent of $r$ — in this limit there also is a simple analytic expression which matches the numerical result closely Leane and Smirnov (2022), such that

In the opposite limit of large Knudsen number, the profile is expected to take on the isothermal distribution derived by Spergel and Press (1985). We leave an analysis of that regime for future work. The energy injection (per unit mass and time) is then given by

where we substitute in the profile (6), $\rho(r)$ is the density of the star at radius $r$, $\langle\sigma v\rangle$ is the annihilation cross section, and the dimensionless factor $f_{\nu}$ accounts for annihilation biproducts that free-stream out of the star. As we are interested in the maximal effect, we will take $f_{\nu}=1$. The DM annihilation luminosity in annihilation equilibrium is given by

from which it is seen that the energy injection must be independent of $\langle\sigma v\rangle$ in annihilation equilibrium Lopes and Lopes (2021). We can now normalise our energy injection profile (7) using (8).

We modified the stellar structure code MESA version 12778 Paxton et al. (2011, 2013, 2015, 2018) to include dark matter injection with the profile (6). This was accomplished by using the other_energy_implict hook. The full details of our implementation can be found in a reproduction package that accompanies this paper: https://zenodo.org/records/10056048. This includes MESA inlists and modifications. We briefly summarize the salient features of our stellar modeling choices here for completeness, referring the reader to Sakstein and Croon (2023) for the full details. Our simulations include convection described by mixing length theory, semiconvection, and overshooting; all with parameters given in the inlists in the reproduction package Sakstein and Croon (2023). We use the mass-loss scheme of Brott et al. (2011), which is a good description of wind-loss for massive stars. Our treatment of the pulses follows the prescription of references Farmer et al. (2019); Croon et al. (2020b), to which we refer the reader for the details. All nuclear reaction rates are set to the MESA defaults described in Farmer et al. (2019); Croon et al. (2020b) with the exception of the ${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ rate. This is known to be the largest source of uncertainty in the location of the black hole mass gap and the shape of the black hole mass spectrum, so we use the state-of-the-art rate table presented in Mehta et al. (2022)¹¹1We are grateful to J. de Boer for providing this table.. Our reproduction package contains a copy of this table Sakstein and Croon (2023).

Our implementation of the DM annihilation is as follows. We write equation (7) (having substituted equation (6)) as

with

With these definitions, equation (8) can be rearranged to find $\varepsilon_{0}$

where we took $f_{\nu}=1$, $\Gamma_{\rm cap}=\pi R^{2}\Phi$, and $\Phi$ is given by equation (2). At each time step, we calculate $\chi$ by integrating over the stellar profile and use equation (12) to calculate $\varepsilon_{0}$. We then use this to calculate the injection in each cell using (9).

In this work, we have investigated the effect of effect of dark matter annihilation on pair-instability supernovae. We modeled the DM injection profile using the local thermal equilibrium assumption (5). We identified two separate phenomena depending on the DM mass. If DM is heavier than $\sim\mathcal{O}(1)$ GeV, it is highly concentrated in the core. The energy injection from its annihilation prolongs the lifetime of core helium burning, reducing the C/O ratio leading to more violent explosions that result in lower mass black holes being formed, and a reduction in the location of the lower edge of the upper black hole mass gap. If DM is lighter than $\sim\mathcal{O}(1)$ GeV then it is diffused throughout the star and its energy injection from annihilation leads to a phase of partial DM support that resists the contraction that usually follows helium depletion. Instead, the star expands, which reduces the core temperature and density. Some stars that would have experienced the pair-instability evade it, leading to heavier black holes being formed and an increase in the location of the lower edge of the upper black hole mass gap. We found that each of these effects — more violent explosions for DM masses heavier than $\sim\mathcal{O}(1)$ GeV and evading the pair-instability for DM masses lower than $\sim\mathcal{O}(1)$ GeV — are exhibited when the ambient DM density $\rho_{\rm DM}\gtrsim 10^{9}\,\rm GeV\,cm^{-3}$. Such densities are only potentially realised in the very centers of DM spikes around galactic nuclei, and giant stars evolving in such backgrounds are expected to be exceedingly rare. We therefore anticipate that the effects of DM will be difficult to observe in current GW catalogs. In practice, the low metallicity ($Z=10^{-5}$) stars we have simulated are unlikely to inhabit galactic centers, but we expect that DM will effect higher metallicity stars in a qualitatively similar manner. The primary effect of metallicity is to alter the rate of mass loss to stellar winds during core helium burning ($\dot{M}\propto Z^{0.85}$). Since DM injections do not alter stellar winds, we do not expect our conclusions to change for higher metallicity objects.

The effects of an additional non-nuclear energy injection on massive stars were studied recently by Ziegler and Freese (2021, 2022). In those works, the authors injected a constant amount of energy per unit mass and time at each point in the star as a proxy for DM annihilation without specifying a profile. The works conclude that the upper black hole mass gap is filled in because stars can evade the instability. In contrast, in this work we find a mass gap for all DM masses and ambient DM densities that we simulated. We attempted to reproduce the results in these works based on the descriptions of the code given in the papers, but all of our models failed to converge (unfortunately, the code used by Ziegler and Freese (2021, 2022) is not publicly-available). The scenario with a uniform heat injection throughout the star implies that the DM profile is diffuse, and this scenario is therefore most similar to our $m_{\rm DM}=0.2$ GeV models. Additionally, if one attempts to do this comparison, the uniform heat injection studied in Ziegler and Freese (2021, 2022) would correspond to ambient DM densities $\rho_{\rm DM}\gg 3\times 10^{9}\,\rm GeV\,cm^{-3}$, in which case more stars would evade the pair-instability but, following our results, one would still expect sufficiently heavy stars to traverse the instability region and experience PPISN/PISN, leading to a gap. An exploration of a larger range of initial masses using their code could test this hypothesis.

As noted above, ambient DM densities $\rho_{\rm DM}>10^{9}{\rm GeV\,cm^{-3}}$ could only be realised deep inside hypothetical DM spikes around SMBHs and, consequentially, one would expect the mass gap-filling objects formed in this way to be extremely rare. In addition, the DM captured by a star is predicted to evaporate when $m_{\rm DM}\lesssim 0.1$ GeV, so a near-uniform injection scenario would require a DM model with an additional attractive force between DM particles with Compton wavelength $\lambda_{C}\sim R_{\odot}$ in order to be viable. Considering our results, we expect that mass-gap filling DM models represent only a portion of the more general space of DM models, namely those with light DM masses, likely requiring an attractive dark force that prevents evaporation, and extremely high ambient DM densities.

Our study constitutes a first systematic exploration of the effect of DM heat injections with a LTE profile on the post-main sequence evolution of heavy stars. DM may have an effect on cooler objects at lower ambient density. We plan to explore this further in future work.

MESA version 12778, MESASDK version 20200325, Asteria version 1, Mathematica version 12.

We thank Rebecca Leane, Juri Smirnov, and Aaron Vincent for useful discussions. DC is supported by the STFC under Grant No. ST/T001011/1. This material is based upon work supported by the National Science Foundation under Grant No. 2207880. Our simulations were run on the University of Hawai ‘ i’s high-performance supercomputer MANA. The technical support and advanced computing resources from University of Hawai ‘ i Information Technology Services – Cyberinfrastructure, funded in part by the National Science Foundation MRI award #1920304, are gratefully acknowledged.