# Dark Black Holes in the Mass Gap

Nicolas Fernandez  Akshay Ghalsasi  Stefano Profumo  Nolan Smyth  and Lillian
Santos-Olmsted

###### Abstract

In the standard picture of stellar evolution, pair-instability – the energy
loss in stellar cores due to electron-positron pair production – is predicted
to prevent the collapse of massive stars into black holes with mass in the
range between approximately 50 and 130 solar masses – a range known as the
“black hole mass gap”. LIGO detection of black hole binary mergers containing
one or both black holes with masses in this mass gap thus challenges the
standard picture, possibly pointing to an unexpected merger history,
unanticipated or poorly understood astrophysical mechanisms, or new physics.
Here, we entertain the possibility that a “dark sector” exists, consisting of
dark electrons, dark protons, and electromagnetic-like interactions, but no
nuclear forces. Dark stars would inevitably form given such dark sector
constituents, possibly collapsing into black holes with masses within the mass
gap. We study in detail the cooling processes necessary for successful stellar
collapse in the dark sector and show that for suitable choices of the particle
masses, we indeed predict populating the mass gap with dark sector black
holes. In particular, we numerically find that the heavier of the two dark
sector massive particles cannot be lighter than, approximately, the visible
sector proton for the resulting dark sector black holes to have masses within
the mass gap. We discuss constraints on this scenario and how to test it with
future, larger black hole merger statistics.

##  1 Introduction

The momentous discovery of gravitational waves [1] ushered a new era in
astronomy and astrophysics with implications that could impact the fundamental
picture of particles and their interactions. As the statistics of binary
mergers detected via their gravitational radiation keep increasing [2],
population studies, especially of massive black holes, are becoming
increasingly statistically meaningful. One of the most puzzling findings at
present is the absence of the expected “mass gap” - the supposed absence of
black holes in the mass range between roughly 50 and 130 solar masses [3]. As
explained below, the mass gap is predicted for black holes with a standard
stellar collapse origin; while per se the absence of a mass gap might well
indicate the breakdown of some of the assumptions leading to its prediction
from stellar evolution, and/or an unexpected merger history of black hole
populations at early times, it might also point to a non-stellar origin for
the black holes populating the predicted gap. If that is the case, such black
holes might be primordial – not originating from stellar collapse, but rather
from the collapse of early density fluctuations dense enough that both the
size of the perturbation is larger than the Jeans scale and the particle
horizon is larger than the gravitational radius [4]. Forming primordial black
holes with masses in the “mass gap” is, however, very problematic, as the
required large density fluctuations would need to arise at very late times. As
such, the existence of “primordial” black holes in the mass gap might point to
an entirely different formation pathway, for instance the one we entertain
here (and that was proposed, in a different context, before, see e.g. [5]):
the black holes in the mass gap might originate from the collapse of “dark
stars” formed by the collapse of dark sector “electrons” and “protons” that
are not supported by any dark nuclear interactions. Before exploring this
admittedly rather unconventional scenario, we briefly review the theoretical
underpinning for the existence of a mass gap in the standard stellar collapse
picture.

The black hole mass gap (BHMG) is a result of the pair instability caused by
the production of electron-positron pairs in stellar cores, which converts a
portion of radiation pressure into massive particles, reducing the overall
outward pressure support the star [6]. The resulting implosion is reversed by
oxygen burning [7], leading to one of two possible outcomes: Stars with an
initial mass $70~{}M_{\odot}\lesssim M_{\rm{in}}\lesssim 140~{}M_{\odot}$
undergo a phase known as “pulsational pair-instability supernova”, consisting
of pulsations in which their cores contract, burn, expand, cool, and then
contract again, ejecting mass in each pulsation. It is generally believed that
these cores eventually collapse into black holes of mass
$35~{}M_{\odot}\lesssim M_{\rm{BH}}\lesssim 50~{}M_{\odot}$ [8]. For stars
with an initial mass $140~{}M_{\odot}\lesssim M_{\rm{in}}\lesssim
260~{}M_{\odot}$, the explosion following the implosion is so violent that the
entire star is disrupted and no remnant is left behind. This explosion is the
“pair-instability supernova” and it accounts for the existence of the mass
gap. Stars with initial masses $\gtrsim 260~{}M_{\odot}$ have insufficient
nuclear energy to counteract the implosion caused by the pair-instability.
Consequently, their cores collapse into black holes of mass $\gtrsim
130~{}M_{\odot}$ [9]. Given that the pair-instability supernova regime is
bounded both from above and from below by the fact that stars form black holes
outside this range, the generic expectation is not to observe any stellar
black holes of mass $50~{}M_{\odot}\lesssim M_{\rm{BH}}\lesssim
130~{}M_{\odot}$ [3].

Intriguingly, some of the merger events observed by LIGO in 2019 [10] have
indicated the existence of black holes in the mass gap with strong statistical
significance. The heavier of the original masses in event GW190521 is
$98.4^{+33.6}_{-21.7}\ ~{}M_{\odot}$ with the lighter mass of
$57.2^{+27.1}_{-21.7}\ ~{}M_{\odot}$; the merger masses for GW200220_061928
are $87^{+40}_{-23}\ ~{}M_{\odot}$ and $61^{+26}_{-25}\ ~{}M_{\odot}$.
Additionally, the primary black hole mass in the events GW190602_175927,
GW190706_222641 and GW190929_012149 are $71.8^{+18.1}_{-14.6}\ ~{}M_{\odot}$,
$74.0^{+20.1}_{-16.9}\ ~{}M_{\odot}$ and $66.3^{+21.6}_{-16.6}\ ~{}M_{\odot}$
respectively [10, 11]. As alluded to above, these events suggest that there
may be some unknown black hole formation mechanism not described by our
current models of stellar evolution. Among credible astrophysical mechanisms
that could fill the mass gap we mention here second-generation BHs emerging
from previous merger of two BHs [12, 13, 14, 15], super-Eddington accretion in
isolated binaries [16], or mergers between an evolved star and a main-sequence
companion [17, 18].

The work in explaining black holes in the BHMG using beyond Standard Model
physics has thus far primarily focused on additional cooling mechanisms in the
star [19, 20, 21, 22], which prevent the pulsation pair instability supernova.
In this paper we present an alternative. We consider the possibility that a
fraction of dark matter consists of a dissipative component, specifically a
“dark atom” with a heavy particle (we refer to as the “dark proton”) with mass
$m_{X}$, a lighter particle (“dark electron”) with mass $m_{c}<m_{X}$, and a
dark fine-structure constant $\alpha_{D}$. Such a complex dark sector is well
motivated by solutions to the hierarchy problem (e.g. Twin Higgs,
$N$-Naturalness) [23, 24]. However we will remain agnostic to the exact origin
of the dark sector.

We consider the cooling and fragmentation of such a dark sector and show that
this can lead to the formation of “dark stars”. We will use the term “stars”
to describe these dark gas clumps since they evolve according to the typical
equations of stellar evolution. However, we assume that there is no nuclear
physics in the dark sector and stress that our use of the word “stars” should
not be understood to imply the presence of nuclear reactions. We will
investigate the parameter space in the $m_{X},m_{c},\alpha_{D}$ plane that
will form black holes in the BHMG.

The likelihood of our atomic dark sector forming black holes within the BHMG
depends on the following considerations. As structure forms hierarchically,
the atomic dark matter falls into virialized halos and shock heats to the
virial temperature. If the virial temperature is above the ionization energy
of atoms, the dark atoms will ionize and cool through Compton scattering,
bremsstrahlung, and atomic transitions. If a significant amount of molecules
form during this process, then molecular processes can efficiently cool the
gas further. As the atomic dark sector cools, it fragments to smaller mass
scales. This continues until the opacity reaches a point where it cannot cool
anymore and the gas stops fragmenting. The smallest scale of fragmentation
reached when the gas reaches thermal equilibrium (TE) is called the “Jeans
Mass”, $M_{J}$.

The Jeans mass of the star is, in turn, controlled by the cooling mechanism
(atomic vs. molecular) that allows for the smallest Jeans mass at a given
point in the parameter space. Finally, once the star becomes opacity limited,
it reaches hydrostatic equilibrium and further evolution into either a black
hole or a degeneracy supported star (“white dwarf”) is determined by its mass.
Assuming the star is heavier than the Chandrasekhar mass, the star contracts
by radiating away energy from its surface until the core temperature becomes
comparable to the electron mass, at which point a heavy enough star will
experience pair instability and undergo dynamic collapse into a BH. We also
verify that our dark stars can form a black hole before $z=2$ and that the
stars will lose pressure support, i.e. the first adiabatic index,
$\Gamma_{1}$, drops below $4/3$ (see Section 4 and Appendix B.2).

We discuss the above considerations in detail in the remainder of this paper,
which is structured as follows: In Section 2 we describe the paradigm of
“atomic dark matter” and define reasonable boundaries on the dark electron
mass $m_{c}$, the dark proton mass $m_{X}$, and the dark fine-structure
constant, $\alpha_{D}$. We also determine the parameter space where cooling
via atomic cooling (including atomic transitions, Compton scattering and
bremsstrahlung) is efficient. In Section 3 we compute the Jeans mass of the
dark stars assuming atomic and molecular cooling, and derive the parameter
space that gives rise to black holes in the mass gap. In Section 4 we consider
the evolution of dark stars, and describe the parameter space where the stars
can cool efficiently within the age of the universe until dynamic collapse
sets in. We also derive the adiabatic index of the core to find the mass of
stars for a given parameter space that experiences core collapse. In Section 5
we discuss other constraints on the dark atom parameter space. We summarize
our findings and discuss future directions in Section 6.

##  2 Atomic Dark Matter

As discussed above, a fraction of dark matter can form atomic states
consisting of a heavy “proton” with mass $m_{X}$ and a light electron with
mass $m_{c}\ll m_{X}$ and we denote the fraction of atomic dark matter (ADM)
$\epsilon=\frac{\Omega_{ADM}}{\Omega_{DM}}$. Atomic dark matter and the
associated phenomenology has been explored extensively [25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 5, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52]. The upper bound on the fraction of atomic dark matter comes
from considering the impact of the inevitable dark acoustic oscillations on
the CMB as well as on the matter power spectrum [28, 34]. While the exact
constraints depend on the ratio of the dark CMB temperature to the visible-
sector CMB temperature, as well as on $m_{X},m_{c},\alpha_{D}$, all of the
parameter space can be accommodated as long as the dark atoms make up at most
$\approx 5\%$ of the total DM density. In Ref. [37] galaxy and star formation
from dark atoms was investigated and constraints on the parameter space were
derived assuming that $5\%$ of the dark matter is atomic. Strong constraints
from measurements of galaxy mass (from rotation curves), as well as dynamic
heating of ultra-faint dwarfs from black holes formed from the dark sector
rule out most of the interesting parameter space of atomic dark matter.
However in this paper we assume that atomic dark matter can generically be
much less than $5\%$ of the dark matter density, potentially avoiding the
above-mentioned constraints.

Let us now discuss the parameter space $(m_{X},m_{c},\alpha_{D})$. The dark
atoms fall into halos formed by dark matter, shock-heating to the virial
temperature of the halo $T_{\mathrm{vir}}$. If the virial temperature is
comparable or larger than the binding energy of the atom,
$B=\alpha^{2}_{D}m_{c}/2$, then the dark hydrogen ionizes and cools via
Compton scattering, bremsstrahlung, and atomic transitions. The cooling time
has to be fast enough that the gas cools before the halo experiences further
mergers which increase its virial temperature. Most of the parameter
constraints we mention below depend on the ability of the gas to cool
efficiently. We discuss the cooling and other relevant timescales in detail in
Appendix A and justify our parameter choice. Our parameter range choices are
as follows:

  * •

$\epsilon\in[5\times 10^{-4},5\times 10^{-2}]$. Here the upper bound of $5\%$
comes from the discussion above (and in [37]) while the lower bound comes from
the ability of the gas to equilibrate within the lifetime of the
universe111Compton scattering does not depend on $\epsilon$ and hence can be
always efficient in a large part of the parameter space, thus the constraint
comes from the ability of gas to equilibrate. See Fig. 1 and Appendix A for
details..

  * •

$\alpha_{D}\in[10^{-3},10^{-1}]$. Here the lower bound on $\alpha_{D}$ comes
from inefficient cooling. The upper bound comes from the fact that we can use
standard perturbation theory results for our hydrogen atoms and molecules as
well as keeping the binding energy low enough to allow the gas to ionize after
virialization.

  * •

$m_{c}\in[10^{-5}~{}\text{GeV},10^{-2}~{}\text{GeV}]$. Here the upper bound
stems from inefficient cooling and the need to keep the binding energy lower
than the virialization temperature. The lower bound comes from the need to
keep the electrons non-relativistic at high virial temperatures.

  * •

$m_{X}\in[1~{}\text{GeV},10~{}\text{GeV}]$. Here the lower bound comes from
requiring that the dark stars we want to form have masses in the mass gap
$\mathcal{O}(100)~{}M_{\odot}$ which is below the Chandrasekhar mass,
$M_{c}\approx 5.6~{}M_{\odot}\left(\frac{~{}\text{GeV}}{m_{X}}\right)^{2}$
222The Chandrasekhar mass is $5.6~{}M_{\odot}$ for pure hydrogen stars. if
$m_{X}\leq 0.1~{}\text{GeV}$ . The upper bound comes from considering the
thermal equilibration timescale between electrons and protons which can be
prohibitively large for $m_{X}\gg 10~{}\text{GeV}$.

  * •

Finally we fix the ratio of dark CMB temperature to CMB temperature to be
$r=\frac{T_{d,0}}{T_{0}}=0.5$ 333Note that $r\leq 0.5$ is required by current
CMB measurement. None of our results depend on the actual value of $r$ except
for a small part of the parameter space that cools by Compton cooling..

###  2.1 Atomic cooling

In order for atomic dark matter to cool efficiently, the gas has to be
partially ionized. Hence, the virial temperature has to satisfy
$T_{\mathrm{vir}}\geq\frac{B_{X}}{10}$ 444It is possible that for
$T_{\mathrm{vir}}<B_{X}/10$ the gas can cool through molecular cooling from
molecules formed during dark recombination. However molecular cooling is
quadrupole ($\alpha^{2}_{D}$) suppressed as well as
$\left(m_{c}/m_{X}\right)^{5}$ suppressed, compared to atomic cooling. Hence
we ignore this possibility.. During hierarchical structure formation, the
virial temperature of the halos grows with decreasing redshift. Cooling is
more efficient at larger redshifts because of the correspondingly larger
number density. Thus the cooling rate is maximal at the redshift $z$ when
$T_{\mathrm{vir}}\approx\frac{B_{X}}{10}$. In Appendix A we provide a semi-
analytical method to calculate the regions of our $\alpha_{D},m_{c}$ parameter
space that will cool at least $10\%$ of atomic dark matter efficiently for a
given $m_{X}$. Since the black holes observed by LIGO have $z<2$, we will
require the gas to cool by $z=2~{}(t=3.3{~{}\rm Gyr})$.

Figure 1: Left: Parameter space in which dark atoms can cool efficiently via
atomic transitions, bremsstrahlung and compton scattering for
$m_{X}=1~{}\text{GeV}$. Right: Parameter space in which dark atoms can cool
efficiently via atomic transitions, bremsstrahlung and compton scattering for
$m_{X}=10~{}\text{GeV}$.

###  2.2 Molecular Cooling

Dark molecules, the bound states of dark atoms, can cool the gas to even lower
temperatures. The ability for dark molecules to cool depends on the fraction
of dark molecules that can form. However, it is not straightforward to predict
the fraction of dark molecules that form in the universe 555For e.g. at low
redshifts hydrogen molecules form on dust grains in our universe.. Even taking
into account equilibrium processes, hydrogen molecule formation is an involved
network of creation and destruction processes that are difficult to
approximate analytically. Moreover $H^{-}$, which catalyzes the formation of
$H_{2}$, has a low binding energy of $B_{X}/10$ and can be easily destroyed by
radiation from the first dark stars. Important progress has been made recently
in calculating the formation of dark molecular hydrogen quantitatively [45,
46, 47]. Such a calculation needs to be embedded into a merger tree or N-body
simulation in order to accurately track the thermal history of dark atoms and
the corresponding formation rate of dark molecules, including ionizing
backgrounds from dark stars. However, such a calculation is beyond the scope
of this paper and will be left for future work.

Molecular cooling proceeds through two mechanisms, rotational cooling and
roto-vibrational cooling. Both mechanisms are quadrupole-suppressed for
$H_{2}$. The Einstein $A$ coefficients (rate of transitions) for both the
processes are given by [45]

$\displaystyle\Gamma_{\rm{rot}}$
$\displaystyle\simeq\alpha^{7}_{D}\frac{m^{6}_{c}}{m^{5}_{X}},$
$\displaystyle\Gamma_{\rm{ro-vib}}$
$\displaystyle\simeq\alpha^{7}_{D}\frac{m^{7/2}_{c}}{m^{5/2}_{X}}\,.$ (2.1)

We will assume the rates above to be equal to the quoted RHS in what follows
666Ideally the $\mathcal{O}(1)$ factor of the rates above can be deduced by
comparing to SM $H_{2}$ data. But for our order of magnitude estimates, this
does not have a large effect.. Note that molecular cooling is quadropole
suppressed as well as suppressed by powers of $\frac{m_{c}}{m_{X}}$ and is not
always efficient. The parameter space in which dark atoms can cool efficiently
will eventually evolve to form dark stars. In the next section we provide
details on the calculations needed to calculate the mass of dark stars.

##  3 Mass of Dark Stars – Jeans Mass

In this section we will discuss the collapse of a dark gas clump into a black
hole. Our method of calculating the mass of the dark star will be similar to
that for SM baryons. If the diffuse gas cools efficiently, it is continually
described by an isothermal distribution. The gas progressively loses pressure
support and collapses to higher densities. As the density increases, sound
waves travel to increasingly smaller distances within a free-fall time
$t_{\rm{dyn}}$ and encompass a smaller mass fraction of the dark clump. Thus
the dark gas clump fragments to larger densities (smaller Jeans mass) until
cooling becomes increasingly inefficient at higher densities. The Jeans radius
and mass are given by

$\displaystyle R_{J}$ $\displaystyle=\left(\frac{15T}{4\pi
Gm^{2}_{X}n_{X}}\right)^{1/2},$ $\displaystyle M_{J}$
$\displaystyle=\frac{4\pi}{3}m_{X}n_{X}R^{3}_{J}.$ (3.1)

Note that the gas fragments to smaller Jeans mass as long as the gas
temperature $T$ scales slower than $n^{1/3}_{X}$. In the absence of nuclear
processes, the final mass of the star, and hence the mass of the remnant black
hole, is likely to be close to the final Jeans mass, absent significant mass
loss mechanisms (i.e. we assume mass of dark star to be equal to the mass of
the black hole $M_{*}=M_{\rm BH}$ ). The final Jeans mass is set by the
temperature and density of the gas when cooling becomes inefficient. This can
happen if the gas reaches local thermal equilibrium or if the cooling becomes
opacity limited (i.e. the photons get absorbed well before they can escape the
gas cloud). We will discuss the Jeans mass for the cases of atomic cooling and
molecular cooling, detailing the area where we get $M_{\rm
BH}=M_{J}=100~{}M_{\odot}$. We choose this value because it is squarely within
the BHMG, and our constraints are calculated assuming a $100M_{\odot}$ dark
star. It should be noted that the following calculation of the Jeans mass
serves as an order of magnitude estimate at an $\mathcal{O}(10)$ level due to
the uncertainties associated with the Jeans mass calculation, the inexact
correspondence between $M_{J}$ and $M_{\rm BH}$ and ignoring accretion that
can increase the mass of the BH. Thus we consider $30-300M_{\odot}$ Jeans mass
as the range where we can reasonably expect to create a $100M_{\odot}$ BH.

###  3.1 Atomic Jeans Mass

As noted above, the Jeans mass is set by the temperature and density of the
gas when cooling becomes inefficient. Cooling can become inefficient if the
gas reaches thermal equilibrium. However for atomic cooling the timescale for
the lowest-energy transition $(2p\rightarrow 1s)$ is extremely short, making
the density at which thermal equilibrium is reached very large. Cooling also
becomes inefficient if the opacity destroys photons escaping the gas cloud, a
phenomenon which transitions the cloud from efficient volume cooling to
inefficient surface cooling. In all of the parameter space we consider,
opacity is set by free-free absorption (inverse process of bremsstrahlung)
777The atomic transition from ($2p\rightarrow 1s$) is fast and opacity due to
this transition is only effective at very large densities and is not relevant
to the Jeans mass calculation for our parameter space.. The Rosseland mean
opacity of an ionized gas for free-free absorption is given by

$\displaystyle\alpha_{R,ff}$
$\displaystyle=0.38\times\frac{\alpha^{3}_{D}n^{2}_{X}}{m^{3/2}_{c}T^{7/2}}$
$\displaystyle=1.3\times
10^{4}\times\frac{n^{2}_{X}}{m_{c}^{5}\alpha^{4}_{D}}\,,$ (3.2)

where in the last line we have assumed that the gas can cool down efficiently
to $T=\frac{B_{X}}{10}$, where $B_{X}=\frac{\alpha^{2}_{D}m_{c}}{2}$ is the
binding energy of dark atoms. We approximate the Jeans mass by setting
$\alpha_{R,ff}R_{J}=1$, i.e. the mean free path of the photons is the size of
the gas cloud. Using Eqs. 3 we get

$\displaystyle M_{J}\simeq
200~{}M_{\odot}\left(\frac{\alpha_{D}}{\alpha}\right)^{2}\left(\frac{m_{p}}{m_{X}}\right)^{7/3}\,.$
(3.3)

Note that if the cooling is atomic then the Jeans mass is independent of the
electron mass. If we had assumed that the gas can cool down efficiently till
$B_{X}$ instead of $0.1B_{X}$ we would have gotten Jeans mass higher by a
factor of $\approx 3$ reflecting the uncertainties associated with calculating
the Jean’s mass. Thus as mentioned above the Jeans calculation is valid to
$\mathcal{O}(10)$ level. In the plots below we show the region that accounts
for black holes with mass $M\in(30~{}M_{\odot}-300~{}M_{\odot})$. This target
range is shown in Fig. 2 (blue band).

###  3.2 Molecular Jeans Mass

Molecular cooling becomes inefficient once the gas reaches local thermal
equilibrium (LTE). Since the quadropole transition timescales for molecular
cooling are much larger compared to atomic cooling, the gas reaches LTE at
relatively low densities. The density at which LTE is achieved can be
calculated by equating the scattering rate $(n_{LTE}\sigma v)$ (where
$v\simeq\sqrt{T/m_{X}}$) to the transition rate of $H_{2}$ molecules given in
Eq. 2.2. For the purposes of our calculation we will assume a geometric cross
section $\sigma=\pi a^{2}_{0}$ where $a_{0}$ is the Bohr radius. For a more
detailed calculation see Ref. [45]. The final temperature that the gas reaches
depends on the dominant cooling mechanism and is given by

$\displaystyle T_{f,rot}$
$\displaystyle=\frac{\alpha^{2}_{D}m^{2}_{c}}{m_{X}},$ $\displaystyle T_{f,ro-
vib}$ $\displaystyle=\frac{\alpha^{2}_{D}m^{3/2}_{c}}{m^{1/2}_{X}}\,.$ (3.4)

Then using Eq. 2.2 and Eq. 3.2, and setting the collisional timescale equal to
the molecular de-excitation timescale we get

$\displaystyle n_{LTE,rot}=10^{3}{~{}\rm
cm^{-3}}\left(\frac{\alpha_{D}}{\alpha}\right)^{8}\left(\frac{m_{c}}{m_{e}}\right)^{7}\left(\frac{m_{p}}{m_{X}}\right)^{4},$
$\displaystyle n_{LTE,ro-vib}=2\times 10^{8}{~{}\rm
cm^{-3}}\left(\frac{\alpha_{D}}{\alpha}\right)^{8}\left(\frac{m_{c}}{m_{e}}\right)^{19/4}\left(\frac{m_{p}}{m_{X}}\right)^{7/4}\,,$
(3.5)

where we have fixed the coefficient to match the SM values for $n_{LTE,rot}$
in [53]. Then using Eq. 3 and Eq. 3.2 we get that

$\displaystyle M_{J,rot}$ $\displaystyle\simeq 5\times
10^{3}~{}M_{\odot}\left(\frac{m_{e}}{m_{c}}\right)^{1/2}\left(\frac{\alpha}{\alpha_{D}}\right)\left(\frac{m_{p}}{m_{X}}\right)^{3/2},$
$\displaystyle M_{J,ro-vib}$ $\displaystyle\simeq 3\times
10^{3}~{}M_{\odot}\left(\frac{m_{e}}{m_{c}}\right)^{1/8}\left(\frac{\alpha}{\alpha_{D}}\right)\left(\frac{m_{p}}{m_{X}}\right)^{15/8}\,.$
(3.6)

The gas cloud will fragment within dynamical time to the smallest allowed
Jeans mass. For $m_{X}=1~{}\text{GeV}$, comparing the Jeans molecular mass Eq.
3.2 to the atomic Jeans mass Eq. 3.3, we can show that the atomic Jeans mass
is the lowest in our desired target range of $30-300~{}M_{\odot}$ and our
allowed parameter space. For $m_{X}=10~{}\text{GeV}$, the minimum between
atomic and molecular Jean’s mass is lower than $30~{}M_{\odot}$. Thus the only
way we get black holes in the target range of $30-300~{}M_{\odot}$ is assuming
molecular cooling is inefficient. Thus we only consider the atomic Jeans mass
in Eq 3.3.

##  4 Stellar Evolution

The dark stars form with masses given by the Jeans mass calculated in the
above section. Under the assumption that the dark gas clump is in hydrostatic
equilibrium, we can use the stellar equations to describe the evolution of the
system. Since we assume no nuclear physics in the dark sector, the dark stars
are purely gravity powered, i.e. the gravitational potential energy lost due
to contraction is radiated away by photons. We assume nuclear reactions are
absent in the dark sector. To begin with we also assume that there is no
neutrino-like energy loss mechanism. However, as we will see, some additional
energy loss mechanism (see e.g. [19, 20]) or the presence of heavier dark
elements is needed if we want a star with $m_{X}=1~{}\text{GeV}$ to collapse
into a BH within the mass gap before $z=2$.

The absence of nuclear reactions makes it possible to study the evolution of
dark stars analytically. In order to understand the stellar evolution and
possible collapse into a BH, we need to describe the central density
$\rho_{c}$ as well as central temperature $T_{c}$ of the star. Once $T_{c}$
becomes comparable to the mass of the dark electron $m_{c}$, pair instability
can set in, triggering the collapse of the star. Since stellar cooling is due
to radiation from the surface, we also need to know the surface temperature
$T_{s}$, so we can confirm that the star can collapse till the central
temperature reaches the pair instability region by $z=2$. We can examine
$\rho_{c},\ T_{c}$ and $T_{s}$ by modelling the star as a polytrope. A
polytropic profile implies a global relation between pressure and density of
the form $P=\rho^{\gamma}$, where $\gamma=(n+1)/n$. This simple treatment for
solving systems in hydrostatic equilibrium is a surprisingly powerful tool for
understanding stellar properties. In particular, massive stars are very close
to a $n=3$ polytrope. This is because they are primarily supported by
radiation pressure and ion pressure. Assuming that relative contribution from
each, parameterized here as $\beta\equiv P_{\rm{gas}}/P$, is constant
throughout the star, the total pressure scales as $P\sim\rho^{4/3}$,
corresponding to a polytrope with $n=3$. We find reasonably good agreement
between the $T_{c}$, $T_{s}$, and $\rho_{c}$ values calculated using the $n=3$
polytrope and those calculated using MESA (see Figure 4 and Appendix B for
details). The central density and temperature of the star evolve as

$T_{c}=\xi\rho^{1/3}_{c}\,,$

where $\xi=\left(\frac{45(1-\beta)}{\pi^{2}(\mu m_{X})\beta}\right)^{1/3}$,
and $\mu$ is the mean molecular weight. The derivation of $T_{c},\xi,T_{s}$ as
well as the derivation of the mass of the stars that reach the pair
instability region is given in App. B. Here we just quote the results. The
potential energy of the dark star that needs to be radiated away for a $n=3$
polytrope is $E_{*}=\frac{3}{2}\frac{GM^{2}_{*}}{R}$ where $R$ is the radius
of the star. The temperature of the star is calculated in Eq. B.10 and is
given by

$T_{s}=9.7\times
10^{-4}~{}\text{GeV}\left(\frac{100}{f}\right)^{2/7}\left(\frac{\alpha}{\alpha_{D}}\right)^{6/7}\left(\frac{m_{c}}{m_{e}}\right)^{3/7}\left(\frac{m_{X}}{m_{p}}\right)^{2/21}\\\
\times\left(\frac{(1-\beta)^{13}}{\beta^{10}}\right)^{1/21}\left(\frac{\rho_{c}}{10^{-15}~{}\text{GeV}^{4}}\right)^{4/21}\,,$
(4.1)

where $f$ is the ratio of bound-free to free-free opacity (see App. B). Thus
the cooling time, defined as the time it takes the star to reach temperatures
where pair instability can set in i.e. $T_{c}=m_{c}/10$, is given by

$\displaystyle t_{\rm{PI}}=\frac{E*}{4\pi
T^{4}_{s}R^{2}}=\frac{\frac{3GM^{2}_{*}}{2R}}{4\pi T^{4}_{s}R^{2}}\,.$ (4.2)

The cooling time scales as $t_{PI}\propto\rho_{c}^{5/21}$ i.e. it gets harder
to cool the star as central density increases. Thus, majority of the total
evolution time spent reaching $T_{c}$ is dominated by the period with the
highest temperature, $T\simeq T_{c}\simeq m_{c}/10$ (i.e.
$\rho_{c}=\left(T_{c}/\xi\right)^{3}$), after which pair instability can
occur. The time until pair instability sets in is given by

$\displaystyle t_{\rm{PI}}=5.3\times
10^{4}\,\mathrm{yrs}\left(\frac{f}{100}\right)^{8/7}\left(\frac{\alpha_{D}}{\alpha}\right)^{24/7}\left(\frac{m_{e}}{m_{c}}\right)\left(\frac{m_{X}}{m_{p}}\right)^{13/7}\left(\frac{M_{*}}{100~{}M_{\odot}}\right)^{2}\left(\frac{\beta^{29/7}}{(1-\beta)^{45/14}}\right)\,.$
(4.3)

Note that $\beta$ also has a dependence on $m_{X}$ and $M_{*}$ (Eq. B.6). In
principle we need to evaluate $t_{\rm{PI}}$ assuming the star mass to be the
Jeans mass i.e. $M_{*}=M_{J}$. However to denote our constraint, we will use
$M_{*}=100~{}M_{\odot}$ i.e. a value inside the BHMG. In Fig. 2 we show the
associated constraint from $t_{\rm{PI}}<3.3\mathrm{Gyr}~{}(z=2)$ (green
region). The constraint should be treated as approximate since it’s a near
linear function of our approximate factor $f$. This constraint does not rule
out our region of interest of atomic Jean’s mass between $30-300M_{\odot}$ for
$m_{X}=1~{}\text{GeV}$. and does not affect any parameter space for
$m_{X}=10~{}\text{GeV}$. Additional cooling mechanisms in the dark atomic
sector can significantly reduce $t_{\rm{PI}}$, however we don’t consider them
here. Even when $T_{c}$ reaches temperatures where pair production occurs,
collapse through pair instability is not guaranteed. The question of dynamical
stability through pair-production motivates introducing the first adiabatic
index, given by

$\Gamma_{1}=\frac{\rho}{P}\Big{(}\frac{\partial
P}{\partial\rho}\Big{)}_{s}\,,$ (4.4)

the logarithmic derivative of pressure with respect to density at fixed
entropy. It can be shown, assuming homologous and adiabatic compression, that
dynamical stability requires $\Gamma_{1}\gtrsim 4/3$. This can be seen as
follows: the pressure on a unit area of a concentric sphere in hydrostatic
equilibrium is given by

$P=\int_{m}^{M}\frac{Gm}{4\pi r^{4}}dm\,.$ (4.5)

If we compress the star adiabatically and assume homology (self-similarity of
stellar evolution at all radii), then a radial shell being compressed from $R$
to $R^{\prime}$ will change the right hand side as
$\sim\Big{(}\frac{R^{\prime}}{R}\Big{)}^{-4}$, whereas the pressure will scale
as
$\Big{(}\frac{\rho^{\prime}}{\rho}\Big{)}^{\Gamma_{1}}\sim\Big{(}\frac{R^{\prime}}{R}\Big{)}^{-3\Gamma_{1}}$.
Thus, if $\Gamma_{1}\lesssim 4/3$, the weight on the sphere will increase
faster than the counteracting pressure and the sphere will collapse. Note that
this condition is approximate because small corrections may be present due to
relativistic effects. However, these corrections are typically only relevant
for neutron stars, hence why we will assume a critical value of
$\Gamma_{1}=4/3$. Here we consider $M_{\rm{crit}}$ to be the lowest mass for
which the core experiences pair instability. We have solved for this
instability in App. B numerically with the result 888Ideally, the collapse
condition should be that the average $\Gamma_{1}$ over the entire star is
below $4/3$. However, we only consider the local value of $\Gamma_{1}$ in the
core. Thus $M_{\rm{crit}}$ derived here is likely smaller than the true value.

$\displaystyle M_{\rm{crit}}\simeq
260~{}M_{\odot}\frac{1~{}\text{GeV}}{m_{X}}\,.$ (4.6)

Note that this gives a lower bound on the mass of the dark proton, if we
require the dark stars to form black holes in the mass gap. For $M_{J}\gtrsim
M_{\rm{crit}}$ the star will experience pair instability and collapse. Note
that most of the region for $m_{X}=1~{}\text{GeV}$ does not experience pair
instability for a $M_{*}=100~{}M_{\odot}$ and hence will not collapse into a
BH within $3.3~{}\rm{Gyr}$. Collapse can be achieved by assuming additional
cooling mechanisms or the presence of a significant fraction of “dark helium”
or heavier “dark elements”, at the expense of minimality of the dark sector.

Figure 2: Left: Target region where we can expect $100M_{\odot}$ BH (blue
band) and constraint for core temperature of star to reach pair instability
temperatures before z = 2 (3.3 Gyr) (green region) for $m_{X}=1\
~{}\text{GeV}$. Right: Target region where we can expect $100M_{\odot}$ BH
(blue band) for $m_{X}=10~{}\text{GeV}$. The $m_{X}=10\ ~{}\text{GeV}$ case is
not affected by this constraint coming from time to reach pair instability.

##  5 Constraints and Results

If dark atoms make up a fraction of dark matter, it is subject to several
constraints from early universe (CMB and matter power spectrum) to late
universe astrophysics (galaxy morphology, ultra-faint dwarfs, lensing etc.).
Most of these constraints are independent of whether the dark atoms collapse
eventually into black holes.

###  5.1 CMB and $N_{eff}$

Dark atoms undergo dark acoustic oscillations, analogous to the baryons in SM
and the dark photons contribute to $N_{eff}$. These constraints have been
derived in [28, 34] and the parameter space we consider here already satisfies
these constraints by assuming $\epsilon\leq 0.05$ and
$r=\frac{T_{d,0}}{T_{0}}\leq 0.5$.

###  5.2 Galaxy Morphology

The primary constraint on the fraction of dark matter constituted by dark
atoms comes from the morphology of dark galaxies. For $\epsilon=0.05$, dark
atoms typically form a bulge-like galaxy at the center of the Milky Way [37].
The mass of our galaxy is well measured with stellar rotation curves as well
as from measurements of stellar luminosity, and the two measurements agree up
to a factor of 2. Thus, atomic dark matter in the bulge cannot contribute more
than the mass of the stars in the center of the Milky Way. This gives the
constraint [37]

$\displaystyle\epsilon\times
f_{bulge}\left(\frac{r_{MW}}{r_{s}}\right)^{3}<0.02\,,$ (5.1)

where $r_{MW}=3\ \rm{kpc}$ is the stellar radius of the Milky Way and $r_{s}$
is the radius of the dark galaxy. For dark galaxies which form a bulge i.e.
$f_{bulge}\simeq 1$ instead of the disk $r_{s}\ll r_{MW}$ providing a
stringent constraint of $\epsilon$. Note that $r_{s}$ is a function of
$\epsilon,\alpha_{D},m_{c},m_{X}$ since it dictates the merger history and
hence morphology of the dark galaxy. For $\epsilon=0.05$, a large part of the
parameter space is ruled out and the constraints derived in [37] are plotted
in Fig. 3 (dotted brown) 999The discretized nature of the constraint is a
limitation arising from the fact that merger tree simulations have been
performed for discretized grid in $\alpha_{D},m_{c}$ parameter space.. It is
possible that $\epsilon=5\times 10^{-4}$ evades this constraint, however an
independent set of merger tree simulations need to be performed with
$\epsilon=5\times 10^{-4}$. We leave this for future work. If the dark atoms
form a disk instead of a bulge, an equivalent constraint of $\epsilon<0.05$
comes from Gaia measurements of the Milky Way surface density [54]. However
since most of the parameter space for $\epsilon=0.05$ does not form a disk, we
do not plot these constraints here.

###  5.3 Other Constraints

We now discuss constraints that don’t affect our scenario or require
additional modelling in order to derive them.

Dynamical Friction causes the heavier dark stars to sink to the center of
dwarf galaxies, imparting their momentum to SM stars resulting in “dynamical
heating” of SM stars and increasing the half light radii of the ultra faint
dwarfs [55]. This gives the constraint

$\epsilon\times f_{cooled}\times M_{*}<10~{}M_{\odot}\,,$

where $f_{cooled}$ is the fraction of gas in the halo that has cooled and
formed stars. However since we consider $\epsilon<0.05$ and $M_{*}\simeq
100~{}M_{\odot}$ in the mass gap, this constrain is not applicable to us.

Microlensing surveys like MACHO [56] and OGLE [57] typically put strong
constraints on sub solar mass black holes and are not relevant for our target
parameter space of $M_{*}\simeq 100~{}M_{\odot}$. See [37] for further
details.

Presence of dark atoms inside a SM star modify the stellar mass luminosity
relation [52]. These constrain the mass fraction of dark atoms inside star to
be $\lesssim 5\%$. However since it’s not straightforward to translate the
total fraction of dark matter that is atomic $(\epsilon)$ to the fraction that
ends up in SM stars we won’t consider this constraint.

Similarly, the presence of baryons inside a dark star can also lead to
observable signatures [48]. However, these constraints depend on the kinetic
mixing between dark and SM photon, which we have ignored here.

Figure 3: Left: Constraint plots including the Galaxy mass and morphology
constraints (starred blue) from [37] for $m_{X}=1~{}\text{GeV}$ and
$\epsilon=0.05$. Right: Constraint plots including the Galaxy mass and
morphology constraints (starred blue) from [37] for $m_{X}=10~{}\text{GeV}$
and $\epsilon=0.05$. Note that the “discretized” nature of the constrains
comes from limited merger tree simulation data in the $\alpha_{D},m_{c}$
plane.

##  6 Discussion, Conclusions, and Future Work

Observations made by LIGO of merging binary black holes with masses within the
BHMG are in direct conflict with our understanding of standard model stars.
Electron-positron pair instability within the cores of SM stars reduces the
supportive radiation pressure and leads to a runaway collapse and subsequent
thermonuclear burning of oxygen, causing a destructive explosion that causes a
star to shed mass or leaves no remnant black hole at all. This motivates a
scenario in which dark analogues of the electron, proton, and photon are able
to dissipate energy through radiation but are unable to generate energy
through fusion. In the absence of nuclear reactions, the radiated energy of a
dark star is sourced by gravitational collapse alone and runaway collapse can
directly form a black hole from the compressed dark stellar matter. In this
case, the instability of a dark star does not lead to a violent explosion due
to a pair instability supernova (PISN). Thus, there is no analogous BHMG for
dark stars. We analytically and numerically explored the best-motivated
regions of parameter space to determine the conditions in which black holes
would form in the mass gap. We showed that under reasonable assumptions about
the dominant atomic and molecular cooling mechanisms, the Jeans mass of dark
stars can fall within the BHMG. Furthermore, assuming a polytropic model and
comparing with state-of-the-art stellar evolution code, we demonstrated that
dark stars can reach a stage of dark electron-positron production so long as
$m_{X}\simeq 10\ ~{}\text{GeV}$, leading to dynamic instability and collapse.
With no fusion to support these dark stars, they can naturally form black
holes with masses in the BHMG before $z=2$, allowing time for these black
holes to merge and be detectable by LIGO. Intriguingly, we discovered that the
black holes that form from dark sector collapse would not have mass within the
mass gap if the dark sector proton were lighter than the visible sector
proton.

While we present a potential explanation for black holes in the BHMG, further
work remains to be done to associate the black holes in the mass gap found by
LIGO to black holes formed by dark atoms. Since the density and merger rates
of such black holes, as well as constraints on dark atoms, depend strongly on
the dark galaxy morphology, it is essential to perform merger tree simulations
similar to [37] for a larger grid of $m_{X},m_{c},\alpha_{D},\epsilon$, and
$\xi$. Specifically, improvements need to be made to the merger tree
simulations by including molecular cooling, which could be important in the
large $m_{c},\alpha_{D}$ parameter space. Effects of stellar feedback could
also be important in the small $m_{c},\alpha_{D}$ parameter space.
Furthermore, it is essential to compare $N$-body simulations to merger tree
simulations in order to have a more accurate understanding of the galactic
astrophysics of atomic dark matter. For simplicity, we assumed in this work
that the mass of the remnant black hole formed from dark star collapse is
approximately equal to the Jeans mass of the star. This is because there is no
supernova expelling mass from the system. However, a full $N$-body simulation
and/or thorough accretion modeling would be needed in order to more accurately
determine the black hole’s final mass. Modifying existing stellar evolution
code, such as MESA, to include different $\alpha_{D},m_{c}$ and $m_{X}$ will
also let us compare our analytical estimates of stellar evolution to the full
numerical solution. If there exists an appreciable population of black holes
formed from the collapse of dark stars, and if one can predict the morphology
of the dark galaxies in which the black holes reside, one can predict the
merger rate of black holes in the mass gap. Current and future observations
from LIGO can then be used to place limits on the abundance of black holes
formed through this mechanism and consequently on the dark atom parameter
space.

Finally one may imagine adding nuclear physics to the dark sector. In this
case, there may exist a dark BHMG with the upper and lower limits being
determined by physics of dark nuclear reactions. It is possible that for
different values of $m_{X}$ and $m_{c}$, there is a dark BHMG that does not
fully or even partially overlap with that of the standard model. In this case,
black holes could form over a full continuum of masses with some being created
by SM stars and some created by dark stars. While we have focused in this work
on the mass range relevant for the pair-instability mass gap of the standard
model, dark stars could form black holes with much larger masses. In
particular, this mechanism could contribute to seeding the supermassive black
holes that are found at the center of large galaxies. On the other end of the
spectrum, since the Chandrasekhar mass limit scales as $\sim 1/m_{X}^{2}$, it
is possible that dark stars formed of atomic dark matter could also form sub-
solar mass black holes, as envisioned in Ref. [5].

## Acknowledgments

We are very grateful to Hiren Patel for his help in the early stages of this
project. We also thank Matt McQuinn for many useful conversations. NS and AG
would like to thank Stan Woosley for help with questions in stellar evolution.
NS and LS would like to thank the MESA community for helping with MESA
modifications as well as troubleshooting. This material is based upon work
supported in part by the National Science Foundation Graduate Research
Fellowship under Grant No. DGE-1842400 to NS. This work is partly supported by
the U.S. Department of Energy grant number de-sc0010107 (SP).The work of AG is
supported by the U.S. Department of Energy under grant No. DE–SC0007914. The
work of NF is supported in part by DOE CAREER grant DESC0017840. NF would like
to thank the Aspen Center for Physics, which is supported by National Science
Foundation grant PHY-1607611, where this work was partly performed. LS
received funding through the Ron Ruby Scholarship at UC Santa Cruz.

##  Appendix A Cooling times

###  A.1 Dark Atomic Cooling

After dark recombination, the atomic dark matter falls into the gravitational
potential wells created by the rest of the non-atomic DM halo during formation
of early galaxies and virializes. If the virial temperature is large enough to
ionize the gas or excite atomic transitions, then the atomic dark matter will
cool through bremsstrahlung, compton scattering, and atomic transitions. The
timescales for each of these processes are given by

$\displaystyle t_{\rm brem}$
$\displaystyle=\frac{9}{2^{5}}\left(\frac{3\pi}{2}\right)^{1/2}\left(\frac{m^{3/2}_{c}T^{1/2}}{\alpha^{3}_{D}n_{c}(z)x^{2}}\right)\,,$
(A.1) $\displaystyle t_{\rm comp}$
$\displaystyle=\frac{135}{64\pi^{2}x}\left(\frac{m^{3}_{c}}{\alpha^{2}_{D}\left(T_{d0}(1+z)\right)^{4}}\right)\,,$
(A.2) $\displaystyle t_{\rm atomic}(T\gtrsim B_{X}/10)$
$\displaystyle=\frac{9}{2^{5}}\left(\frac{3\pi}{2}\right)^{1/2}\left(\frac{T^{2}}{\alpha^{6}_{D}n(z)}\right),$
(A.3)

where $x$ is the ionization fraction, $T$ is the gas temperature, $n$ is the
atomic DM number density and $T_{d0}$ is the dark CMB temperature today. The
atomic cooling timescale here is only referenced for $T\gtrsim B_{X}/10$, a
more general formula is given in Appendix A of [37]. The number density at
redshift $z$ is given by

$\displaystyle n(z)=\frac{\epsilon\rho_{m,0}\,\Delta(1+z)^{3}}{m_{X}},$ (A.4)

where $\rho_{m,0}$ is the dark matter density today, $\epsilon$ is the
fraction of DM made up of atomic dark matter, and $\Delta\approx 180$ is the
overdensity associated with a virialized halo compared to the background dark
matter density. The transfer of energy (and hence the thermal equilibrium
between the dark protons and electrons) is controlled by the timescale
associated with particle-particle collisions and is given by

$\displaystyle t_{\rm
eq}=\frac{m_{X}m_{c}}{2\sqrt{3}\alpha^{2}_{D}n}\left(\frac{3T}{m_{c}}\right)^{3/2}\left(\log\left(1+\frac{T^{2}}{4\alpha^{2}_{D}n^{2/3}}\right)\right)^{-1}\,.$
(A.5)

The maximum temperature that the gas reaches is the virial temperature of the
halo given by

$\displaystyle
T_{\mathrm{vir}}=\frac{1}{3}\left(\frac{4\pi}{3}\right)^{1/3}GM^{2/3}_{halo}m^{4/3}_{X}\epsilon^{-1/3}n(z)^{1/3}\,.$
(A.6)

Note that the virial temperature is independent of $\epsilon$ but scales as
$T_{\mathrm{vir}}\propto m_{X}$. Thus the equilibration timescale goes as
$t_{\rm eq}\propto m^{7/2}_{X}$ for the halo of the same mass. Atomic dark
matter $m_{X}\gg 10\hskip 2.84526pt~{}\text{GeV}$ cannot cool in equilibrium,
and we therefore ignore that possibility. Note that for $T>B_{X}$,
$\frac{t_{\rm{atomic}}}{t_{\rm{brem}}}\simeq\left(\frac{T}{B_{X}}\right)^{3/2}$
and thus atomic cooling is more efficient for $T\approx B_{X}/10$, with
bremsstrahlung becoming more efficient for $T\gg B_{X}$. We will require that
a significant fraction of the gas is able to virialize and cool by
$z=2$101010The BH mergers for BH in the mass gap are estimated to be at
$z\approx 1$. However since the gas needs to cool, form stars, evolve to
collapse into a BH and then undergo a binary merger. Thus we take the
benchmark of $z=2$ for our calculations.. We require that
$T_{\mathrm{vir}}>\frac{B_{X}}{10}$ for the majority of dark atoms. Imposing
the above condition and using Eq. A.6 we get

$\displaystyle M_{\rm halo}$ $\displaystyle\geq 2\times
10^{20}\left(\frac{m_{c}}{m_{X}}\right)^{3/2}\alpha^{3}_{D}~{}M_{\odot}=5.7\times
10^{20}\left(\frac{B_{X}}{m_{X}}\right)^{3/2}M_{\odot}\,.$ (A.7)

The cumulative mass fraction of halos above a given mass is given by [58, 59,
60]

$\displaystyle
P(>M,z=2)=\mathrm{erfc}\left(\frac{\delta_{c}(z=2)}{\sqrt{2}\sigma(M)}\right)\geq
0.1\,,$ (A.8)

where we assume that at least $10\%$ of the dark matter should be in halos
that have temperatures higher than $B_{X}/10$ and hence can potentially cool.
Note that the constraints we will derive are sensitive to the fractions we
assume. Throughout, we use a conservative value of $10\%$. Above,
$\delta_{c}(z)=1.686\frac{gf(0)}{gf(z)}$, where $gf(z)$ is the growth function
of matter perturbations given as

$\displaystyle
g(z)=\frac{5}{2}\Omega_{m}\frac{H(z)}{100h}\int_{0}^{1/(1+z)}\left(\frac{100h}{aH(a)}\right)^{3}da\,,$
(A.9)

and where $\sigma(M)$ is the standard deviation of matter overdensities, given
by

$\displaystyle\sigma^{2}(M)=\int
dk\frac{k^{2}}{2\pi^{2}}P(k)\left(\frac{j_{1}(kR)}{kR}\right)^{2},$ (A.10)

where $P(k)$ is the linear power spectrum at $z=0$ and $j_{1}$ is the first
spherical Bessel function and $R=\left(\frac{3M}{4\pi\rho_{m}}\right)^{1/3}$
is the Lagrange radius of a halo with mass $M$. Using Eqs. A.8, A.9, A.10 we
get that

$\sigma(M)=\frac{\delta_{c}(z=2)}{\sqrt{2}\times 1.15}=2.37,$

where we have used $\Omega_{m}=0.26,h=0.7$ in the above calculation. Using Eq.
A.10, this corresponds to a halo mass of $M\simeq 3\times
10^{11}~{}M_{\odot}$. If the lower bound on halo mass given in Eq. A.7 is
larger than halo mass of $3\times 10^{11}~{}M_{\odot}$ derived above, the
fraction of gas that can cool efficiently is less than $10\%$. This gives us
the constraint

$\displaystyle B_{X}=\frac{\alpha^{2}_{D}m_{c}}{2}<5.4\times 10^{-7}m_{X}\,,$
(A.11)

which is plotted in Fig. 1. The parameter space not bounded by the above
constraint can cool efficiently as long as

${\rm
min}(t_{\rm{brehm}},t_{\rm{compton}},t_{\rm{atomic}})<t_{\rm{dyn}}\qquad\text{and}\qquad
t_{equil}<t_{\rm{dyn}}$

for a given halo. Here, $t_{\rm dyn}=\frac{1}{\sqrt{8\pi G\rho_{m}}}$ is the
approximate timescale associated with doubling of the halo mass through
mergers and accretion. Different processes are efficient at cooling the halos
at different times and the exact fraction of halos that cool has to be
understood through a merger tree [37] or a N-body simulation. We will however
try to get an analytic understanding of the parametric space that can cool
efficiently. From Eqs. A.1 we understand that cooling is most efficient for
high $z$ as long as the gas can be ionized i.e. $T\gtrsim B_{X}/10$. Cooling
is also efficient for low virial temperature (except for compton cooling which
does not depend on the gas temperature). Thus once $T_{\mathrm{vir}}\approx
B_{X}/10$, the gas can cool. Since virial temperature increases with
decreasing redshift $z$, one can find a $z_{\rm max}$ for which at least
$10\%$ of the gas has a temperature above $B_{X}/10$. Then we can compare
cooling and equilibration times to the dynamical time of the halo to find out
if cooling is efficient. Here is the algorithm we follow:

  * •

Solve $M(z)$ s.t. $P(>M,z)=0.1$. This corresponds to solving

$\displaystyle\sigma(M(z))=\frac{\delta_{c}(z)}{\sqrt{2}\times 1.15}\,,$
(A.12)

  * •

with $M_{\rm halo}=M(z)$ solve for $z_{\rm max}$ using

$\displaystyle T_{\mathrm{vir}}(z)=\frac{B_{X}}{10}\,.$ (A.13)

  * •

Check if the cooling and equilibration times are smaller than the dynamical
timescale of the halo.

Following the above steps, the parameter space that can cool efficiently is
shown in Fig. 1.

##  Appendix B $n=3$ polytrope evolution and collapse

In this section we will discuss the evolution of the star assuming the star is
a $n=3$ polytrope. The $n=3$ polytropic solution will be applicable when the
core temperature $T_{c}<m_{c}$ i.e. we will investigate the early stages of
our star’s evolution. In order for a polytropic solution to apply, we will
assume that our dark star is in hydrostatic equilibrium. We will also assume
that we are in a regime (the star is heavy enough) that the pressure inside
the star is a combination of non-relativistic gas pressure and radiation (i.e.
we avoid electron degeneracy pressure experienced by lighter stars closer to
the Chandrasekhar mass)111111The model described below is also known as the
Eddington $\beta$ model in literature.. Since the lower limit of the mass gap
is much higher than the Chandrasekhar mass for $m_{X}\gtrsim
0.1~{}\text{GeV}$, this assumption will hold for all cases studied here. The
pressure inside the star can be written as (note everything is in natural
units with $\hbar=c=k_{B}=1$)

$\displaystyle P$ $\displaystyle=P_{gas}+P_{rad}$ (B.1)
$\displaystyle=\frac{\rho}{\mu m_{X}}T+\frac{\pi^{2}}{45}T^{4}\,,$ (B.2)

where $\mu=0.5$ is the mean molecular weight in units of $m_{X}$. Assuming
$\beta=P_{gas}/P$ and $T=\xi\rho^{1/3}$ we get

$\displaystyle\xi=\left(\frac{45(1-\beta)}{\pi^{2}(\mu
m_{X})\beta}\right)^{1/3}\,.$ (B.3)

The assumption of $T=\xi\rho^{1/3}$ then gives us a $n=3$ polytrope with
$P=K\rho^{4/3}$ with

$K=\frac{\xi}{\beta(\mu
m_{X})}=\left(\frac{45}{\pi^{2}}\right)^{1/3}\left(\frac{(1-\beta)}{(\beta\mu
m_{X})^{4}}\right)^{1/3}\,.$

For an $n=3$ polytrope we can define change of variable into a dimensionless
parameter $z=Ar$ where

$A=\left(\frac{\pi G}{K}\right)^{1/2}\rho_{c}^{1/3},$

where $\rho_{c}$ is the central density of the star. The polytrope solution to
the density is given by the Lane-Emden equation of n = 3 given by

$\displaystyle\frac{d^{2}w}{dz^{2}}+\frac{2}{z}\frac{dw}{dz}+w^{3}=0\,.$ (B.4)

The solution has the first zero at $z_{3}=6.89$ which corresponds to the
radius of the star $R=\frac{z_{3}}{A}$. The density at a given radius,
parameterized by $w$, is given by $\rho=\rho_{c}w^{3}$. The mass of the star
$M_{*}$ can be written as

$\displaystyle M_{*}=\int_{0}^{R}4\pi
r^{2}\rho(r)dr=\int_{0}^{z_{3}}4\pi\frac{\rho_{c}}{A^{3}}w^{3}(z)z^{2}dz\simeq
1225.39\times\frac{M^{3}_{\rm pl}}{(\mu
m_{X})^{2}}\left(\frac{1-\beta}{\beta^{4}}\right)^{1/2}\,,$ (B.5)

where $M_{\rm pl}=2.435\times 10^{18}\ ~{}\text{GeV}$ is the reduced Planck
mass. Note that $\rho_{c}$ dependence disappears for $M_{*}$ for n = 3
polytrope. This uniquely determines $\beta$ in terms of $M_{*}$.

$\displaystyle\frac{\beta^{2}}{(1-\beta)^{1/2}}=15\times\left(\frac{~{}M_{\odot}}{M_{*}}\right)\left(\frac{~{}\text{GeV}}{\mu
m_{X}}\right)^{2}\,.$ (B.6)

###  B.1 Cooling timescale for $n=3$ Polytrope

In order to figure out if the star cools within the lifetime of the universe
we need to find out the surface temperature of the star. It’s true that
$T(R)=0$ by definition of a $n=3$ polytrope, however $R$ is not the surface of
the star. Here we will define the surface of the star $r_{0}$ as the radius
from which an emitted photon has an optical depth of $\frac{2}{3}$121212Using
the Eddington gray approximation the effective temperature is equal to the
temperature of the star at optical depth $\tau=\frac{2}{3}$. . Here we assume
Kramer’s opacity law for the gas at the surface parametrized by
$\kappa=\kappa_{0}\rho T^{-7/2}$. Note that Kramer’s opacity law assumes free-
free and bound-free scattering as dominant where $\kappa_{0}$ is determined
the dominant process. Here we ignore $H^{-}$ opacity and electron scattering
opacity because the surface temperatures we will deal with are higher than
$H^{-}$ binding energy and lower than the temperatures at which Compton
scattering dominates. We will determine if these assumptions are justified in
what follows. The scattering rate is then given by $\alpha_{K}=\kappa\rho$ and
the optical depth at radius r is given by

$\displaystyle\tau(r)=\int_{r}^{R}\alpha_{K}(r)dr=\int_{r}^{R}\kappa_{0}\rho(r)^{2}T(r)^{-7/2}dr=\int_{z}^{z_{3}}\kappa_{0}\frac{1}{A\xi^{7/2}}\rho_{c}^{5/6}w(z)^{5/2}dz\,.$
(B.7)

Let us assume $\kappa_{0}=f\times
0.38\times\frac{\alpha^{3}_{D}}{m_{c}^{3/2}m^{2}_{X}}$ where $f=1$ corresponds
to free free opacity. Let $z_{0}$ be the $z$ at which the optical depth is
$\frac{2}{3}$ which we define to be the surface of the star. We will assume
$\delta z_{0}=z_{3}-z_{0}\ll 1$ 131313This is done to get an analytical
approximation for surface and it will be justified in what follows.. We can
Taylor expand $w(z)=-w^{\prime}(z)|_{z=z_{3}}\delta z_{0}=0.042\,\delta
z_{0}$, then

$\displaystyle\tau=2.5\times 10^{-5}\times f\times\left(\frac{M_{\rm
pl}\alpha^{2}_{D}\sqrt{\mu\rho_{c}}}{\left(m_{c}m_{X}\right)^{3/2}}\right)\left(\frac{\sqrt{\beta}}{1-\beta}\right)\left(\delta
z_{0}\right)^{7/2}=\frac{2}{3}\,.$ (B.8)

This gives us

$\displaystyle\delta z_{0}=8.7\times
10^{-3}\left(\frac{100}{f}\right)^{2/7}\left(\frac{m_{c}}{m_{e}}\right)^{3/7}\left(\frac{m_{X}}{m_{p}}\right)^{3/7}\left(\frac{\alpha}{\alpha_{D}}\right)^{6/7}\left(\frac{1-\beta}{\beta^{1/2}}\right)^{2/7}\left(\frac{10^{3}\rm{gm~{}cm^{-3}}}{\rho_{c}}\right)^{1/7}\,.$
(B.9)

Note that the benchmark value $\left(\rho_{c}=10^{3}\rm{gm~{}cm^{-3}}\right)$
we use above is the density around at which core collapse begins for SM
parameters. Here we assume $f=100$ (which is the case in SM) i.e. the
Rosseland mean for the atomic opacity is assumed to be $10^{2}$ times larger
than free-free opacity 141414$f=100$ is a conservative estimate as larger $f$
implies lower temperature and hence a longer cooling time.. Then we get

$T_{s}=T(r_{0})=\xi\rho^{1/3}_{c}w(z_{0})=1.29\times 10^{-8}\
~{}\text{GeV}\left(\frac{100}{f}\right)^{2/7}\left(\frac{m_{c}}{m_{e}}\right)^{3/7}\\\
\times\left(\frac{m_{X}}{m_{p}}\right)^{2/21}\left(\frac{\alpha}{\alpha_{D}}\right)^{6/7}\left(\frac{(1-\beta)^{13}}{\beta^{10}}\right)^{1/21}\left(\frac{\rho_{c}}{10^{3}\rm{gm~{}cm^{-3}}}\right)^{4/21}\,.$
(B.10)

For SM parameters for a $M_{*}=60~{}M_{\odot}$ we get $\beta=0.75$. Using the
benchmark density and $f=100$, we get $T_{s}=72,000K$. This compares well to
the value from a MESA simulation which gives $10^{5}K$.The small difference
can be accounted for by the several approximations we made 151515The
photosphere of a star is also effectively a $n=3.25$ polytrope, not a $n=3$
polytrope as we assumed. We also assumed the opacity was defined at a hard
cutoff (Eq. B.8) of the optical depth to which our final answer is sensitive..
As shown in Fig. 4, our analytical approximation of the surface density
matches reasonably well in slope as well as magnitude when compared to the
exact results from MESA.

Figure 4: Comparison of temperature-density evolution between MESA (blue) and
equations B.3 and B.10 (red) for $M_{*}=100~{}M_{\odot}$. The dashed lines
indicate a comparison between the surface temperatures, and the full lines
represent the core temperatures.We have only shown the mesa evolution up to
$T_{c}=m_{c}/10$.

###  B.2 Equation of State and Dynamical Instability

The calculation of $\Gamma_{1}$ requires computation of the complete equation
of state. We note that the adiabatic index is, by definition, a local
quantity. This means that even if a particular region of the ”star” is
dynamically unstable, it is possible that the system as a whole is stable.
Ideally, one should use a global quantity to describe the system, such as the
pressure-weighted averaged adiabatic index (see e.g. [61])

$\Big{\langle}\Gamma_{1}\Big{\rangle}\equiv\frac{\int\Gamma_{1}PdV}{\int
PdV}.$ (B.11)

However, in this work, we consider the local adiabatic index, evaluated at the
center of the star for simplicity. We consider a star experiences pair
instability and collapse for any $M>M_{\rm crit}$, where $M_{\rm crit}$ is the
lowest mass star whose adiabatic index in the core drops below $4/3$.

Figure 5: Dynamic instability regions and polytropic evolution lines for
$m_{X}=1~{}\text{GeV}$.

The calculation of the ion and radiation contribution to the total pressure is
straightforward for a given density and temperature as follows. The
contribution to the total pressure and entropy from radiation is

$\displaystyle P_{\mathrm{rad}}$ $\displaystyle=\frac{a}{3}T^{4}$ (B.12)
$\displaystyle s_{\mathrm{rad}}$
$\displaystyle=\frac{4}{3}a\frac{T^{3}}{\rho}\,,$ (B.13)

while the contribution to the total pressure and entropy from ions is

$\displaystyle P_{\mathrm{ion}}$
$\displaystyle=\Big{\langle}\frac{1}{A}\Big{\rangle}\frac{\rho}{m_{X}}k_{B}T$
(B.14) $\displaystyle s_{\mathrm{ion}}$
$\displaystyle=\Big{\langle}\frac{1}{A}\Big{\rangle}\frac{k_{B}}{m_{X}}\Big{(}\frac{5}{2}+\log(\frac{T^{3/2}}{\rho})\Big{)}\,,$
(B.15)

where $\Big{\langle}\frac{1}{A}\Big{\rangle}$ is the mean reciprocal atomic
weight, which for ionized (dark) hydrogen is $\approx 1$. Note that the
electron contribution will be dealt with separately to allow for relativistic
and quantum effects. The contribution from dark electrons and positrons ,
$P_{c}$, requires integration over the density of states but is numerically
tractable. We outline this procedure here for completeness, following [6, 20].
We define the following quantities,

$C=\frac{1}{\pi^{2}}\Big{(}\frac{m_{c}c}{\hbar}\Big{)}^{3},\hskip
8.53581pt\beta=\frac{m_{c}c^{2}}{k_{B}T},\hskip
8.53581pt\phi=\frac{\mu}{k_{B}T}\,,$ (B.16)

where $\mu$ is the chemical potential of the negatively charged charge
carrier. The chemical potential is determined by imposing charge neutrality
and solving for the roots of the excess charge carrier number density,
satisfying

$n_{c}(\phi,\beta)=n_{c}^{-}-n_{c}^{+}=\Big{\langle}\frac{Z}{A}\Big{\rangle}\frac{\rho}{m_{X}}.$
(B.17)

With these definitions, the number density, pressure, and entropy are given by

$\displaystyle n_{c}(\phi,\beta)$
$\displaystyle=m_{c}CF_{2}^{+}(\phi,\beta)\,,$ (B.18) $\displaystyle
P_{c}(\phi,\beta)$ $\displaystyle=m_{c}c^{2}CF_{1}(\phi,\beta)\,,$ (B.19)
$\displaystyle s_{c}(\phi,\beta)$
$\displaystyle=\frac{k_{B}C\beta}{\rho}\Big{(}F_{1}(\phi,\beta)+F_{3}(\phi,\beta)-\frac{\phi}{\beta}F_{2}^{-}(\phi,\beta)\Big{)}\,,$
(B.20)

where the modified incomplete Fermi-Dirac integrals that appear in the
previous definitions are given by:

$\displaystyle F_{1}(\phi,\beta)$
$\displaystyle=\int_{\varepsilon=\beta}^{\infty}\Gamma\Big{(}\frac{\varepsilon}{\beta}\Big{)}D^{+}(\varepsilon,\phi)\frac{d\varepsilon}{\beta}\,,$
(B.21) $\displaystyle F_{2}^{+}(\phi,\beta)$
$\displaystyle=\int_{\varepsilon=\beta}^{\infty}\Gamma^{\prime}\Big{(}\frac{\varepsilon}{\beta}\Big{)}D^{+}(\varepsilon,\phi)\frac{d\varepsilon}{\beta}\,,$
(B.22) $\displaystyle F_{2}^{-}(\phi,\beta)$
$\displaystyle=\int_{\varepsilon=\beta}^{\infty}\Gamma^{\prime}\Big{(}\frac{\varepsilon}{\beta}\Big{)}D^{-}(\varepsilon,\phi)\frac{d\varepsilon}{\beta}\,,$
(B.23) $\displaystyle F_{3}(\phi,\beta)$
$\displaystyle=\int_{\varepsilon=\beta}^{\infty}\varepsilon\Gamma^{\prime}\Big{(}\frac{\varepsilon}{\beta}\Big{)}D^{+}(\varepsilon,\phi)\frac{d\varepsilon}{\beta^{2}}\,,$
(B.24)

where $\Gamma(x)\equiv\frac{1}{3}(x^{2}-1)^{3/2}$ and
$\varepsilon\equiv\frac{E}{k_{B}T}$. The adiabatic index (4.4) can be written
as

$\frac{\rho}{P}\Big{(}\frac{\partial
P}{\partial\rho}\Big{)}_{s}=\frac{\rho}{P}\Big{[}\Big{(}\frac{\partial
P}{\partial\rho}\Big{)}_{T}+\Big{(}\frac{\partial P}{\partial
T}\Big{)}_{\rho}\Big{(}\frac{\partial T}{\partial\rho}\Big{)}_{s}\Big{]}\,,$
(B.25)

which can then be calculated numerically term-by-term from the above
expressions.

Figure 6: Dynamic instability regions and polytropic evolution lines for
$m_{X}=10~{}\text{GeV}$.

Under what conditions will a dark star enter the dynamic instability region
and subsequently collapse? It turns out that, assuming an $n=3$ polytropic
evolution, whether or not a dark star will enter the region of dynamic
instability is independent of the charge carrier mass and dark fine-structure
constant, $m_{c},\alpha_{D}$ and only depends on the dark proton mass,
$m_{X}$. Although the region where $\Gamma_{1}<4/3$ varies with $m_{c}$, a
single line corresponding to $T\sim\rho^{1/3}$ is tangent to all regions. In
other words, even though the density and temperature at which a star would
enter the $\Gamma_{1}<4/3$ region depends on $m_{c}$, the ultimate fate of the
dark star is independent of the charge carrier mass, given sufficient time for
the system to evolve. In contrast, the final fate of a dark star is strongly
dependent on the dark proton mass, $m_{X}$. For larger $m_{X}$, lighter stars
will pass through the dynamic instability region throughout their evolution,
causing them to collapse. For $m_{X}=10\ ~{}\text{GeV}$, a star with $M\gtrsim
30~{}M_{\odot}$ will become dynamically unstable. Empirically, the critical
mass scales as

$M_{\mathrm{crit}}\approx
260~{}M_{\odot}\Big{(}\frac{\text{GeV}}{m_{X}}\Big{)}\,.$ (B.26)

Here, we assume that all dark stars with sufficiently high masses will
collapse.

## References

  * [1] LIGO Scientific Collaboration and Virgo Collaboration collaboration, B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley et al., _Observation of gravitational waves from a binary black hole merger_ , _Phys. Rev. Lett._ 116 (2016) 061102.
  * [2] LIGO Scientific, VIRGO, KAGRA collaboration, R. Abbott et al., _The population of merging compact binaries inferred using gravitational waves through GWTC-3_ , 2111.03634.
  * [3] S. E. Woosley and A. Heger, _The pair-instability mass gap for black holes_ , _The Astrophysical Journal Letters_ 912 (2021) L31.
  * [4] B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _New cosmological constraints on primordial black holes_ , _Phys. Rev. D_ 81 (2010) 104019 [0912.5297].
  * [5] S. Shandera, D. Jeong and H. S. G. Gebhardt, _Gravitational Waves from Binary Mergers of Subsolar Mass Dark Black Holes_ , _Phys. Rev. Lett._ 120 (2018) 241102 [1802.08206].
  * [6] G. Rakavy and G. Shaviv, _Instabilities in Highly Evolved Stellar Models_ , _apj_ 148 (1967) 803.
  * [7] Z. Barkat, G. Rakavy and N. Sack, _Dynamics of supernova explosion resulting from pair formation_ , _Phys. Rev. Lett._ 18 (1967) 379.
  * [8] S. E. Woosley, _Pulsational pair-instability supernovae_ , _The Astrophysical Journal_ 836 (2017) 244.
  * [9] A. Heger and S. E. Woosley, _The nucleosynthetic signature of population III_ , _The Astrophysical Journal_ 567 (2002) 532.
  * [10] LIGO Scientific, VIRGO collaboration, R. Abbott et al., _GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run_ , 2108.01045.
  * [11] LIGO Scientific, VIRGO, KAGRA collaboration, R. Abbott et al., _GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run_ , 2111.03606.
  * [12] D. Gerosa and E. Berti, _Are merging black holes born from stellar collapse or previous mergers?_ , _Phys. Rev. D_ 95 (2017) 124046 [1703.06223].
  * [13] M. Fishbach, D. E. Holz and B. Farr, _Are LIGO’s Black Holes Made From Smaller Black Holes?_ , _Astrophys. J. Lett._ 840 (2017) L24 [1703.06869].
  * [14] C. L. Rodriguez, M. Zevin, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio et al., _Black holes: The next generation—repeated mergers in dense star clusters and their gravitational-wave properties_ , _Phys. Rev. D_ 100 (2019) 043027 [1906.10260].
  * [15] LIGO Scientific, Virgo collaboration, R. Abbott et al., _Properties and Astrophysical Implications of the 150 M ⊙ Binary Black Hole Merger GW190521_, _Astrophys. J. Lett._ 900 (2020) L13 [2009.01190].
  * [16] L. A. C. van Son, S. E. de Mink, F. S. Broekgaarden, M. Renzo, S. Justham, E. Laplace et al., _Polluting the pair-instability mass gap for binary black holes through super-Eddington accretion in isolated binaries_ , _Astrophys. J._ 897 (2020) 100 [2004.05187].
  * [17] M. Spera, M. Mapelli, N. Giacobbo, A. A. Trani, A. Bressan and G. Costa, _Merging black hole binaries with the SEVN code_ , 1809.04605.
  * [18] U. N. Di Carlo, N. Giacobbo, M. Mapelli, M. Pasquato, M. Spera, L. Wang et al., _Merging black holes in young star clusters_ , _Mon. Not. Roy. Astron. Soc._ 487 (2019) 2947 [1901.00863].
  * [19] D. Croon, S. D. McDermott and J. Sakstein, _Missing in axion: Where are XENON1T’s big black holes?_ , _Phys. Dark Univ._ 32 (2021) 100801 [2007.00650].
  * [20] D. Croon, S. D. McDermott and J. Sakstein, _New physics and the black hole mass gap_ , _Phys. Rev. D_ 102 (2020) 115024 [2007.07889].
  * [21] J. Sakstein, D. Croon, S. D. McDermott, M. C. Straight and E. J. Baxter, _Beyond the Standard Model Explanations of GW190521_ , _Phys. Rev. Lett._ 125 (2020) 261105 [2009.01213].
  * [22] J. Ziegler and K. Freese, _Filling the black hole mass gap: Avoiding pair instability in massive stars through addition of nonnuclear energy_ , _Phys. Rev. D_ 104 (2021) 043015 [2010.00254].
  * [23] Z. Chacko, H.-S. Goh and R. Harnik, _The Twin Higgs: Natural electroweak breaking from mirror symmetry_ , _Phys. Rev. Lett._ 96 (2006) 231802 [hep-ph/0506256].
  * [24] N. Arkani-Hamed, T. Cohen, R. T. D’Agnolo, A. Hook, H. D. Kim and D. Pinner, _Solving the Hierarchy Problem at Reheating with a Large Number of Degrees of Freedom_ , _Phys. Rev. Lett._ 117 (2016) 251801 [1607.06821].
  * [25] H. Goldberg and L. J. Hall, _A New Candidate for Dark Matter_ , _Phys. Lett. B_ 174 (1986) 151.
  * [26] D. E. Kaplan, G. Z. Krnjaic, K. R. Rehermann and C. M. Wells, _Atomic Dark Matter_ , _JCAP_ 05 (2010) 021 [0909.0753].
  * [27] D. E. Kaplan, G. Z. Krnjaic, K. R. Rehermann and C. M. Wells, _Dark Atoms: Asymmetry and Direct Detection_ , _JCAP_ 10 (2011) 011 [1105.2073].
  * [28] F.-Y. Cyr-Racine and K. Sigurdson, _Cosmology of atomic dark matter_ , _Phys. Rev. D_ 87 (2013) 103515 [1209.5752].
  * [29] J. M. Cline, Z. Liu and W. Xue, _Millicharged Atomic Dark Matter_ , _Phys. Rev. D_ 85 (2012) 101302 [1201.4858].
  * [30] J. M. Cline, Z. Liu, G. Moore and W. Xue, _Scattering properties of dark atoms and molecules_ , _Phys. Rev. D_ 89 (2014) 043514 [1311.6468].
  * [31] J. Fan, A. Katz, L. Randall and M. Reece, _Double-Disk Dark Matter_ , _Phys. Dark Univ._ 2 (2013) 139 [1303.1521].
  * [32] J. Fan, A. Katz, L. Randall and M. Reece, _Dark-Disk Universe_ , _Phys. Rev. Lett._ 110 (2013) 211302 [1303.3271].
  * [33] J. Fan, A. Katz and J. Shelton, _Direct and indirect detection of dissipative dark matter_ , _JCAP_ 06 (2014) 059 [1312.1336].
  * [34] F.-Y. Cyr-Racine, R. de Putter, A. Raccanelli and K. Sigurdson, _Constraints on Large-Scale Dark Acoustic Oscillations from Cosmology_ , _Phys. Rev. D_ 89 (2014) 063517 [1310.3278].
  * [35] L. Randall and M. Reece, _Dark matter as a trigger for periodic comet impacts_ , _Physical Review Letters_ 112 (2014) .
  * [36] E. Rosenberg and J. Fan, _Cooling in a Dissipative Dark Sector_ , _Phys. Rev. D_ 96 (2017) 123001 [1705.10341].
  * [37] A. Ghalsasi and M. McQuinn, _Exploring the astrophysics of dark atoms_ , _Phys. Rev. D_ 97 (2018) 123018 [1712.04779].
  * [38] M. I. Gresham, H. K. Lou and K. M. Zurek, _Astrophysical Signatures of Asymmetric Dark Matter Bound States_ , _Phys. Rev. D_ 98 (2018) 096001 [1805.04512].
  * [39] R. Essig, S. D. Mcdermott, H.-B. Yu and Y.-M. Zhong, _Constraining Dissipative Dark Matter Self-Interactions_ , _Phys. Rev. Lett._ 123 (2019) 121102 [1809.01144].
  * [40] G. Alvarez and H.-B. Yu, _Astrophysical probes of inelastic dark matter with a light mediator_ , _Phys. Rev. D_ 101 (2020) 043002 [1911.11114].
  * [41] J.-S. Roux and J. M. Cline, _Constraining galactic structures of mirror dark matter_ , _Phys. Rev. D_ 102 (2020) 063518 [2001.11504].
  * [42] F.-Y. Cyr-Racine, F. Ge and L. Knox, _Symmetry of Cosmological Observables, a Mirror World Dark Sector, and the Hubble Constant_ , _Phys. Rev. Lett._ 128 (2022) 201301 [2107.13000].
  * [43] J. M. Cline, _Dark atoms and composite dark matter_ , 2108.10314.
  * [44] Z. Chacko, D. Curtin, M. Geller and Y. Tsai, _Direct detection of mirror matter in Twin Higgs models_ , _JHEP_ 11 (2021) 198 [2104.02074].
  * [45] M. Ryan, J. Gurian, S. Shandera and D. Jeong, _Molecular Chemistry for Dark Matter_ , 2106.13245.
  * [46] J. Gurian, D. Jeong, M. Ryan and S. Shandera, _Molecular Chemistry for Dark Matter II: Recombination, Molecule Formation, and Halo Mass Function in Atomic Dark Matter_ , 2110.11964.
  * [47] M. Ryan, S. Shandera, J. Gurian and D. Jeong, _Molecular Chemistry for Dark Matter III: DarkKROME_ , 2110.11971.
  * [48] A. Howe, J. Setford, D. Curtin and C. D. Matzner, _How to search for mirror stars with Gaia_ , _JHEP_ 07 (2022) 059 [2112.05766].
  * [49] N. Blinov, G. Krnjaic and S. W. Li, _Realistic model of dark atoms to resolve the Hubble tension_ , _Phys. Rev. D_ 105 (2022) 095005 [2108.11386].
  * [50] S. Bansal, J. H. Kim, C. Kolda, M. Low and Y. Tsai, _Mirror twin Higgs cosmology: constraints and a possible resolution to the H 0 and S8 tensions_, _JHEP_ 05 (2022) 050 [2110.04317].
  * [51] A. Cruz and M. McQuinn, _Astrophysical Plasma Instabilities induced by Long-Range Interacting Dark Matter_ , 2202.12464.
  * [52] G. Peled and T. Volansky, _Constraining Dark Matter Inside Stars Using Spectroscopic Binaries and a Modified Mass-Luminosity Relation_ , 2203.09522.
  * [53] V. Bromm, P. S. Coppi and R. B. Larson, _The Formation of the First Stars. I. The Primordial Star Forming Cloud_ , _The Astrophysical Journal_ 564 (2002) 23.
  * [54] K. Schutz, T. Lin, B. R. Safdi and C.-L. Wu, _Constraining a Thin Dark Matter Disk with Gaia_ , _Phys. Rev. Lett._ 121 (2018) 081101 [1711.03103].
  * [55] T. D. Brandt, _Constraints on MACHO Dark Matter from Compact Stellar Systems in Ultra-Faint Dwarf Galaxies_ , _Astrophys. J. Lett._ 824 (2016) L31 [1605.03665].
  * [56] MACHO collaboration, C. Alcock et al., _The MACHO project: Microlensing results from 5.7 years of LMC observations_ , _Astrophys. J._ 542 (2000) 281 [astro-ph/0001272].
  * [57] Ł. Wyrzykowski, A. E. Rynkiewicz, J. Skowron, S. Kozłowski, A. Udalski, M. K. Szymański et al., _OGLE-III Microlensing Events and the Structure of the Galactic Bulge_ , _ApJS_ 216 (2015) 12 [1405.3134].
  * [58] W. H. Press and P. Schechter, _Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation_ , _ApJ_ 187 (1974) 425.
  * [59] J. R. Bond, S. Cole, G. Efstathiou and N. Kaiser, _Excursion Set Mass Functions for Hierarchical Gaussian Fluctuations_ , _ApJ_ 379 (1991) 440.
  * [60] C. G. Lacey and S. Cole, _Merger rates in hierarchical models of galaxy formation_ , _Mon. Not. Roy. Astron. Soc._ 262 (1993) 627.
  * [61] M. Renzo, R. Farmer, S. Justham, Y. Götberg, S. E. de Mink, E. Zapartas et al., _Predictions for the hydrogen-free ejecta of pulsational pair-instability supernovae_ , _Astronomy & Astrophysics_ 640 (2020) A56.