# Dark Matter Annihilation and Pair-Instability Supernovae

Djuna Croon<EMAIL_ADDRESS>Institute for Particle Physics
Phenomenology, Department of Physics, Durham University, Durham DH1 3LE, U.K.
Jeremy Sakstein<EMAIL_ADDRESS>Department of Physics & Astronomy,
University of Hawai ‘ i, Watanabe Hall, 2505 Correa Road, Honolulu, HI, 96822,
USA

###### Abstract

We study the evolution of heavy stars ($M\geq 40{\rm M}_{\odot}$) undergoing
pair-instability in the presence of annihilating dark matter. Focusing on the
scenario where the dark matter is in capture-annihilation equilibrium, we
model the profile of energy injections in the local thermal equilibrium
approximation. We find that significant changes to masses of astrophysical
black holes formed by (pulsational) pair-instability supernovae can occur when
the ambient dark matter density $\rho_{\rm DM}\gtrsim
10^{9}\rm\,GeV\,cm^{-3}$. There are two distinct outcomes, depending on the
dark matter mass. For masses $m_{\rm DM}\gtrsim 1$ GeV the DM is primarily
confined to the core. The annihilation increases the lifetime of core helium
burning, resulting in more oxygen being formed, fueling a more violent
explosion during the pair-instability-induced contraction. This drives
stronger pulsations, leading to lighter black holes being formed than
predicted by the standard model. For masses $m_{\rm DM}\lesssim 0.5$ GeV there
is significant dark matter in the envelope, leading to a phase where the star
is supported by the energy from the annihilation. This reduces the core
temperature and density, allowing the star to evade the pair-instability
allowing heavier black holes to be formed. We find a mass gap for all models
studied.

††preprint: IPPP/23/62

##  I Introduction

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 16O
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.

##  II Dark matter energy injections

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

$\begin{split}\tau_{\rm th}(r)=&\frac{m_{\rm DM}}{\rho_{\rm
SM}(r)\sigma_{0}}\sqrt{\frac{r}{2GM}}\\\ =&1.9\times 10^{2}{\rm
yr}\left(\frac{m_{\rm DM}}{1\rm GeV}\right)\left(\frac{\sigma_{0}}{10^{-44}\rm
cm^{2}}\right)^{-1}\\\
&\times\left(\frac{R_{\star}}{0.5R_{\odot}}\right)^{7/2}\left(\frac{M_{\star}}{10M_{\odot}}\right)^{-3/2}.\end{split}$
(1)

where the last equality approximates the star as having uniform density, for
reference.

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)

$\Phi=v_{\rm DM}\sqrt{\frac{8}{3\pi}}\left[1+\frac{3}{2}\left(\frac{v_{\rm
esc}}{v_{\rm DM}}\right)^{2}\right]\frac{\rho_{\rm DM}f_{\rm cap}}{m_{\rm
DM}}$ (2)

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

$\begin{split}\tau_{\rm eq}\sim&\left(\Gamma_{\rm ann}\Gamma_{\rm
cap}\right)^{-1/2}\\\ \sim&3\times 10^{4}{\rm yr}\sqrt{\frac{10^{-5}}{f_{\rm
cap}}\frac{m_{\chi}\,10^{-30}\text{ cm}^{3}/\text{s}}{\langle\sigma_{\rm
ann}v\rangle\text{ GeV}}}\end{split}$ (3)

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

$\Phi_{\rm J}(T)=\frac{nv}{2\sqrt{\pi}}\left(1+\frac{v_{\rm
esc}^{2}}{v^{2}}\right)\exp\left(-\frac{v_{\rm esc}^{2}}{v^{2}}\right)$ (4)

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 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):

$\left(\frac{n_{\rm DM}(r)}{n_{\rm
DM}(0)}\right)=\left(\frac{T(r)}{T(0)}\right)^{3/2}e^{-\int_{0}^{r}d\widetilde{r}\frac{\alpha(\widetilde{r})dT/d\widetilde{(}r)+m_{\rm
DM}g(\widetilde{r})}{T(\widetilde{r})}}$ (5)

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

$\left(\frac{n_{\rm DM}(r)}{n_{\rm
DM}(0)}\right)=\left(\frac{T(r)}{T(0)}\right)^{-1+\frac{1}{2}\left(1+\frac{m_{\rm
DM}}{m_{\rm SM}}\right)^{-3/2}}e^{-\int_{0}^{r}d\widetilde{r}\frac{m_{\rm
DM}g(\widetilde{r})}{T(\widetilde{r})}}.$ (6)

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

$\begin{split}\epsilon_{\rm DM}=f_{\nu}\frac{\langle\sigma v\rangle n_{\rm
DM}^{2}(r)m_{\rm DM}}{\rho(r)}\end{split}$ (7)

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

$L_{\rm DM}=4\pi\int_{0}^{R_{\star}}\rho(r)\epsilon_{\rm DM}(r)r^{2}dr=m_{\rm
DM}\Gamma_{\rm cap}$ (8)

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).

Figure 1: Results of DM injections described in text for two different DM
masses.

##  III Stellar Structure Code

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)111We 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

$\varepsilon(r)=\varepsilon_{0}\frac{f(r)^{2}}{\rho(r)}$ (9)

with

$\displaystyle\varepsilon_{0}$ $\displaystyle=\langle\sigma v\rangle m_{\rm
DM}n(0)^{2};\quad\textrm{and}$ (10) $\displaystyle f(r)$
$\displaystyle=\left(\frac{T(r)}{T(0)}\right)^{-1+\frac{1}{2}\left(1+\frac{m_{\rm
DM}}{m_{\rm SM}}\right)^{-3/2}}e^{-\int_{0}^{r}d\widetilde{r}\frac{m_{\rm
DM}g(\widetilde{r})}{T(\widetilde{r})}}.$ (11)

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

$\varepsilon_{0}=\frac{R^{2}\Phi m_{\rm
DM}}{4\chi};\quad\chi=\int_{0}^{R_{*}}r^{2}f(r)^{2}\mathrm{d}r,$ (12)

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).

##  IV Results

###  IV.1 Theoretical Considerations

Before proceeding to the numerical results, one can gain some insight into the
effects of DM annihilation from purely theoretical considerations. The initial
question of interest is: can the injection alter the location of the tracks in
the $\rho$–$T$ plane? The answer reveals whether the star will encounter the
pair-instability. The shape of the stellar tracks in the $\rho$–$T$ plane can
be found by assuming that the star is radiation-supported with equation of
state (EOS) $P\propto T^{4}$, which is reasonable for massive objects. In this
case, the equation for hydrostatic equilibrium gives Straight  _et al._ (2020)

$\log(\rho)=\frac{1}{3}\log(T)+\frac{1}{6}\log(M)+c,$ (13)

where $c$ is a constant that depends on Newton’s constant and the radiation
constant. A similar relation but with a different slope can be derived for the
case of gas pressure domination where $P\propto\rho T$, with the true EOS
being a combination of the two. Thus we see that the slope of the tracks
depend only on the equation of state, and therefore the DM injection only
alters the track if it is sufficiently strong that the EOS is significantly
modified. In terms of the DM mass, the profile (6) implies that for light DM
the injection is diffuse throughout the star and can inject energy into
regions where there is no nuclear burning and therefore dominate the EOS. For
heavier DM, the injection is concentrated in the core, and we observe that the
stellar tracks are unmodified before the onset of pair instability.

###  IV.2 Numerical Results

We simulated a grid of stars with masses in the range $40{\rm M}_{\odot}\leq
M\leq 90{\rm M}_{\odot}$ from the onset of helium burning through the regime
of pair-instability to their ultimate core collapse or PISN, depending on the
mass. We took the metallicity to be $Z=10^{-5}$ corresponding to population-
III stars. Two different DM masses were investigated, $m_{\rm DM}=1$GeV and
$m_{\rm DM}=0.2$GeV for a range of ambient DM densities $\rho_{\rm DM}$. We
fixed the DM velocity to be $v_{\rm DM}=220$ km/s corresponding to the
circular velocity of the Solar System through the Milky Way. Our reproduction
package includes the option of varying this Sakstein and Croon (2023), but we
chose not to do so since this is nearly degenerate with $\rho_{\rm DM}$. The
degeneracy is only lifted by the $\mathcal{O}(1)$ gravitational focusing
correction to the flux in equation (2). In addition, the DM density is
expected to vary over many orders-of-magnitude, whereas we expect only an
$\mathcal{O}(1)$ variation of $50\textrm{km/s}\lesssim v_{\rm DM}\lesssim
300\textrm{km/s}$ given the distribution of circular velocities of galaxies.

We varied $0.42\textrm{GeV cm}^{-3}\leq\rho_{\rm DM}\leq 5\times
10^{9}\textrm{GeV cm}^{-3}$ i.e., spanning DM densities from that of the solar
neighborhood to the largest densities one would expect for stars near the
galactic center in scenarios with a DM spike Gondolo and Silk (1999) (see Shen
_et al._ (2023) for recent observational constraints on such a spike in the
Milky Way based on the trajectories of the nearest stars). Such a large
exploration was necessary because the DM energy injection is a complicated
function of the stellar structure, temperature, and density, and, similarly,
the energy released from nuclear burning is a strong function of temperature
and density, making a theoretical comparison difficult. We found deviations
from the SM for $\rho_{\rm DM}\gtrsim 10^{9}$ GeV/cm3. Our results are shown
in figure 1, where we plot the final black hole mass as a function of zero-age
helium burning (ZAHB) mass for the SM and both DM cases we investigated. The
two masses studied show opposite effects. We elucidate each of these in turn
below.

Figure 2:  Core elemental abundances for a $59{\rm M}_{\odot}$ star. The 16O
abundance at helium depletion is indicated by the horizontal lines. As
explained in text, the DM heat injection slows down the evolution and leads to
a higher core oxygen abundance, exacerbating pair instability.

The inclusion of annihilating DM with mass $m_{\rm DM}=1$GeV results in
stronger PPISN, and a reduction in the ZAHB mass at which the PPISN/PISN
transition occurs. This leads to lower black hole masses, and a reduction in
the location of the lower edge of the upper black hole mass gap. DM with this
mass is concentrated in the core. The energy released from its annihilation
prolongs the lifetime of core helium burning because less nuclear fuel needs
to be burned to stave off gravitational collapse. During this time, the
${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ reaction reprocesses some of the
carbon into oxygen. The increased lifetime provides more time for this
process, ultimately enhancing the amount of oxygen produced at core helium
depletion. Both of these effects are exemplified in 2 for a $59{\rm
M}_{\odot}$ star. The increase in core oxygen levels drives in a more violent
explosion. This results in enhanced mass loss for stars that undergo PPISN, as
seen in Fig. 1.

The case where the DM mass is $m_{\rm DM}=0.2$GeV corresponds to a more
diffuse DM profile where there is a significant amount of energy being
injected in the star’s outer layers. As the temperature and density in these
layers are too low for significant nuclear burning, DM heat can overpower SM
heat at certain stages of the evolution. After core helium depletion, this
injection supports the star, preventing the contraction that would lead to
C-burning in the SM. The injected energy causes the star to expand, reducing
the core temperature and density, moving the location of the track in the
$\rho$–$T$ plane. Examples are shown in figure 3. The new track is akin to
that of a lower mass star. For sufficiently light objects, the pair-
instability is avoided, and hence no PPISN/PISN occurs. Nuclear burning
eventually proceeds through the usual channels, and the star ultimately core
collapses to a heavier black hole than predicted by the SM. Heavier objects
cannot evade the instability, and experience PPISN/PISN depending on the mass.
Therefore, a black hole mass gap still forms, but at higher masses than
predicted by the SM.

Figure 3: Example trajectory in core temperature and density of a $59{\rm
M}_{\odot}$ star in the SM and in the presence of DM energy injections, and
the corresponding evolution of the radius with model number in the simulation.
The rainbow shading of the DM graph indicates the model number, for ease of
comparison of the two panels. As expected, the trajectory is unaltered before
core helium depletion, indicated by “HD” (see text for details). After HD, the
injected energy can support the star against collapse, altering the trajectory
in such a way that the star avoids the pulsational pair instability
encountered in the SM.

##  V Discussion

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 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 significant deviations from the SM 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.

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.

## Software

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

## Acknowledgements

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.

## References

  * Farmer  _et al._ (2019) R. Farmer, M. Renzo, S. de Mink, P. Marchant, and S. Justham,  (2019), 10.3847/1538-4357/ab518b, arXiv:1910.12874 [astro-ph.SR] .
  * Farmer  _et al._ (2020) R. Farmer, M. Renzo, S. de Mink, M. Fishbach, and S. Justham, (2020), arXiv:2006.06678 [astro-ph.HE] .
  * Mehta  _et al._ (2022) A. K. Mehta, A. Buonanno, J. Gair, M. C. Miller, E. Farag, R. J. deBoer, M. Wiescher, and F. X. Timmes, Astrophys. J. 924, 39 (2022), arXiv:2105.06366 [gr-qc] .
  * Croon  _et al._ (2020a) D. Croon, S. D. McDermott, and J. Sakstein, (2020a), arXiv:2007.00650 [hep-ph] .
  * Croon  _et al._ (2020b) D. Croon, S. D. McDermott, and J. Sakstein, (2020b), arXiv:2007.07889 [gr-qc] .
  * Straight  _et al._ (2020) M. C. Straight, J. Sakstein, and E. J. Baxter, Phys. Rev. D 102, 124018 (2020), arXiv:2009.10716 [gr-qc] .
  * Sakstein  _et al._ (2020) J. Sakstein, D. Croon, S. D. McDermott, M. C. Straight, and E. J. Baxter, Phys. Rev. Lett. 125, 261105 (2020), arXiv:2009.01213 [gr-qc] .
  * Ziegler and Freese (2021) J. Ziegler and K. Freese, Phys. Rev. D 104, 043015 (2021), arXiv:2010.00254.
  * Sakstein  _et al._ (2022) J. Sakstein, D. Croon, and S. D. McDermott, Phys. Rev. D 105, 095038 (2022), arXiv:2203.06160 [hep-ph] .
  * Fernandez  _et al._ (2022) N. Fernandez, A. Ghalsasi, S. Profumo, N. Smyth, and L. Santos-Olmsted, (2022), arXiv:2208.08557.
  * Ziegler and Freese (2022) J. Ziegler and K. Freese, (2022), arXiv:2212.13903.
  * Croon and Sakstein (2023) D. Croon and J. Sakstein, Phys. Rev. D 108, 015034 (2023), arXiv:2208.01110 [astro-ph.HE] .
  * Iocco  _et al._ (2008) F. Iocco, A. Bressan, E. Ripamonti, R. Schneider, A. Ferrara, and P. Marigo, Mon. Not. Roy. Astron. Soc. 390, 1655 (2008), arXiv:0805.4016 [astro-ph] .
  * Leane and Smirnov (2022) R. K. Leane and J. Smirnov, (2022), arXiv:2209.09834 [hep-ph] .
  * Croon and Smirnov (2023) D. Croon and J. Smirnov, (2023), arXiv:2309.02495 [hep-ph] .
  * Leane and Smirnov (2023) R. K. Leane and J. Smirnov, (2023), arXiv:2309.00669 [hep-ph] .
  * Acevedo  _et al._ (2023) J. F. Acevedo, R. K. Leane, and J. Smirnov, (2023), arXiv:2303.01516 [hep-ph] .
  * Gould and Raffelt (1990) A. Gould and G. Raffelt, Astrophys. J. 352, 654 (1990).
  * Spergel and Press (1985) D. N. Spergel and W. H. Press, Astrophys. J. 294, 663 (1985).
  * Lopes and Lopes (2021) J. Lopes and I. Lopes, Astron. Astrophys. 651, A101 (2021), arXiv:2107.13885 [astro-ph.SR] .
  * Paxton  _et al._ (2011) B. Paxton, L. Bildsten, A. Dotter, F. Herwig, P. Lesaffre, and F. Timmes, Astrophys. J. Suppl. 192, 3 (2011), arXiv:1009.1622 [astro-ph.SR] .
  * Paxton  _et al._ (2013) B. Paxton _et al._ , Astrophys. J. Suppl. 208, 4 (2013), arXiv:1301.0319 [astro-ph.SR] .
  * Paxton  _et al._ (2015) B. Paxton _et al._ , Astrophys. J. Suppl. 220, 15 (2015), arXiv:1506.03146 [astro-ph.SR] .
  * Paxton  _et al._ (2018) B. Paxton _et al._ , Astrophys. J. Suppl. 234, 34 (2018), arXiv:1710.08424 [astro-ph.SR] .
  * Sakstein and Croon (2023) J. Sakstein and D. Croon, “Dark Matter Annihilation and Pair Instability Supernovae,” (2023).
  * Brott  _et al._ (2011) I. Brott, S. E. de Mink, M. Cantiello, N. Langer, A. de Koter, C. J. Evans, I. Hunter, C. Trundle, and J. S. Vink, Astron. Astrophys. 530, A115 (2011), arXiv:1102.0530 [astro-ph.SR] .
  * Gondolo and Silk (1999) P. Gondolo and J. Silk, Phys. Rev. Lett. 83, 1719 (1999), arXiv:astro-ph/9906391 .
  * Shen  _et al._ (2023) Z.-Q. Shen, G.-W. Yuan, C.-Z. Jiang, Y.-L. S. Tsai, Q. Yuan, and Y.-Z. Fan, (2023), arXiv:2303.09284 [astro-ph.GA] .