# Measuring dark matter spikes around primordial black holes with Einstein
Telescope and Cosmic Explorer

Philippa S. Cole<EMAIL_ADDRESS>Gravitation Astroparticle Physics Amsterdam
(GRAPPA),
Institute for Theoretical Physics Amsterdam and Delta Institute for
Theoretical Physics,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Adam Coogan<EMAIL_ADDRESS>Gravitation Astroparticle Physics
Amsterdam (GRAPPA),
Institute for Theoretical Physics Amsterdam and Delta Institute for
Theoretical Physics,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
Département de Physique, Université de Montréal, 1375 Avenue Thérèse-Lavoie-
Roux, Montréal, QC H2V 0B3, Canada  Mila – Quebec AI Institute, 6666 St-
Urbain, #200, Montreal, QC, H2S 3H1  Bradley J. Kavanagh
<EMAIL_ADDRESS>Instituto de Física de Cantabria (IFCA, UC-CSIC), Av.
de Los Castros s/n, 39005 Santander, Spain  Gianfranco Bertone
<EMAIL_ADDRESS>Gravitation Astroparticle Physics Amsterdam (GRAPPA),
Institute for Theoretical Physics Amsterdam and Delta Institute for
Theoretical Physics,
University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

###### Abstract

Future ground-based gravitational wave observatories will be ideal probes of
the environments surrounding black holes with masses
$110\text{\,}\mathrm{\mathrm{M}_{\odot}}$. Binary black hole mergers with mass
ratios of order $q=m_{2}/m_{1}\lesssim 10^{-3}$ can remain in the frequency
band of such detectors for months or years, enabling precision searches for
modifications of their gravitational waveforms with respect to vacuum
inspirals. As a concrete example of an environmental effect, we consider here
a population of binary primordial black holes which are expected to be
embedded in dense cold dark matter spikes. We provide a viable formation
scenario for these systems compatible with all observational constraints, and
predict upper and lower limits on the merger rates of small mass ratio pairs.
Given a detected signal of one such system by either Einstein Telescope or
Cosmic Explorer, we show that the properties of the binary and of the dark
matter spike can be measured to excellent precision with one week’s worth of
data, if the effect of the dark matter spike on the waveform is taken into
account. However, we show that there is a risk of biased parameter inference
or missing the events entirely if the effect of the predicted dark matter
overdensity around these objects is not properly accounted for.

dark matter — gravitational waves — primordial black holes

##  I Introduction

Future ground-based gravitational wave observatories will be capable of
detecting binary black hole mergers at unprecedented distances Hall (2019).
Whilst event rates will be dominated by close-to-equal mass mergers, the
enhanced sensitivity and wider frequency range of the proposed Einstein
Telescope (ET) Punturo  _et al._ (2010); Maggiore  _et al._ (2020) and Cosmic
Explorer (CE) Evans  _et al._ (2021) observatories will also open the door to
observing intermediate mass ratio mergers.

The planned frequency ranges of Cosmic Explorer, $f\sim[5,5000]\,{\rm Hz}$,
and Einstein Telescope $f\sim[1,5000]\,{\rm Hz}$ mean that the observatories
will be sensitive to long duration signals from light systems with mass ratios
of order $q=m_{2}/m_{1}\sim 10^{-3}$, where $m_{1}$ is the mass of the central
compact object and $m_{2}$ is the mass of its lighter companion. Binaries with
$m_{1}=1-10\,{\rm M_{\odot}}$ will remain in band for many months or years.
Detecting such systems would have exciting consequences for gravitational wave
astronomy as well as fundamental physics Barack  _et al._ (2019); Bertone  _et
al._ (2020).

Systems with small mass-ratios are particular interesting because they are
more likely to be influenced by environmental effects. The gravitational
waveform of a binary inspiralling through an environment will be different to
that of the equivalent system merging in vacuum. In the case of an environment
of collisionless matter, dynamical friction, accretion, and a varying mass
enclosed within the orbit alter the dynamics of the binary Eda  _et al._
(2013, 2015); Macedo  _et al._ (2013); Barausse  _et al._ (2014, 2015); Yue
and Han (2018); Cardoso and Maselli (2020); Hannuksela  _et al._ (2020);
Kavanagh  _et al._ (2020); Coogan  _et al._ (2022). This appears as a gradual
change in the cumulative phase of the waveform with respect to the system in
vacuum, detectable if a very large number of cycles are observed. Dense
environments are more likely to survive around small mass-ratio systems,
unlike in equal-mass binaries, where any environments are likely to be
disrupted on a timescale of a few orbits Merritt  _et al._ (2002); Kavanagh
_et al._ (2018a). Given the amount of time these systems will spend in band,
this places Einstein Telescope and Cosmic Explorer well to detect this
dephasing effect in small mass-ratio binaries.

In this work, we explore how well ET and CE can measure the properties of
primordial black hole (PBH) binaries embedded in cold dark matter (DM) spikes.
PBHs may be formed shortly after the end of inflation from large density
fluctuations, contributing to the non-baryonic content of the Universe Green
and Kavanagh (2021). However, PBHs cannot make up all of the dark matter by
themselves in the stellar mass range (see e.g. Clesse and García-Bellido
(2017); Young and Byrnes (2020a) for caveats), and therefore must be
accompanied by another dark component. If the remaining DM is made up of cold,
collisionless particles, it will form dense spikes around primordial black
holes with a well-defined density profile Bertschinger (1985); Mack  _et al._
(2007); Ricotti (2007); Boudaud  _et al._ (2021), which will have an effect on
a PBH binary’s dynamics. Intermediate or extreme mass ratio PBH binaries can
therefore only be found and correctly interpreted if the effect of the dark
matter spike is taken into account. Moreover, electromagnetic signatures of DM
spikes around PBHs are largely ruled out for $f_{\mathrm{PBH}}\gtrsim 10^{-8}$
Lacki and Beacom (2010); Adamek  _et al._ (2019); Bertone  _et al._ (2019);
Carr  _et al._ (2021), making gravitational wave searches a particularly
promising avenue for discovery of such mixed DM scenarios.

We show that Einstein Telescope and Cosmic Explorer are ideally positioned to
observe gravitational wave signals from dark matter-dressed PBH binaries. In
particular, the range of GW frequencies accessible to both experiments makes
them sensitive to solar and sub-solar mass PBH binaries in the local Universe
(Section II). We provide a concrete formation scenario for such PBH binaries
and predict upper and lower limits on the merger rates of intermediate-mass-
ratio pairs (Section III). We specify the properties of the DM spikes which
are expected to form around PBHs and describe how the influence of these DM
spikes on the gravitational waveform can be modelled (Section IV). With these
tools, we show that it will be possible to distinguish these systems from GR-
in-vacuum inspirals and to measure the properties of the binary and the dark
matter spike (Sections V and VI). This is only possible if parameter
estimation is conducted with waveform templates that take into account the
effects of dynamical friction; we risk missing these signals in the data if it
is assumed that all binaries are inspiralling in vacuum. We conclude by
discussing these challenges associated with realistic search and inference
strategies, as well as relevance of our results to other environmental effects
around BH binaires (Section VII).

##  II Einstein Telescope and Cosmic Explorer reach

Firstly, we assess which systems are best-placed to be detected by various
observatories.

The improvement in sensitivity and frequency reach of Einstein Telescope and
Cosmic Explorer beyond aLIGO is shown by the noise power spectral density
curves $S_{n}(t)$ in Fig. 1. Throughout this paper, we choose our ‘benchmark’
PBH binary to have masses $m_{1}=1\,{\rm M_{\odot}}$ and $m_{2}=10^{-3}\,{\rm
M_{\odot}}$. The characteristic strain Moore  _et al._ (2015) of the
gravitational wave signal caused by the inspiralling of black holes with these
masses is shown by the rightmost black line in Fig. 1, for one week’s worth of
frequency evolution. This system lies very nicely in the frequency range of
the future ground-based detectors. The presence of a dense DM spike around the
heavier PBH will influence the precise frequency evolution of the system,
while inspiral of the binary will lead to a gradual depletion of the spike due
to feedback effects Kavanagh  _et al._ (2020).

The black dot on each trajectory in Fig. 1 marks the ‘break frequency’
$f_{b}$; at frequencies below $f_{b}$ the timescale for the depletion of the
DM spike is much shorter than the timescale for inspiral due to GW emission,
while at frequencies above $f_{b}$ the evolution of the DM spike becomes
negligible. This parametrization in terms of $f_{b}$ is relevant for
calculating the effects of feedback on the spike and will be discussed later
in Section IV. As argued in Ref. Coogan  _et al._ (2022), the important point
is that in order to detect the dephasing of the gravitational waveform, the
break frequency should lie in a region of high sensitivity of the detector. In
addition, the system should remain in band for a long time period so that as
many dephased cycles as possible are observed. For example, for this benchmark
system with a year of observations, approximately $10^{7}$ cycles occur in
band.

For comparison, we also show the LISA noise curve which lies at far lower
frequencies, making it more suited to observing the higher-mass systems which
have previously been studied in the context of DM-induced dephasing Eda  _et
al._ (2013, 2015); Yue and Han (2018); Yue  _et al._ (2019); Hannuksela  _et
al._ (2020); Cardoso and Maselli (2020); Kavanagh  _et al._ (2020); Dai  _et
al._ (2021); Li  _et al._ (2021); Coogan  _et al._ (2022); Becker  _et al._
(2022); Speeney  _et al._ (2022). We also plot the characteristic strain for
five years worth of frequency evolution of the system with $m_{1}=10^{3}\,{\rm
M_{\odot}}$, $m_{2}=1.4\,{\rm M_{\odot}}$ embedded in a dark matter spike,
which was explored in Coogan  _et al._ (2022). It was shown that the
properties of the dark matter spike could be reconstructed to very good
accuracy with five years of observations, and motivates this investigation
into what future ground-based observatories can do for dark matter spike
searches.

Figure 1: Characteristic strain as a function of frequency for the noise in
various detectors as well as two primordial dark dress systems. The heavy and
light systems “trajectories” (thick black lines) begin five years and one week
before the innermost stable circular orbit frequency. The systems’ distances
were respectively chosen to given signal-to-noise ratios of 15 at LISA and 12
at Cosmic Explorer over these observing durations. The vertical dashed lines
indicate the frequencies at a year, month, week and day before coalescence for
the light binary. The black dots indicate the frequency above which depletion
of the spike via dynamical friction from the compact object becomes negligible
(see Section IV). The system in the LISA band was studied in detail in Ref.
Coogan  _et al._ (2022).

For each detector, the distance at which our benchmark system will be
detectable varies. We calculate the SNR averaged over sky position,
orientation and polarization angle as a function of chirp mass
$\mathcal{M}=(m_{1}m_{2})^{\frac{3}{5}}/(m_{1}+m_{2})^{\frac{1}{5}}$ and
luminosity distance to source $d_{L}$ via

$\mathrm{SNR}=\frac{1}{d_{L}}\sqrt{4\int\frac{4}{5}\frac{h_{0}^{2}(f)}{S_{n}(f)}\mathrm{d}f}\,,$
(1)

where $S_{n}(f)$ is the power spectral density for a detector and the
amplitude of the Fourier transform of the gravitational waveform $h_{0}$ is
constructed from the following expressions for the phase $\Phi$ and its
derivatives with respect to time Maggiore (2007):

$\displaystyle\begin{split}h_{0}(f)&=\frac{1}{2}\frac{4\pi^{2/3}G^{5/3}\mathcal{M}^{5/3}f^{2/3}}{c^{4}}\sqrt{\frac{2\pi}{\ddot{\Phi}}}\,,\\\
\ddot{\Phi}&=4\pi^{2}f\quantity(\derivative{\Phi}{f})^{-1}\,,\\\
\Phi(f)&=\int_{f}^{f_{\mathrm{ISCO}}}\frac{\mathrm{d}t}{\mathrm{d}f^{\prime}}f^{\prime}\,\mathrm{d}f^{\prime}\,,\\\
t(f)&=\frac{1}{2\pi}\int^{f}\frac{\differential{f^{\prime}}}{f^{\prime}}\derivative{\Phi}{f}\,.\end{split}$
(2)

Here $G$ is Newton’s gravitational constant, $c$ is the speed of light, and
$f_{\mathrm{ISCO}}$ is the GW frequency at the innermost stable circular orbit
(ISCO), which we take to be the point of merger.

The SNRs for aLIGO, ET and CE are shown in Fig. 2 where we highlight an SNR
threshold of 12 as our definition of detectability, although note that this
threshold could be lowered if multiple detectors are online simultaneously.
With aLIGO, our benchmark system would only be detectable within a small
volume of $<300\,{\rm Mpc^{3}}$. However, for ET this increases to
$10^{6}\,{\rm Mpc^{3}}$ and for CE to $\sim 3\times 10^{7}\,{\rm Mpc^{3}}$.

We will now estimate the merger rates that we can expect for these light
systems, and use the distances calculated above for an SNR threshold of 12 to
predict the number of events we can hope to detect per year.

Figure 2: SNR s for benchmark system observed by terrestrial detectors. The
vertical red line indicates the chirp mass of our benchmark
$(m_{1},m_{2})=(1,10^{-3})$\mathrm{\mathrm{M}_{\odot}}$$ system. The
horizontal lines indicate the distances at which this system would have an
(optimal) SNR of 12 in each detector, with their colors matching the ones used
in Figure 1. The diagonal contours indicate curves of constant SNR, with the
red one highlighting the special value of 12 which we use to define a system
as detectable. Note that the SNR for vacuum binaries and the dark dresses we
consider is the same since it is not substantially altered by the DM spike.
Since a vacuum binary’s frequency at merger is not purely a function of its
chirp mass, in this plot we assume each system is observed for
$1\text{\,}\mathrm{week}$ up to the high-frequency cutoff for each detector.

##  III PBH merger rate

We now calculate how many PBH mergers with small mass ratios $q<10^{-2.5}$ we
can hope to detect with ET and CE. We first calculate the distribution of PBHs
that can be expected to form from a primordial power spectrum which is boosted
by $6-7$ orders of magnitude with respect to the observed amplitude on CMB
scales. We then calculate the merger rate for this population, and isolate
those with mass ratios where a dark matter spike is expected to have survived
around the larger of the two black holes.

We do not account for the effect of the DM spikes on the PBH merger rate, as
was done in for example Kavanagh  _et al._ (2020) for equal-mass mergers. In
principle, dynamical friction from the DM spike may reduce the merger time of
the binary, affecting the predicted merger rate for these systems. However,
detailed calculations have not yet been performed for the large mass ratio PBH
binaries we consider here (though see Ref. Pilipenko  _et al._ (2022) for
recent work on circumbinary accretion of DM). We also neglect the effects of
baryonic accretion on the PBH merger rate, which is relevant only for PBH
masses larger than a few solar masses De Luca  _et al._ (2020a).

###  III.1 PBH mass function

As a concrete example of a realisation of PBH formation, we assume that PBHs
can form from large overdensities which collapse shortly after single-field
inflation. We assume an initially Gaussian density field, and that there is no
clustering, however see Ünal  _et al._ (2021); Clesse and García-Bellido
(2017); Young and Byrnes (2020a); Ballesteros  _et al._ (2018); Bringmann  _et
al._ (2019); Atal  _et al._ (2020) for other possibilities. Note that PBHs can
also form from e.g. the collapse of cosmic strings Hawking (1989); Polnarev
and Zembowicz (1991) or during phase transitions in the early universe
Crawford and Schramm (1982); Hawking  _et al._ (1982); Kodama  _et al._
(1982).

A reference primordial power spectrum (PPS) that satisfies all current
constraints can be modelled by a piece-wise power law spectrum Atal and
Germani (2019); Byrnes  _et al._ (2019):

$\displaystyle\begin{split}\mathcal{P}_{\mathcal{R}}(k)&=A\begin{cases}\left(\frac{k}{k_{p}}\right)^{n_{g}}&k\leqslant
k_{p}\,,\\\
\left(\frac{k}{k_{p}}\right)^{-n_{d}}&k>k_{p}\,.\end{cases}\end{split}$ (3)

We choose $n_{g}=4$ because it is representative of peaks produced via ultra-
slow-roll models of inflation Atal and Germani (2019); Byrnes  _et al._
(2019), and we choose $n_{d}=0.05$ as being the slowest decay possible (in
order to enhance production of light PBHs) without conflicting with
observational constraints on scales smaller than the peak Gow  _et al._
(2021). We then compute the mass function of PBHs that would be formed, shown
in Fig. 3, using the Press-Schechter formalism as laid out in Ref. Gow  _et
al._ (2021) and accounting for the effect of the evolving equation of state in
the early universe Carr  _et al._ (2019); Byrnes  _et al._ (2018). This
results in a mass function with multiple peaks owing to the fact that PBHs
form more easily at times when the equation of state was lower. The most
prominent peak is around $1\,{\rm M_{\odot}}$, because scales with this
horizon mass were collapsing at the time of the QCD phase transition, where
the equation of state decreases by a factor of $\sim 1/3$ Byrnes  _et al._
(2018). Since we are interested in stellar-mass black holes, we choose
$k_{p}=10^{5}\,{\rm Mpc^{-1}}$ and verify that changes to $k_{p}$ within a
factor of 10 are always dominated by the QCD effect and still preferentially
produce solar-mass PBHs. This further motivates our choice of benchmark mass
$m_{1}=1\,{\rm M_{\odot}}$.

We do not take into account various considerations that may affect the
relationship between the amplitude of the primordial power spectrum and the
PBH abundance Kalaja  _et al._ (2019), for example non-Gaussianity, the shape
of the curvature perturbation Musco (2019) or the non-linear relationship
between the density and curvature perturbation Young  _et al._ (2019);
Kawasaki and Nakatsuka (2019); Luca  _et al._ (2019) because we are not trying
to constrain $f_{\rm PBH}$ with a specific amplitude of the power spectrum.
Instead, we fix the shape of the spectrum and then have the freedom to vary
$A$ in order to obey $f_{\rm PBH}$ constraints or to absorb subtleties in the
calculation of $f_{\rm PBH}$ from the initial distribution of densities.
However, we note that if the initial conditions are non-Gaussian, the initial
distribution of PBHs will be clustered, and this will effect the merger rate,
although it is uncertain by how much as the level of non-Gaussianity increases
Ballesteros  _et al._ (2018); Young and Byrnes (2020b).

We find that for $k_{p}=10^{5}\,\mathrm{Mpc^{-1}}$, an amplitude of
$A=3.5\times 10^{-3}$ saturates the constraints as described below in Section
III.3 and leads to $f_{\mathrm{PBH}}=0.066$.

Figure 3: The mass function of PBHs adopted in this work.  We normalize the
mass function according to $\int f(m)\,\mathrm{d}\ln m=f_{\mathrm{PBH}}$. The
multiple peak structure arises from the evolving equation of state in the
early universe Carr  _et al._ (2019). The peak of the primordial power
spectrum is at $k_{p}=10^{5}\,\rm{Mpc}^{-1}$, and the abundance of PBHs in
this case is $f_{\mathrm{PBH}}=0.066$.

###  III.2 Differential merger rate calculation

Two neighbouring PBHs may form a binary in the early universe when their self-
gravity dominates over the Hubble flow Nakamura  _et al._ (1997); Ioka  _et
al._ (1998); Ali-Haïmoud  _et al._ (2017). Two further requirements are that
the radial tidal forces from all other PBHs and matter fluctuations are weaker
than the attraction of the pair and tidal torques large enough to prevent a
head-on collision Raidal  _et al._ (2019); Liu  _et al._ (2019).

The differential merger rate for PBH binaries formed in the early universe
(early-forming binaries (EB) being dominant over those formed by tidal capture
in PBH clusters in the late Universe Bird  _et al._ (2016); Korol  _et al._
(2020); Franciolini  _et al._ (2022a), at least for small values of
$f_{\mathrm{PBH}}$) with component masses $m_{1}$ and $m_{2}$ is given by
Raidal  _et al._ (2019); Hütsi  _et al._ (2021); Phukon  _et al._ (2021):

$\displaystyle\begin{split}\mathrm{d}R_{EB}&=f_{\rm SUP}\frac{1.6\times
10^{6}}{\mathrm{Gpc}^{3}\,\mathrm{yr}}f_{\mathrm{PBH}}^{\frac{53}{37}}\eta^{-\frac{34}{37}}\left(\frac{m_{1}}{M_{\odot}}\right)^{-\frac{32}{37}}\\\
&\quad\times\left(\frac{t}{t_{0}}\right)^{-\frac{34}{37}}\psi(m_{1})\psi(m_{2})\,\mathrm{d}m_{1}\,\mathrm{d}m_{2}\,,\end{split}$
(4)

where $f_{\rm SUP}$ is the so-called suppression factor (details below),
$\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$, $t$ is the cosmic time of merger, $t_{0}$
is the cosmic time today, and $\psi(m)$ is the mass function defined by:

$\int m\psi(m)\,\mathrm{d}\ln
m=\int\frac{f(m)}{f_{\mathrm{PBH}}}\,\mathrm{d}\ln m=1,$ (5)

where the second equality demonstrates the relationship with the notation for
the mass function $f(m)$ as plotted in Fig. 3 for clarity. The differential
merger rate calculated for the mass function in Fig. 3 including the
suppression factor, is shown in Fig. 4.

Figure 4: Differential merger rate of binaries formed in the early universe,
with $m_{2}<10^{-2.5}m_{1}$, including suppression factor (and therefore an
underestimate).

However, it is very important to note that this merger rate is only reliable
for $f_{\rm PBH}\ll 1$, which we satisfy, because otherwise PBH clusters may
form that would perturb binaries Vaskonen and Veermäe (2020); De Luca  _et
al._ (2020b); Tkachev  _et al._ (2020). It is also only valid for relatively
narrow mass functions. Since we are using a very broad mass function, and are
especially reporting merger rates for low mass ratios because those are the
binaries where the DM spike will survive, caution must be used. The reason for
this is that the suppression factor assumes that a population of lower mass
black holes will cause third-body disruption to binaries. However, a large
population of very low mass black holes are unlikely to take part in this
process very effectively. It is not clear by how much the rate is
underestimated, and that study is beyond the scope of this work, so we instead
report upper and lower bounds based on using $f_{\rm SUP}$ as in Phukon  _et
al._ (2021), which approximates the formulations in Raidal  _et al._ (2019);
Hütsi  _et al._ (2021)

$f_{\rm SUP}=S_{1}\times S_{2}\,,$ (6)

with

$\displaystyle S_{1}$ $\displaystyle=1.42\left(\frac{\langle m_{\rm
PBH}^{2}\rangle/\langle m_{\rm
PBH}\rangle^{2}}{\overline{N}+C}+\frac{\sigma_{M}^{2}}{f_{\rm
PBH}^{2}}\right)^{-\frac{21}{74}}e^{-\overline{N}}\,,$ (7) $\displaystyle
S_{2}$ $\displaystyle\approx{\rm min}\left(1,9.6\times 10^{-3}f_{\rm
PBH}^{-0.65}e^{0.03\log^{2}(f_{\rm PBH})}\right)\,,$ (8)

where $\langle m_{\rm PBH}^{2}\rangle$ and $\langle m_{\rm PBH}\rangle^{2}$
are the variance and squared mean of the PBH mass respectively,
$\sigma_{M}^{2}=0.005$, and $\overline{N}$ (for which we use the full
expression rather than the upper and lower limits given in Phukon  _et al._
(2021)), and $C$ are given by:

$\overline{N}=\frac{m_{1}+m_{2}}{\langle m_{\rm PBH}\rangle}\frac{f_{\rm
PBH}}{f_{\rm PBH}+\sigma_{M}}\,,$ (9)
$\displaystyle\begin{split}C&=\frac{f_{\rm PBH}^{2}\langle m_{\rm
PBH}^{2}\rangle}{\sigma_{M}^{2}\langle m_{\rm PBH}\rangle^{2}}\times\\\
&\left[\left(\frac{\Gamma(\frac{29}{37})}{\sqrt{\pi}}U\left(\frac{21}{74},\frac{1}{2},\frac{5f_{\rm
PBH}^{2}}{6\sigma_{M}^{2}}\right)\right)^{-\frac{74}{21}}-1\right]^{-1}\,,\end{split}$
(10)

with $\Gamma$ the gamma function and $U$ the confluent hypergeometric
function. The expression for $\overline{N}$ determines the effect of nearby
PBHs on the merger rate and using the full expression over-suppresses the
merger rate. This gives us our lower bound in all of the merger rate plots,
whilst $f_{\rm SUP}=1$ (i.e. no effect), gives us our upper bound. Given that
we have produced the largest possible merger rate by maximising the amplitude
of the PPS, a detection in the observing runs of Einstein Telescope (ET) and
Cosmic Explorer (CE) would put tight constraints on $f_{\rm SUP}$, modulo
effects on the merger rate caused by the dark matter spikes. Non-detection
would constrain combinations of $f_{\rm SUP}$ and power spectrum amplitude
$A$, but would not break that degeneracy.

###  III.3 Saturating observational constraints

In order to determine the most optimistic merger rate in this scenario, we
saturate direct constraints on the PBH abundance. We do so by increasing the
amplitude of the PPS $A$ until we hit the strongest constraints. At the time
of writing, LIGO/Virgo stochastic gravitational wave background (SGWB)
spectrum constraints due to BBH mergers111Constraints from the direct
observations of individual mergers with LIGO/Virgo/Kagra are typically weaker
than those from the SGWB in the sub-solar mass range (see e.g. Gow  _et al._
(2020); De Luca  _et al._ (2020c); Hall  _et al._ (2020); García-Bellido  _et
al._ (2021); Wong  _et al._ (2021); De Luca  _et al._ (2021); Franciolini  _et
al._ (2022b); Chen  _et al._ (2022); Hütsi  _et al._ (2021)). and the envelope
of microlensing constraints from various probes on the sub-solar mass range
are the strongest on the relevant scales Green and Kavanagh (2021); Kavanagh
(2022).

We calculate the resulting stochastic gravitational wave signal due to all
mergers by inputting Eq. 4 into the expression from Phinney (2001); Zhu  _et
al._ (2011); Raidal  _et al._ (2019)

$\displaystyle\begin{split}\frac{\mathrm{d}\Omega_{\rm
GW}(f)}{\mathrm{d}m_{1}\,\mathrm{d}m_{2}\,\mathrm{d}z}&=\frac{f}{\rho_{c}c^{2}H(z)(1+z)}\frac{\mathrm{d}E_{GW}}{\mathrm{d}f}(f_{r})\\\
&\qquad\qquad\times\frac{\mathrm{d}R_{\rm
EB}(m_{1},m_{2})}{\mathrm{d}m_{1}\,\mathrm{d}m_{2}}\,,\end{split}$ (11)

with $H(z)$ the Hubble factor and $\rho_{c}$ the critical density. The
observed GW frequency $f$ is related to the source-frame frequency as
$f_{r}=(1+z)f$, and the GW spectrum $\mathrm{d}E_{\mathrm{GW}}/\mathrm{d}f$
from merging BHs is given in App. A.

We then integrate over all $m_{1}$ and $m_{2}$, dividing by 2 to avoid double-
counting, and we integrate between $z=0$ and $z_{\rm eq}$. We increase the
amplitude of the primordial power spectrum (which filters through to $f(m)$)
until the spectrum without suppression factor hits the LIGO/Virgo O5
sensitivity curve. This means that the merger rates we now report will be
viable by the time ET and CE come online in terms of being compatible with
other gravitational wave constraints that will be realised in the meantime.

The microlensing constraints do not depend directly on the merger rate, since
they search for distinctive variatons in stellar light curves due to
individual black holes passing in front of stars in the centre of the Milky
Way or in nearby galaxies Allsman  _et al._ (2001); Tisserand  _et al._
(2007); Niikura  _et al._ (2019a, b). Therefore, using our results for the
mass function $f(m)$ and the prescription in Ref. Carr  _et al._ (2017), we
check that the values of $f_{\rm PBH}$ corresponding to $k_{\rm
p}=10^{5}\,{\rm Mpc^{-1}}$ and $A=3.5\times 10^{-3}$ which saturates the SGWB
constraints is below the threshold allowed by the microlensing constraints
from figure 3 of Green and Kavanagh (2021). Indeed, $f_{\rm PBH}=0.066$ whilst
the constraints require $f_{\rm PBH}<0.78$.

###  III.4 Binned merger rates

The merger rates of a central black hole of mass $m_{1}$, within a bin of
width $\pm 0.5\,{\rm M_{\odot}}$, merging with all masses below
$10^{-2.5}\,m_{1}$ are shown in Fig. 5. All rates are averaged between
redshifts 0 and 0.1, and the upper and lower limits of the filled regions are
calculated without and with the suppression factor respectively. The maximum
merger rate we find is $R_{\rm EB}\sim 40\,{\rm Gpc^{-3}yr^{-1}}$ for
$m_{1}=1\,{\rm M_{\odot}}$ without suppression factor, and $R_{\rm EB}\sim
2\,{\rm Gpc^{-3}yr^{-1}}$ with the suppression factor which acts as a lower
bound.

Figure 5: PBH Merger rate for $m_{1}\pm\,0.5\,M_{\odot}$, including mergers
with all secondary masses below $m_{2}<10^{-2.5}\,m_{1}$. The lower value of
the filled region is calculated including the suppression factor of Eq. 6 (and
therefore an underestimate), whilst the upper value is calculated without. All
rates are averaged between redshift 0 and 0.1.

In order to report an event rate for a specific benchmark system, we can
combine the binned merger rate for $m_{1}=1\pm 0.5\,{\rm M_{\odot}}$ and
$m_{2}=10^{-3}\,m_{1}$, shown in figure Fig. 6, with the observable volumes of
each detector from Section II. We find that aLIGO, ET and CE can expect to
observe $2\text{\times}{10}^{-6}$, $8\text{\times}{10}^{-3}$ and 0.3 events
with an SNR of at least 12 per year respectively.

Figure 6: Merger rate of black hole binaries with
$q=10^{-3},\,10^{-2},\,10^{-1}$, without suppression factor. All averaged
between redshift 0 and 0.1. The merger rate for our benchmark system is given
by the left-most bin of the $q=10^{-3}$ case.

Furthermore, we would expect to see larger mass ratio mergers (where the spike
will have been disrupted) occurring in vacuum. For example, for mass ratios of
$q=10^{-1}$ and $q=10^{-2}$, the expected merger rates for our most optimistic
model are also shown in Fig. 6. LIGO/Virgo searches such as Refs. Nitz and
Wang (2021a, b) can probe these mass ranges but are not sensitive enough to
probe these merger rates yet. However, ET and CE will be sensitive enough to
search for these vacuum systems Barsanti  _et al._ (2021). If the merger rate
for this region of the parameter space was observed to be in line with the
rates presented in Fig. 6, this would be a strong hint that the more extreme-
mass ratio systems, with dark matter spikes, exist.

###  III.5 Spike survival

Finally, in order to ensure that the spike will survive the first few
encounters of the inspiral, we check that the radius of closest approach is
larger than the separation of the binary at the lowest value of the detectors’
frequency range, $f_{\rm low}$, for our benchmark masses. The radius of
closest approach is

$r_{\rm min}=a(1-\sqrt{1-j^{2}}),$ (12)

where $a$ is the semi-major axis of the binary and $j$ is the angular momentum
(related to the eccentricity by $j=\sqrt{1-e^{2}}$). We fix the time of
coalescence of the binary system to be $13\,{\rm Gyr}$,

$t_{\rm coa}=\frac{3a^{4}j^{7}c^{5}}{85m_{1}m_{2}(m_{1}+m_{2})G^{3}}=13\,{\rm
Gyr}$ (13)

and then select the most probable values of $a$ and $j$ that satisfy this
using equation (5) of Kavanagh  _et al._ (2018b) to compare with $r_{\rm
low}=(G(m_{1}+m_{2})/(f_{\rm low}^{2}\pi^{2}))^{\frac{1}{3}}$. We plot the
probability distribution in Fig. 7, highlighting $t_{\rm coa}=13\,{\rm Gyr}$
with the green line. For a relevant range of values of $a$ and $j$, we find
$r_{\rm min}/r_{\rm low}>100$. This means that the first encounter while the
orbit is still elliptical should not destroy the spike within the radius that
corresponds to the lower end of the frequency band of ET and CE, $f_{\rm
low}=1\,{\rm Hz}$.

Figure 7: Probability distribution of PBH binaries that decouple in the early
universe as a function of semi-major axis and angular momentum. The green line
represents a coalescence time of 13 Gyr.

Then, using the definition of the tidal radius of the companion object from
equation (4) of van den Bosch  _et al._ (2017) which is a refinement of the
Roche limit that does not treat the central black hole as a point mass, we
more stringently check that for a separation given by the radius of closest
approach $r_{\rm min}$, the tidal radius is still greater than $r_{\rm low}$.
This is indeed the case with $r_{\rm tidal}/r_{\rm low}>10$ for the most
probable range of semi-major axes.

##  IV Dephasing due to dark matter spike

Now that we have seen that there could be a number of detectable events per
year, we turn our focus to how we will need to go about detecting them.

We have seen in the previous section that PBHs cannot explain the entire DM
budget in the mass range relevant for ET and CE searches. However, if they
account for a sub-dominant fraction of the DM, they must be accompanied by a
component of particle DM making up the remainder. Both analytical calculations
and verifying simulations Bertschinger (1985); Mack  _et al._ (2007); Adamek
_et al._ (2019); Boudaud  _et al._ (2021) have shown that in this scenario,
(cold, collisionless) DM particles will form dense spikes around PBHs with a
distinct power law density profile

$\rho(r)=\rho_{\mathrm{sp}}\left(\frac{r_{\mathrm{sp}}}{r}\right)^{\gamma_{s}}\,,$
(14)

where $\gamma_{s}=9/4$, $\rho_{\rm sp}=\frac{1}{2}\rho_{\rm eq}$ and $r_{\rm
sp}=\left(2Gm_{1}t_{\rm eq}^{2}\right)^{1/3}$ Adamek  _et al._ (2019). Here,
$r$ is the radial distance from the central PBH of mass $m_{1}$, with
$\rho_{\rm eq}$ and $t_{\rm eq}$ being the background density and cosmic time
at matter-radiation equality. We neglect relativistic corrections to the
simple power-law which occur for radii $r\sim r_{\mathrm{isco}}$ Sadeghian
_et al._ (2013); Ferrer  _et al._ (2017); Speeney  _et al._ (2022). At these
small radii, GW emission (rather than dynamical friction) is expected to
dominate. A self-consistent description of dynamical friction and feedback in
the relativistic regime is not yet available, though see Ref. Speeney  _et
al._ (2022) for estimates of the impact of post-Newtonian corrections in these
systems.

It is useful to reparametrize the density profile in terms of the DM density
$\rho_{6}$ at a distance $r_{6}\equiv 10^{-6}\,\mathrm{pc}$ from the BH, as in
Ref. Coogan  _et al._ (2022):

$\rho(r)=\rho_{6}\left(\frac{r_{6}}{r}\right)^{\gamma_{s}}\,.$ (15)

For the purposes of parameter estimation, we treat $\rho_{6}$ as a free
parameter which controls the overall normalization of the spike density. Of
course, for a specific formation scenario, the benchmark value of $\rho_{6}$
will depend on $m_{1}$; equating Equation 14 and Equation 15, we see that for
a PBH mass $m_{1}$, the benchmark density normalization should be
$\rho_{6}=1.396\times
10^{13}\,(m_{1}/M_{\odot})^{3/4}\,M_{\odot}/\mathrm{pc}^{3}$.222This value of
$\rho_{6}$ differs slightly from that given for PBHs in Ref. Coogan  _et al._
(2022), as we here use slightly updated cosmological parameters. In this work,
we take $\rho_{\mathrm{eq}}=2605\,M_{\odot}/\mathrm{pc}^{3}$ and
$t_{\mathrm{eq}}=1.617\times 10^{12}\,\mathrm{s}$, calculated using CAMB Lewis
_et al._ (2000); Howlett  _et al._ (2012), assuming Planck-2018 cosmology
Aghanim  _et al._ (2020).

If the dominant DM component is made up of canonical WIMP-like DM, then the
dense DM spikes around PBHs would give rise to a large gamma-ray flux due to
WIMP annihilation. Constraints from point source searches and the diffuse
gamma-ray flux therefore severely limit the possibility of mixed PBH-WIMP DM
scenarios Lacki and Beacom (2010); Adamek  _et al._ (2019); Bertone  _et al._
(2019). However, if the remaining DM is composed of another candidate, such as
axion-like particles or asymmetric WIMPs, these constraints can be evaded and
these dense DM spikes may have an impact on GW searches for PBHs Kavanagh  _et
al._ (2018).

Gravitational wave searches in the LIGO/Virgo data have predominantly used
matched filtering, where a bank of gravitational waveform templates which are
calculated for a particular set of parameters are compared to the data to
search for a ‘match’ with a pre-specified signal-to-noise threshold. Thus far,
these gravitational waveforms have been modelled assuming that the inspiral
and merger have occurred in vacuum. However, the gravitational waveform looks
different if the inspiral instead occurs in non-empty space, namely a DM
spike. We must, therefore, use non-vacuum waveforms to avoid missing signals
in the data due to SNR loss, and to avoid mischaracterizing signals, where the
use of vacuum waveforms may lead to biased parameter reconstruction.

As the smaller BH moves through the spike of DM particles, they will form a
wake which imparts a drag force on the BH, reducing its orbital velocity. This
in turn causes the smaller BH to drop into a lower orbit more quickly than it
would in vacuum. This effect is known as dynamical friction Chandrasekhar
(1943a, b, c). The inspiral happens in fewer GW cycles, causing the
gravitational waveform to gradually go out of phase with respect to the vacuum
case, an effect known as ‘dephasing’. Full details of the state-of-the-art in
calculating the DM dephasing effect are given in Ref. Kavanagh  _et al._
(2020); Coogan  _et al._ (2022); here we briefly summarize the physics and our
numerical approach.

Working at Newtonian order, the evolution of the binary separation $r$ for
quasi-circular orbits can be described by:

$\begin{split}&\dot{r}=-\frac{64\,G^{3}\,M\,m_{1}\,m_{2}}{5\,c^{5}\,r^{3}}\\\
&-\frac{8\pi\,G^{1/2}\,m_{2}\,r^{5/2}\,\rho_{\mathrm{DM}}(r,t)\,\xi(r,t)\,\log\Lambda}{\sqrt{M}m_{1}}\,,\end{split}$
(16)

where $M=m_{1}+m_{2}$ is the total mass of the binary. The first term on the
right hand side corresponds to GW emission, while the second corresponds to
dynamical friction.333We neglect contributions from accretion onto the smaller
BH and the varying enclosed mass due to the DM spike; in the absence of
feedback, these effects are expected to be sub-dominant Macedo  _et al._
(2013); Cardoso and Maselli (2020). The dynamical friction force traces the
density profile of the DM within the spike $\rho_{\mathrm{DM}}(r,t)$, while
the factor $\xi(r,t)$ corresponds to the fraction of DM particles at a given
radius $r$ moving more slowly that the local orbital speed, which are those
relevant for the calculation of the dynamical friction force (Binney and
Tremaine, 2008, Sec. 8.1). The Coulomb logarithm $\log\Lambda$ incorporates
information about the range of distances from the smaller BH at which
gravitational scattering with DM particles is effective. Following Ref.
Kavanagh  _et al._ (2020), we take $\Lambda=\sqrt{m_{1}/m_{2}}$.

During the inspiral, the motion of the BH binary will inject energy into the
DM spike, altering its density profile. In Ref. Kavanagh  _et al._ (2020), the
DM spike was described in terms of a spherically symmetric, isotropic
distribution function, which is altered by the gravitational scattering of DM
particles with the orbiting BH. Following the evolution of this distribution
allows us to calculate the time evolution of the DM density and velocity
distribution, $\rho_{\mathrm{DM}}(r,t)\xi(r,t)$, while simultaneously solving
for the binary separation $r$ through Equation 16. Given a trajectory $r(t)$
the GW frequency is calculated as:

$f(t)=\frac{1}{\pi}\sqrt{\frac{GM}{r(t)^{3}}}\,,$ (17)

from which the strain and phase to coalescence can be calculated using Eq. 2.

The feedback formalism described above was implemented in the publicly-
available HaloFeedback code Kavanagh (2020). The most prominent effect is that
the spike will be locally depleted as some particles are ejected, reducing the
size of the dephasing effect. Taking into account this effect requires time-
consuming numerical simulations using HaloFeedback. However, for the parameter
estimation that we would like to conduct, we need to produce waveforms for a
densely sampled region of the parameter space, and therefore we will rely on
the analytic model that was proposed and used in Coogan  _et al._ (2022) to
describe the dephasing including feedback. In this model, the phase to
coalescence $\Phi(f)$ approximately follows a broken power-law, with the break
frequency $f_{b}$ a function of $m_{1},m_{2}$ and $\gamma_{s}$. The functional
form for $f_{b}(m_{1},m_{2},\gamma_{s})$ was fit to HaloFeedback simulations
of binaries with total masses in the range $M\sim 10^{3}-10^{4}\,M_{\odot}$.
However, as we show in Appendix B, this parametrization also provides an
accurate fit to the behaviour of the much lighter systems we consider here. We
therefore rely on this parametrization to calculate the waveforms of light,
dressed PBH binaries in the following sections.

##  V Assessing detectability, discoverability and measurability

In order to assess the prospects for concretely measuring DM spikes around
PBHs with ground-based GW observatories, we follow the approach of Ref. Coogan
_et al._ (2022). Assuming Gaussian noise, the likelihood can be written as

$p\quantity(d|h_{\bm{\theta}})\propto\exp\quantity[\innerproduct{h_{\bm{\theta}}}{d}^{2}-\frac{1}{2}\innerproduct{h_{\bm{\theta}}}{h_{\bm{\theta}}}]\,,$
(18)

where $d(t)$ is the strain time series data measured by the detector (which we
assume to match the our benchmark signal $s(t)$) and $h_{\bm{\theta}}$ is a
model waveform $h_{\bm{\theta}}(t)$ with parameters $\bm{\theta}$. In Equation
18, the noise-weighted inner product is defined as:

$\innerproduct{a}{b}=4\real\int_{0}^{\infty}\differential{f}\frac{\tilde{a}(f)^{*}\,\tilde{b}(f)}{S_{n}(f)}\,,$
(19)

where $S_{n}(f)$ is the noise power spectral density (effectively the
sensitivity curve of a given detector as a function of GW frequency). For
Einstein Telescope, we assume the “ET-D” configuration Hild  _et al._ (2011),
making use of sensitivity estimates provided by the ET
collaboration444http://www.et-gw.eu/index.php/etsensitivities. For Cosmic
Explorer, we assume the sensitivity as given in Ref. Evans  _et al._ (2016).
In all cases, we adopt the sensitivity averaged over sky locations,
polarizations and binary orientations.

The parameters
$\bm{\theta}=\bm{\theta}_{\mathrm{int}}\cup\bm{\theta}_{\mathrm{ext}}$ which
describe the model waveform $h_{\bm{\theta}}$ can be divided into intrinsic
parameters $\bm{\theta}_{\mathrm{int}}$, which describe the properties of the
source, and extrinsic parameters $\bm{\theta}_{\mathrm{ext}}$, which depend on
the observer. The intrinsic parameters describing the systems are:

$\displaystyle\text{Vacuum:}\quad\bm{\theta}_{V,\mathrm{int}}$
$\displaystyle=\\{\mathcal{M}\\}$ (20) $\displaystyle\text{Dark
Dress:}\quad\bm{\theta}_{D,\mathrm{int}}$
$\displaystyle=\\{\gamma_{s},\rho_{6},\mathcal{M},\log_{10}q\\}\,,$ (21)

where $q=m_{2}/m_{1}$ is the mass ratio of the binary.555We define the masses
in the detector frame $m_{\mathrm{det}}$ related to the source-frame mass
$m_{\mathrm{src}}$ by $m_{\mathrm{det}}=m_{\mathrm{src}}(1+z)$. Note that at a
luminosity distance of $d_{\mathrm{L}}\sim 200\,\mathrm{Mpc}$ (roughly the
detectability horizon for these systems with ET and CE, as shown in Fig. 2),
the redshift is $z\sim 0.05$, so we expect only a small correction
distinguishing between detector- and source-frame masses. Assuming the
detector’s response is constant over the duration of the waveform, the
extrinsic parameters $\bm{\theta}_{\mathrm{ext}}$ we consider are the
luminosity distance to the system and the phase and time at coalescence:

$\bm{\theta}_{\mathrm{ext}}\equiv\quantity{d_{L},\phi_{c},\tilde{t}_{c}}\,.$
(22)

In practice, at each point in parameter space, we maximise the likelihood over
the extrinsic parameters using a fast Fourier transform, as described in Ref.
Coogan  _et al._ (2022). This maximised likelihood is denoted
$p_{\mathrm{max}}(d|h_{\bm{\theta}})$.

To compute the noise-weighted inner product, we integrate over the frequency
range beginning a week before the system being observed coalesces and ending
at the ISCO frequency. The use of this fixed frequency window introduces a
subtle issue: the waveform being compared with the observation may last for
more or less than a week over this frequency window. This means that it is not
quite correct to maximize the noise-weighted inner product using a Fourier
transform. Fixing this issue requires either numerically maximizing over
$\tilde{t}_{c}$ or including it in the parameter estimation. We postpone
further investigation to future work.

We explore the posterior distribution $p(\bm{\theta}|d)\propto
p_{\mathrm{max}}(d|h_{\bm{\theta}})p(\bm{\theta})$ using nested sampling
Skilling (2004); Higson  _et al._ (2019), implemented in the code dynesty
Speagle (2020). The priors we take for the intrinsic parameters are summarized
in Table 1. The initial slope of the density profile around PBHs has been
analytically predicted and numerically verified (e.g. Adamek  _et al._ (2019))
as $\gamma_{s}=2.25$, hence the reasonably narrow prior on this parameter,
given that for these masses, we must be observing primordial black holes. We
allow for a wide prior on the density normalisation which includes the vacuum
value, and a range of mass ratios that have an upper bound of $10^{-2.5}$
where the assumptions on the survival of the spike break down. Mapping out the
posteriors allows us to assess the question of measurability: how well the
properties of binary with a DM dress can be measured or constrained.

Parameter | Prior range
---|---
$\gamma_{s}$ | $[2.0,2.5]$
$\rho_{6}$ [$10^{15}$ $\mathrm{\mathrm{M}_{\odot}}\mathrm{/}\mathrm{pc}^{3}$] | $[0,2]$
$\mathcal{M}$ [$\mathrm{\mathrm{M}_{\odot}}$] | –
$\log_{10}q$ | $[-4.5,-2.5]$
Table 1: Prior ranges on intrinsic parameters. In each case, we assume a
uniform prior over the range given above. Since we use the same wide prior on
the chirp mass for the dark dress and vacuum systems, its precise value
cancels out in the Bayes factor and does not matter.

In order to address the question of discoverability (i.e. whether a given dark
dress waveform can be distinguished from a GR-in-vacuum waveform), we compute
the Bayes factor for the dark dress and vacuum models, defined as the evidence
ratio:

$\operatorname{BF}(d)\equiv\frac{p\quantity(d|\mathrm{D})}{p\quantity(d|\mathrm{V})}\,.$
(23)

The evidence for a given model is defined as:

$p\quantity(d)=\int\differential{\bm{\theta}}p_{\mathrm{max}}\quantity(d|h_{\bm{\theta}})\,p(\bm{\theta})\,,$
(24)

where $\bm{\theta}=\bm{\theta}_{D}$ or $\bm{\theta}=\bm{\theta}_{V}$ for the
dark dress and vacuum models respectively. Large values of the Bayes Factor
($\operatorname{BF}>100$) correspond to strong evidence in favor of a dark
dress system, compared to a GR-in-vacuum system Jeffreys (1998); Kass and
Raftery (1995).

Figure 8: Dephasing and SNR loss between best-fit vacuum template and dark
dresses with various masses. The rows correspond to ET and CE. Assumes one
week of observing with the system at a distance such that it is detected with
an SNR of 12 by each detector. Figure 9: One- and two-dimensional marginal
posteriors for the intrinsic parameters of PBH dark dresses observed with
Einstein Telescope for one week (orange), one month (green) or one year (red).
The dashed black lines indicate the true parameter values. The black hole
masses are set to the benchmark values of
$1\text{\,}\mathrm{\mathrm{M}_{\odot}}$ and
${10}^{-3}\text{\,}\mathrm{\mathrm{M}_{\odot}}$, and $\Delta\mathcal{M}$ is
the difference between the measured and true chirp mass
($\sim$1.585\text{\times}{10}^{-2}\text{\,}\mathrm{\mathrm{M}_{\odot}}$$). The
distances for the systems are set so each is detected with an SNR of 12. The
shaded contours show the $65\text{\,}\mathrm{\char 37\relax}$,
$95\text{\,}\mathrm{\char 37\relax}$ and $99.7\text{\,}\mathrm{\char
37\relax}$ credible regions. Figure 10: Marginal posteriors for a PBH dark
dress as in Figure 9, but instead observed with Cosmic Explorer. The system’s
distance was set to give an SNR of 12 at CE. Figure 11: The 1D marginal
posteriors for $\rho_{6}$ and $\gamma_{s}$ obtained using Einstein Telescope
and Cosmic Explorer. Note the target dark dress systems in the two rows are
_not_ the same. Their distances differ since each is observed with an SNR of
12 in its respective detector.

##  VI Results

Parameter | Einstein Telescope | Cosmic Explorer
---|---|---
$\gamma_{s}$ | $2.22_{-0.05}^{+0.07}$ | $2.20_{-0.06}^{+0.07}$
$\rho_{6}$ [$10^{13}$ $\mathrm{\mathrm{M}_{\odot}}\mathrm{/}\mathrm{pc}^{3}$] | $1.7_{-0.6}^{+0.8}$ | $2.0_{-0.8}^{+1}$
$\mathcal{M}$ [$\mathrm{\mathrm{M}_{\odot}}$] | $$1\,584\,577$_{-3}^{+3}\times 10^{-8}$ | $$1\,584\,578$_{-3}^{+3}\times 10^{-8}$
$\log_{10}(q)$ | $-3.1_{-0.3}^{+0.4}$ | $-3.2_{-0.3}^{0.4}$
Table 2: Parameters inferred for benchmark system with one week of
observations by different detectors. The error bars indicate the 68% credible
intervals.

#### Detectability.

We find that aLIGO, ET and CE will be able to detect our benchmark system with
a SNR of 12 out to a distance of 6.5, 78 and $286\text{\,}\mathrm{Mpc}$
respectively (see Fig. 2). These distances together with the upper bound on
the (binned) merger rate for systems with $m_{1}=$1\pm
0.5\text{\,}\mathrm{\mathrm{M}_{\odot}}$$ and $m_{2}=10^{-3}\,m_{1}$
corresponds to event rates of $2\text{\times}{10}^{-6}$,
$8\text{\times}{10}^{-3}$ and 0.3 per year. Increasing the mass of the systems
would increase the distance out to which ET and CE will be able to detect the
system, but the price to pay is that the merger rate decreases with larger
$m_{1}$ as shown in Fig. 6.

#### Discoverability.

We have evaluated the Bayes factor for a PBH binary with
$(m_{1},m_{2})=(10,10^{-2})\,$\text{\,}\mathrm{\mathrm{M}_{\odot}}$$, assuming
1 week of observation with ET and CE. We find $\mathrm{BF}\approx 10^{18}$ for
ET and $\mathrm{BF}\approx$6\text{\times}{10}^{6}$$ for CE, indicating
overwhelming evidence in favor of the presence of a dark dress in the system
in both cases. For this system we also find that the dephasing with respect to
the best-fit vacuum system is $\mathcal{O}(100)$ cycles. As shown in the left
panels of Fig. 8, lighter systems (including our benchmark system) than this
will have an even larger dephasing (for example, exceeding 1000 cycles for
$(m_{1},m_{2})=(10,10^{-3})\,$\text{\,}\mathrm{\mathrm{M}_{\odot}}$$). We
therefore conclude that over the parameter space of interest for light PBHs,
the presence of a dark dress should be discoverable with observation times of
1 week or more.

We also show in Fig. 8 the percentage SNR loss between the dephased system and
the best-fitting vacuum system. For the system with
$(m_{1},m_{2})=(10,10^{-2})\,$\text{\,}\mathrm{\mathrm{M}_{\odot}}$$ described
above, the SNR loss is roughly 40%. This indicates that searching for dark
dresses with vacuum templates will result in a large number of missed
detections relative to a search based on dark dress templates. Given that the
optimal matched-filtered SNR scales as $d_{L}^{-1}$ (see Eq. 1), this SNR loss
would correspond to a reduction of the detectable volume for these systems
$\sim d_{L}^{3}$ by a factor of $(0.6)^{3}\sim 4.6$. The use of vacuum
waveforms would substantially reduce the observable rate of large mass ratio
PBH mergers. Looking again at lighter systems, the SNR loss increases further.
For binaries with
$(m_{1},m_{2})=(1,10^{-3})\,$\text{\,}\mathrm{\mathrm{M}_{\odot}}$$, the SNR
loss with CE approaches $80\%$, reducing the observable volume by a factor
$\sim 100$, highlighting the importance of using dephased waveforms to
effectively search for such systems.

#### Measurability.

The one- and two-dimensional marginal posteriors for the intrinsic parameters
of a dressed PBH system with benchmark masses of
$m_{1}=$1\text{\,}\mathrm{\mathrm{M}_{\odot}}$$ and
$m_{2}=${10}^{-3}\text{\,}\mathrm{\mathrm{M}_{\odot}}$$ observed with Einstein
Telescope are shown in Fig. 9. The true values of the parameters are shown by
the dashed black lines, and the posteriors for week, month and year long
signals are shown by the orange, green and red contours respectively. With
even one week’s worth of data, it will be possible to measure the slope of the
density profile to precision $\gamma_{s}=2.22_{-0.05}^{+0.07}$,666 We use
error bars indicating the 68% credible interval.  and with a year’s worth of
data, $\gamma_{s}=2.25_{-0.01}^{0.01}$. The measured values for the other
parameters are given in Table 2.

We run the same parameter estimation for Cosmic Explorer, and find comparable
errors on the parameters to Einstein Telescope in Figure 10. We note, however,
that increasing the duration of the signal from one month to one year from
merger does not improve the size of the contours for Cosmic Explorer because
of the low-frequency cut-off of the noise curve - see Figure 1.

In Figure 11 we highlight the one-dimensional marginal posteriors for
$\rho_{6}$ and $\gamma_{s}$ to show the precision with which we can measure
these parameters, as well as the improvement from increasing the duration of
the signal. We stress, however, that one week’s worth of data will be enough
to measure these parameters with a few percent accuracy, and presents a less
daunting data analysis challenge than year-long signals.

So far, we have included all of the intrinsic parameters of the system in the
parameter estimation pipeline, specifically
$\\{\gamma_{s},\rho_{6},\mathcal{M},\log_{10}q\\}$. However, for PBH systems,
the value of $\gamma_{s}$ is fixed and the density normalization $\rho_{6}$ is
uniquely determined by the central PBH mass, as described in Section IV.
Assuming a PBH origin for the dephasing signal therefore allows us to reduce
the parameter estimation problem to two dimensions –
$\\{\rho_{6},\mathcal{M}\\}$ – which substantially strengthens constraints, as
we describe in Appendix C.

##  VII Discussion

We have explored the prospects for detecting the presence of Dark Matter (DM)
spikes around light primordial black hole (PBH) binaries with future ground-
based gravitational wave (GW) observatories, such as Einstein Telescope (ET)
and Cosmic Explorer (CE). We focus on a broad, well-motivated PBH mass
function with $f_{\mathrm{PBH}}\approx 0.066$, which is still allowed by
current constraints from microlensing and the stochastic GW background.

We can expect between 1 and 100 PBH mergers per $\mathrm{Gpc}^{3}$ volume per
year, and an event rate for our specific benchmark system of $m_{1}=1\pm
0.5\,{\mathrm{M}}_{\odot}$, $m_{2}=10^{-3}\,m_{1}$, of
$8\text{\times}{10}^{-3}$ and 0.3 events per year, observable with an SNR of
at least 12 in each of ET and CE respectively. We note that if ET and CE are
online simultaneously, the SNR threshold for each individual detector could be
lowered, and hence the distance at which PBH mergers of a given mass are
observable increases along with the event rate.

We find that for our benchmark system, the SNR lost if the waveform is matched
with a vacuum template could be as large as $80\%$. This suggests that search
strategies will need to take into account the effect of the DM spikes so as to
avoid missing these signals. For larger $m_{1}$ than our benchmark system, the
SNR loss will be smaller, so the system may be detectable with vacuum
templates but would lead to biased parameter inference.

We show that using the correct model for parameter estimation, i.e. including
the effects of the dark matter spike on the waveform, we can reconstruct the
intrinsic parameters of the binary, the chirp mass and the mass ratio, as well
as the parameters of the spike, the density normalisation and the power law of
the dark matter density profile, to very good precision with one week’s worth
of data, as summarised in Table 2.

We find that Einstein Telescope can measure the dark matter spike parameters
with better precision than Cosmic Explorer, owing to the lower frequency
reach, which allows more cycles to be observed and hence more dephasing to
accumulate. We note that there may be opportunities for multi-band
observations of systems that overlap with the frequency ranges of both LISA
and ET or CE Cutler  _et al._ (2019).

Since PBH binaries of these masses and mass ratios must be embedded in DM
spikes (assuming they cannot make up all of the dark matter themselves), we
conclude that in order to find these systems and measure their parameters
correctly, the effect of dark matter on the phase of the inspiral must be
taken into account in at least the parameter estimation process to avoid
biased parameter inference. For some ranges of the parameter space (for
example our benchmark system), the effect of the dark matter must be taken
into account in order to detect the signal in the first place, as the SNR loss
incurred by assuming the system is inspiralling in vacuum could be
catastrophic.

We also emphasise that less extreme mass ratio mergers will also exist for our
PBH formation scenario, which should be inspiralling in vacuum, since we
expect their dark matter spikes to have been disrupted. Observations of these
systems (potentially even with aLIGO/Virgo Nitz and Wang (2021a, b)) would be
a very strong indicator that the more extreme mass ratio systems are out there
to be found, and would provide strong motivation for conducting a search for
these exotic waveforms in the data.

These conclusions are drawn assuming that the search and inference will be
conducted using matched filtering with template banks. However, considering
the duration of the signals expected, techniques from continuous wave searches
Aasi (2014) may be more suitable for searching for these signals that go
through millions of cycles. Recently, Ref. Guo and Miller (2022) proposed the
use of the Hough Transform to search for ‘mini-EMRIs’, systems similar to
those we consider here. The authors estimate that for a strain sensitivity
similar to that of LIGO, a binary with masses $(10,\,10^{-2})\,M_{\odot}$ may
be detectable out to a distance of a few Mpc with this technique. This is
roughly a factor of two less than the detectable distance we estimate in Fig.
2, suggesting that the application of more realistic search strategies need
not substantially degrade the detectability of the signals.

In any case, it will be vital to understand the evolution of the frequency as
a function of time for these systems and how it differs from the vacuum case.
There is also potential for the use of machine learning to search the data for
such long duration signals, in order to decrease the expense of computing the
(many thousands of) waveforms of these systems directly.

While we have focused here on DM spikes around PBH binaries, many of the tools
and conclusions apply also to other environmental effects. These include
binaries embedded in accretion disks Graham  _et al._ (2020) or ‘gravitational
atoms’ (clouds of light scalar fields bound to BHs) Baumann  _et al._ (2019,
2022a, 2022b). For these systems too, future ground-based observatories will
provide exquisite sensitivity to slowly-accumulating dephasing effects. Our
results suggest that $\mathcal{O}(\mathrm{weeks})$ of data would be required
to extract useful physical information from such systems, though a detailed
study of this – and of whether the environmental effects from these different
sources can be distinguished – we leave for future work.

## Acknowledgements

We thank Thomas Edwards, Andrew Gow, Samaya Nissanke and Ville Vaskonen for
useful conversations.

P.C. acknowledges funding from the Institute of Physics, University of
Amsterdam. A.C. received funding from the Netherlands eScience Center (grant
number ETEC.2019.018) and the Schmidt Futures Foundation. B.J.K. thanks the
Spanish Agencia Estatal de Investigación (AEI, Ministerio de Ciencia,
Innovación y Universidades) for the support to the Unidad de Excelencia María
de Maeztu Instituto de Física de Cantabria, ref. MDM-2017-0765. We acknowledge
Santander Supercomputing support group at the University of Cantabria who
provided access to the supercomputer Altamira at the Institute of Physics of
Cantabria (IFCA-CSIC), member of the Spanish Supercomputing Network, for
performing simulations/analyses.

##  Appendix A GW spectrum from coalescing BHs

For the interested reader, we here give explicitly expressions for the GW
spectrum of coalescing black holes, which appears in the calculation of the
stochastic GW background in Section III.3.

The energy emitted by coalescing BHs with masses $m_{1}$ and $m_{2}$ in the GW
frequency range $[f,f+\mathrm{d}f]$ is given by Raidal  _et al._ (2017);
Chernoff and Finn (1993); Zhu  _et al._ (2011):

$\displaystyle\frac{\mathrm{d}E_{GW}}{\mathrm{d}f}=\mathcal{C}\begin{cases}f^{-\frac{1}{3}}&\text{for
}f<f_{1}\,,\\\ \frac{f^{\frac{2}{3}}}{f_{1}}&\text{for }f_{1}\leq
f<f_{2}\,,\\\
\frac{ff_{3}^{4}}{(f_{1}f_{2}^{\frac{4}{3}}(4(f-f_{2})^{2}+f_{4})^{2})^{2}}&\text{for
}f_{2}\leq f<f_{3}\,,\\\ 0&\text{for }f_{3}\leq f<f_{4}\,,\end{cases}$

where the prefactor is $\mathcal{C}=(\pi G)^{\frac{2}{3}}M^{\frac{5}{3}}\eta$.
The frequency window limits are given numerically by Ajith  _et al._ (2008):

$f_{j}=\frac{(a_{j}\eta^{2}+b_{j}\eta+c_{j})c^{3}}{\pi G(m_{1}+m_{2})}$ (25)

where e.g. $a_{j}$ denotes the $j$-th index of the following arrays:

$\displaystyle a$ $\displaystyle=(0.2971,0.59411,0.84845,0.50801),$ (26)
$\displaystyle b$ $\displaystyle=(0.04481,0.089794,0.12848,0.077515),$ (27)
$\displaystyle c$ $\displaystyle=(0.095560,0.19111,0.27299,0.022369),$ (28)

and $\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$.

##  Appendix B HaloFeedback Validation

Here, we present some validation tests which were performed for the $f(t)$
parametrization presented in Ref. Coogan  _et al._ (2022) and used in Sec. V
of this work.

We begin by checking that the formalism developed in Ref. Kavanagh  _et al._
(2020) behaves well for the low BH masses we consider here. In Figure 12, we
show the effective density profile extracted from simultaneously evolving the
separation of the PBH binary and the DM spike distribution function (as
summarized in Section IV and implemented in the HaloFeedback code Kavanagh
(2020)). The effective density is the instantaneous DM density experienced by
the smaller inspiraling BH when it reaches an orbital separation $r$. For
reference, the orbital reference corresponding to the break frequency $f_{b}$
for this system is shown as a vertical dotted line.

The qualitative behaviour of the effective density matches that observed in
the heavier systems presented in Ref. Coogan  _et al._ (2022). Because of
energy injected by the BH into the spike, the spike is rapidly depleted and
thus the effective density is smaller than the initial, unperturbed density
(grey dashed line). As the GW inspiral continues (from large $r$ to small
$r$), the timescale for depletion of the spike eventually becomes longer than
the timescale for the GW inspiral, and the effective density converges to the
initial unperturbed profile. Lines of different colors in Figure 12 correspond
to simulations which were started at different initial BH separations. After
an initial depletion phase, each of these density profiles converges to the
same behaviour, indicating that we do not need to explicitly specify the
initial separation of the binary at formation (as long as this is larger than
the radius at which the binary enters the GW observing band).

Figure 12: Effective density profile of the DM spike, obtained using
HaloFeedback. Each colored line corresponds to a simulation beginning at a
different initial PBH binary separation, while the diagonal dashed line shows
the unperturbed DM density profile. The behaviour of the system at small radii
is independent of the initial separation.

Having verified that the behaviour of HaloFeedback is sensible for these light
systems, we can now compare these simulations against the phase
parametrization $\Phi(f)$ which we use for parameter estimation in this work.
From the trajectories obtained in these simulations, we can calculate the
evolution of the GW phase with frequency and compare this with the
corresponding predictions from the analytic phase parametrization. In Figure
13, we plot the dephasing with respect to the vacuum case for PBH binaries
with masses $(m_{1},\,m_{2})=(1.0,\,10^{-3})\,M_{\odot}$. The full
HaloFeedback simulations are shown in solid blue, using a maximum timestep of
100 orbits (corresponding to a phase given by the horizonal dotted line). We
use this simulation to extract the effective DM density profile and resimulate
assuming a static spike with this density profile. This allows us to
substantially reduce the maximum timestep and therefore resolve the dephasing
much closer to the merger (dashed blue line). These results are matched
closely by the output of the analytical parametrization, as implemented in
pydd Coogan (2021), down to the level of a few cycles.

Despite being calibrated on BH masses a few orders of magnitude larger than
those being considered here, the analytic parametrization maintains percent-
level accuracy over many years of the inspiral. We have also verified this
behaviour for various mass ratios and total masses of the light PBH systems.

Figure 13: DM-induced dephasing with respect to the vacuum for different
modeling approaches. The solid blue line shows the dephasing modeled using the
HaloFeedback code Kavanagh (2020), solving simultaneously the inspiral and
feedback on the DM spike. For computational ease, the maximum step size is set
to 100 orbits (corresponding to a phase given by the horizontal dotted line).
The dashed blue line is obtained assuming a static DM spike with effective
density profile extracted from the full HaloFeedback run, allowing for a
shorter timestep. This matches closely the output of pydd Coogan (2021) (solid
orange line), which implements the analytic fit to the phase $\Phi(f)$ used in
Sec. V of this work and first presented in Ref. Coogan  _et al._ (2022).

##  Appendix C Posteriors for fixed $\gamma_{s}$

Figure 14: One- and two-dimensional marginal posteriors for a subset of the
intrinsic parameters of PBH dark dresses, keeping $\gamma_{s}$ and $q$ fixed,
observed with Einstein Telescope for times ranging from one week to one year.
The dashed black lines indicate the true parameter values. The black hole
masses are set to the benchmark values of
$1\text{\,}\mathrm{\mathrm{M}_{\odot}}$ and
${10}^{-3}\text{\,}\mathrm{\mathrm{M}_{\odot}}$, and $\Delta\mathcal{M}$ is
the difference between the measured and true chirp mass
($\sim$1.585\text{\times}{10}^{-2}\text{\,}\mathrm{\mathrm{M}_{\odot}}$$). The
distances for the systems are set so each is detected with an SNR of 12. The
contours indicate the $65\text{\,}\mathrm{\char 37\relax}$,
$95\text{\,}\mathrm{\char 37\relax}$ and $99.7\text{\,}\mathrm{\char
37\relax}$ credible regions. Figure 15: Marginal posteriors for a PBH dark
dress with $\gamma_{s}$ and $q$ fixed as in Figure 14, but instead observed
with Cosmic Explorer.

In the main body of the paper, we have performed parameter estimation for all
four intrinsic parameters of the dark dress systems: $\gamma_{s}$, $\rho_{6}$,
$\mathcal{M}$ and $q$. However, it is also interesting in the case of PBHs to
make use of the fact that there is a concrete prediction for the slope of the
spike $\gamma_{s}=9/4$, as well as that $\rho_{6}$ is directly related to the
mass of the primary black hole $m_{1}$. With this in mind, we fix
$\gamma_{s}=9/4$ and perform inference over the parameters $\rho_{6}$ and
$\mathcal{M}$, using Equation 14 and Equation 15 to fix the mass ratio $q$.
This is a strong assumption on two fronts: firstly it assumes that the system
is definitely primordial (reasonable if the primary mass is below
$1.4\text{\,}\mathrm{\mathrm{M}_{\odot}}$) and that the predictions from
analytic calculations and simulations are correct and have no scatter, and
secondly it assumes that the spike has not been at all disrupted down to
redshifts less than one, such that the density profile slope and relationship
between the normalisation and $m_{1}$ remains intact. While we do not foresee
this type of parameter estimation being implemented for the discovery of the
first system of this type, it demonstrates how well we can predict the
parameters of the system if we are confident that we are observing one of a
(perhaps previously confirmed) population of primordial black holes with pre-
existing evidence that the density profile slope of these systems is
$\gamma_{s}=9/4$.

The resulting posteriors are shown in Figure 14 for an analysis with Einstein
Telescope and in Figure 15 for an analysis with Cosmic Explorer.
Unsurprisingly, in both cases the posteriors are extremely narrow even with
just one day of observations before coalescence. The chirp mass can be
measured to within approximately $3\times 10^{-9}$ of the true value, and the
density normalisation $\rho_{6}$ to better than $1\text{\,}\mathrm{\char
37\relax}$ precision with $99.7\text{\,}\mathrm{\char 37\relax}$ confidence.

## References

  * Hall (2019) E. Hall, “Horizon plots and noise curves for second- and third-generation detectors,” LIGO Document T1800084-v5, https://dcc.ligo.org/LIGO-T1800084-v5/public (2019).
  * Punturo  _et al._ (2010) M. Punturo _et al._ , “The Einstein Telescope: A third-generation gravitational wave observatory,” Class. Quant. Grav. 27, 194002 (2010).
  * Maggiore  _et al._ (2020) Michele Maggiore _et al._ , “Science Case for the Einstein Telescope,” JCAP 03, 050 (2020), arXiv:1912.02622 [astro-ph.CO] .
  * Evans  _et al._ (2021) Matthew Evans _et al._ , “A Horizon Study for Cosmic Explorer: Science, Observatories, and Community,” (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Barack  _et al._ (2019) Leor Barack _et al._ , “Black holes, gravitational waves and fundamental physics: a roadmap,” Class. Quant. Grav. 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Bertone  _et al._ (2020) Gianfranco Bertone _et al._ , “Gravitational wave probes of dark matter: challenges and opportunities,” SciPost Phys. Core 3, 007 (2020), arXiv:1907.10610 [astro-ph.CO] .
  * Eda  _et al._ (2013) Kazunari Eda, Yousuke Itoh, Sachiko Kuroyanagi, and Joseph Silk, “New Probe of Dark-Matter Properties: Gravitational Waves from an Intermediate-Mass Black Hole Embedded in a Dark-Matter Minispike,” Phys. Rev. Lett. 110, 221101 (2013), arXiv:1301.5971 [gr-qc] .
  * Eda  _et al._ (2015) Kazunari Eda, Yousuke Itoh, Sachiko Kuroyanagi, and Joseph Silk, “Gravitational waves as a probe of dark matter minispikes,” Phys. Rev. D 91, 044045 (2015), arXiv:1408.3534 [gr-qc] .
  * Macedo  _et al._ (2013) Caio F. B. Macedo, Paolo Pani, Vitor Cardoso, and Luís C. B. Crispino, “Into the lair: gravitational-wave signatures of dark matter,” Astrophys. J. 774, 48 (2013), arXiv:1302.2646 [gr-qc] .
  * Barausse  _et al._ (2014) Enrico Barausse, Vitor Cardoso, and Paolo Pani, “Can environmental effects spoil precision gravitational-wave astrophysics?” Phys. Rev. D 89, 104059 (2014), arXiv:1404.7149 [gr-qc] .
  * Barausse  _et al._ (2015) Enrico Barausse, Vitor Cardoso, and Paolo Pani, “Environmental Effects for Gravitational-wave Astrophysics,” J. Phys. Conf. Ser. 610, 012044 (2015), arXiv:1404.7140 [astro-ph.CO] .
  * Yue and Han (2018) Xiao-Jun Yue and Wen-Biao Han, “Gravitational waves with dark matter minispikes: the combined effect,” Phys. Rev. D 97, 064003 (2018), arXiv:1711.09706 [gr-qc] .
  * Cardoso and Maselli (2020) Vitor Cardoso and Andrea Maselli, “Constraints on the astrophysical environment of binaries with gravitational-wave observations,” Astron. Astrophys. 644, A147 (2020), arXiv:1909.05870 [astro-ph.HE] .
  * Hannuksela  _et al._ (2020) Otto A. Hannuksela, Kenny C. Y. Ng, and Tjonnie G. F. Li, “Extreme dark matter tests with extreme mass ratio inspirals,” Phys. Rev. D 102, 103022 (2020), arXiv:1906.11845 [astro-ph.CO] .
  * Kavanagh  _et al._ (2020) Bradley J. Kavanagh, David A. Nichols, Gianfranco Bertone, and Daniele Gaggero, “Detecting dark matter around black holes with gravitational waves: Effects of dark-matter dynamics on the gravitational waveform,” Phys. Rev. D 102, 083006 (2020), arXiv:2002.12811 [gr-qc] .
  * Coogan  _et al._ (2022) Adam Coogan, Gianfranco Bertone, Daniele Gaggero, Bradley J. Kavanagh, and David A. Nichols, “Measuring the dark matter environments of black hole binaries with gravitational waves,” Phys. Rev. D 105, 043009 (2022), arXiv:2108.04154 [gr-qc] .
  * Merritt  _et al._ (2002) David Merritt, Milos Milosavljevic, Licia Verde, and Raul Jimenez, “Dark matter spikes and annihilation radiation from the galactic center,” Phys. Rev. Lett. 88, 191301 (2002), arXiv:astro-ph/0201376 .
  * Kavanagh  _et al._ (2018a) Bradley J. Kavanagh, Daniele Gaggero, and Gianfranco Bertone, “Merger rate of a subdominant population of primordial black holes,” Phys. Rev. D 98, 023536 (2018a), arXiv:1805.09034 [astro-ph.CO] .
  * Green and Kavanagh (2021) Anne M. Green and Bradley J. Kavanagh, “Primordial Black Holes as a dark matter candidate,” J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Clesse and García-Bellido (2017) Sébastien Clesse and Juan García-Bellido, “The clustering of massive Primordial Black Holes as Dark Matter: Measuring their mass distribution with advanced LIGO,” Physics of the Dark Universe 15, 142–147 (2017), arXiv:1603.05234 [astro-ph.CO] .
  * Young and Byrnes (2020a) Sam Young and Christian T. Byrnes, “Initial clustering and the primordial black hole merger rate,” Journal of Cosmology and Astroparticle Physics 2020, 004–004 (2020a).
  * Bertschinger (1985) E. Bertschinger, “Self - similar secondary infall and accretion in an Einstein-de Sitter universe,” Astrophys. J. Suppl. 58, 39 (1985).
  * Mack  _et al._ (2007) Katherine J. Mack, Jeremiah P. Ostriker, and Massimo Ricotti, “Growth of structure seeded by primordial black holes,” Astrophys. J. 665, 1277–1287 (2007), arXiv:astro-ph/0608642 .
  * Ricotti (2007) Massimo Ricotti, “Bondi accretion in the early universe,” Astrophys. J. 662, 53–61 (2007), arXiv:0706.0864 [astro-ph] .
  * Boudaud  _et al._ (2021) Mathieu Boudaud, Thomas Lacroix, Martin Stref, Julien Lavalle, and Pierre Salati, “In-depth analysis of the clustering of dark matter particles around primordial black holes. Part I. Density profiles,” JCAP 08, 053 (2021), arXiv:2106.07480 [astro-ph.CO] .
  * Lacki and Beacom (2010) Brian C. Lacki and John F. Beacom, “Primordial Black Holes as Dark Matter: Almost All or Almost Nothing,” Astrophys. J. Lett. 720, L67–L71 (2010), arXiv:1003.3466 [astro-ph.CO] .
  * Adamek  _et al._ (2019) Julian Adamek, Christian T. Byrnes, Mateja Gosenca, and Shaun Hotchkiss, “WIMPs and stellar-mass primordial black holes are incompatible,” Phys. Rev. D 100, 023506 (2019), arXiv:1901.08528 [astro-ph.CO] .
  * Bertone  _et al._ (2019) Gianfranco Bertone, Adam M. Coogan, Daniele Gaggero, Bradley J. Kavanagh, and Christoph Weniger, “Primordial Black Holes as Silver Bullets for New Physics at the Weak Scale,” Phys. Rev. D 100, 123013 (2019), arXiv:1905.01238 [hep-ph] .
  * Carr  _et al._ (2021) Bernard Carr, Florian Kuhnel, and Luca Visinelli, “Black holes and WIMPs: all or nothing or something else,” Mon. Not. Roy. Astron. Soc. 506, 3648–3661 (2021), arXiv:2011.01930 [astro-ph.CO] .
  * Moore  _et al._ (2015) C. J. Moore, R. H. Cole, and C. P. L. Berry, “Gravitational-wave sensitivity curves,” Class. Quant. Grav. 32, 015014 (2015), arXiv:1408.0740 [gr-qc] .
  * Yue  _et al._ (2019) Xiao-Jun Yue, Wen-Biao Han, and Xian Chen, “Dark matter: an efficient catalyst for intermediate-mass-ratio-inspiral events,” Astrophys. J. 874, 34 (2019), arXiv:1802.03739 [gr-qc] .
  * Dai  _et al._ (2021) Ning Dai, Yungui Gong, Tong Jiang, and Dicong Liang, “Intermediate Mass-Ratio Inspirals with Dark Matter Minispike,” (2021), arXiv:2111.13514 [gr-qc] .
  * Li  _et al._ (2021) Gen-Liang Li, Yong Tang, and Yue-Liang Wu, “Probing Dark Matter Spikes via Gravitational Waves of Extreme Mass Ratio Inspirals,” (2021), arXiv:2112.14041 [astro-ph.CO] .
  * Becker  _et al._ (2022) Niklas Becker, Laura Sagunski, Lukas Prinz, and Saeed Rastgoo, “Circularization versus eccentrification in intermediate mass ratio inspirals inside dark matter spikes,” Phys. Rev. D 105, 063029 (2022), arXiv:2112.09586 [gr-qc] .
  * Speeney  _et al._ (2022) Nicholas Speeney, Andrea Antonelli, Vishal Baibhav, and Emanuele Berti, “The impact of relativistic corrections on the detectability of dark-matter spikes with gravitational waves,” (2022), arXiv:2204.12508 [gr-qc] .
  * Maggiore (2007) Michele Maggiore,  _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford Master Series in Physics (Oxford University Press, 2007).
  * Pilipenko  _et al._ (2022) Sergey Pilipenko, Maxim Tkachev, and Pavel Ivanov, “Evolution of a primordial binary black hole due to interaction with cold dark matter and the formation rate of gravitational wave events,” Phys. Rev. D 105, 123504 (2022), arXiv:2205.10792 [astro-ph.CO] .
  * De Luca  _et al._ (2020a) V. De Luca, G. Franciolini, P. Pani, and A. Riotto, “Constraints on Primordial Black Holes: the Importance of Accretion,” Phys. Rev. D 102, 043505 (2020a), arXiv:2003.12589 [astro-ph.CO] .
  * Ünal  _et al._ (2021) Caner Ünal, Ely D. Kovetz, and Subodh P. Patil, “Multimessenger probes of inflationary fluctuations and primordial black holes,” Physical Review D 103 (2021), 10.1103/physrevd.103.063519.
  * Clesse and García-Bellido (2017) Sébastien Clesse and Juan García-Bellido, “The clustering of massive primordial black holes as dark matter: Measuring their mass distribution with advanced ligo,” Physics of the Dark Universe 15, 142–147 (2017).
  * Ballesteros  _et al._ (2018) Guillermo Ballesteros, Pasquale D. Serpico, and Marco Taoso, “On the merger rate of primordial black holes: effects of nearest neighbours distribution and clustering,” JCAP 10, 043 (2018), arXiv:1807.02084 [astro-ph.CO] .
  * Bringmann  _et al._ (2019) Torsten Bringmann, Paul Frederik Depta, Valerie Domcke, and Kai Schmidt-Hoberg, “Towards closing the window of primordial black holes as dark matter: The case of large clustering,” Phys. Rev. D 99, 063532 (2019), arXiv:1808.05910 [astro-ph.CO] .
  * Atal  _et al._ (2020) Vicente Atal, Albert Sanglas, and Nikolaos Triantafyllou, “LIGO/Virgo black holes and dark matter: The effect of spatial clustering,” JCAP 11, 036 (2020), arXiv:2007.07212 [astro-ph.CO] .
  * Hawking (1989) S.W. Hawking, “Black Holes From Cosmic Strings,” Phys. Lett. B 231, 237–239 (1989).
  * Polnarev and Zembowicz (1991) Alexander Polnarev and Robert Zembowicz, “Formation of Primordial Black Holes by Cosmic Strings,” Phys. Rev. D 43, 1106–1109 (1991).
  * Crawford and Schramm (1982) Matt Crawford and David N. Schramm, “Spontaneous Generation of Density Perturbations in the Early Universe,” Nature 298, 538–540 (1982).
  * Hawking  _et al._ (1982) S.W. Hawking, I.G. Moss, and J.M. Stewart, “Bubble Collisions in the Very Early Universe,” Phys. Rev. D 26, 2681 (1982).
  * Kodama  _et al._ (1982) Hideo Kodama, Misao Sasaki, and Katsuhiko Sato, “Abundance of Primordial Holes Produced by Cosmological First Order Phase Transition,” Prog. Theor. Phys. 68, 1979 (1982).
  * Atal and Germani (2019) Vicente Atal and Cristiano Germani, “The role of non-gaussianities in Primordial Black Hole formation,” Phys. Dark Univ. 24, 100275 (2019), arXiv:1811.07857 [astro-ph.CO] .
  * Byrnes  _et al._ (2019) Christian T. Byrnes, Philippa S. Cole, and Subodh P. Patil, “Steepest growth of the power spectrum and primordial black holes,” JCAP 06, 028 (2019), arXiv:1811.11158 [astro-ph.CO] .
  * Gow  _et al._ (2021) Andrew D. Gow, Christian T. Byrnes, Philippa S. Cole, and Sam Young, “The power spectrum on small scales: robust constraints and comparing PBH methodologies,” JCAP 2021, 002 (2021), arXiv:2008.03289 [astro-ph.CO] .
  * Carr  _et al._ (2019) Bernard Carr, Sebastien Clesse, Juan Garcia-Bellido, and Florian Kuhnel, “Cosmic Conundra Explained by Thermal History and Primordial Black Holes,” arXiv e-prints , arXiv:1906.08217 (2019), arXiv:1906.08217 [astro-ph.CO] .
  * Byrnes  _et al._ (2018) Christian T. Byrnes, Mark Hindmarsh, Sam Young, and Michael R. S. Hawkins, “Primordial black holes with an accurate QCD equation of state,” JCAP 2018, 041 (2018), arXiv:1801.06138 [astro-ph.CO] .
  * Kalaja  _et al._ (2019) Alba Kalaja, Nicola Bellomo, Nicola Bartolo, Daniele Bertacca, Sabino Matarrese, Ilia Musco, Alvise Raccanelli, and Licia Verde, “From Primordial Black Holes Abundance to Primordial Curvature Power Spectrum (and back),” JCAP 10, 031 (2019), arXiv:1908.03596 [astro-ph.CO] .
  * Musco (2019) Ilia Musco, “Threshold for primordial black holes: Dependence on the shape of the cosmological perturbations,” Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127 [gr-qc] .
  * Young  _et al._ (2019) Sam Young, Ilia Musco, and Christian T. Byrnes, “Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation,” JCAP 11, 012 (2019), arXiv:1904.00984 [astro-ph.CO] .
  * Kawasaki and Nakatsuka (2019) Masahiro Kawasaki and Hiromasa Nakatsuka, “Effect of nonlinearity between density and curvature perturbations on the primordial black hole formation,” Physical Review D 99 (2019), 10.1103/physrevd.99.123501.
  * Luca  _et al._ (2019) V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto, and C. Ünal, “The ineludible non-gaussianity of the primordial black hole abundance,” Journal of Cosmology and Astroparticle Physics 2019, 048–048 (2019).
  * Young and Byrnes (2020b) Sam Young and Christian T. Byrnes, “Initial clustering and the primordial black hole merger rate,” JCAP 03, 004 (2020b), arXiv:1910.06077 [astro-ph.CO] .
  * Nakamura  _et al._ (1997) Takashi Nakamura, Misao Sasaki, Takahiro Tanaka, and Kip S. Thorne, “Gravitational waves from coalescing black hole MACHO binaries,” Astrophys. J. Lett. 487, L139–L142 (1997), arXiv:astro-ph/9708060 .
  * Ioka  _et al._ (1998) Kunihito Ioka, Takeshi Chiba, Takahiro Tanaka, and Takashi Nakamura, “Black hole binary formation in the expanding universe: Three body problem approximation,” Phys. Rev. D 58, 063003 (1998), arXiv:astro-ph/9807018 .
  * Ali-Haïmoud  _et al._ (2017) Yacine Ali-Haïmoud, Ely D. Kovetz, and Marc Kamionkowski, “Merger rate of primordial black-hole binaries,” Phys. Rev. D 96, 123523 (2017), arXiv:1709.06576 [astro-ph.CO] .
  * Raidal  _et al._ (2019) Martti Raidal, Christian Spethmann, Ville Vaskonen, and Hardi Veermäe, “Formation and evolution of primordial black hole binaries in the early universe,” Journal of Cosmology and Astroparticle Physics 2019, 018–018 (2019).
  * Liu  _et al._ (2019) Lang Liu, Zong-Kuan Guo, and Rong-Gen Cai, “Effects of the surrounding primordial black holes on the merger rate of primordial black hole binaries,” Phys. Rev. D 99, 063523 (2019), arXiv:1812.05376 [astro-ph.CO] .
  * Bird  _et al._ (2016) Simeon Bird, Ilias Cholis, Julian B. Muñoz, Yacine Ali-Haïmoud, Marc Kamionkowski, Ely D. Kovetz, Alvise Raccanelli, and Adam G. Riess, “Did LIGO detect dark matter?” Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Korol  _et al._ (2020) Valeriya Korol, Ilya Mandel, M. Coleman Miller, Ross P. Church, and Melvyn B. Davies, “Merger rates in primordial black hole clusters without initial binaries,” Mon. Not. Roy. Astron. Soc. 496, 994–1000 (2020), arXiv:1911.03483 [astro-ph.HE] .
  * Franciolini  _et al._ (2022a) Gabriele Franciolini, Konstantinos Kritos, Emanuele Berti, and Joseph Silk, “Primordial black hole mergers from three-body interactions,” (2022a), arXiv:2205.15340 [astro-ph.CO] .
  * Raidal  _et al._ (2019) Martti Raidal, Christian Spethmann, Ville Vaskonen, and Hardi Veermäe, “Formation and evolution of primordial black hole binaries in the early universe,” JCAP 2019, 018 (2019), arXiv:1812.01930 [astro-ph.CO] .
  * Hütsi  _et al._ (2021) Gert Hütsi, Martti Raidal, Ville Vaskonen, and Hardi Veermäe, “Two populations of ligo-virgo black holes,” Journal of Cosmology and Astroparticle Physics 2021, 068 (2021).
  * Phukon  _et al._ (2021) Khun Sang Phukon, Gregory Baltus, Sarah Caudill, Sebastien Clesse, Antoine Depasse, Maxime Fays, Heather Fong, Shasvath J. Kapadia, Ryan Magee, and Andres Jorge Tanasijczuk, “The hunt for sub-solar primordial black holes in low mass ratio binaries is open,” (2021), arXiv:2105.11449 [astro-ph.CO] .
  * Vaskonen and Veermäe (2020) Ville Vaskonen and Hardi Veermäe, “Lower bound on the primordial black hole merger rate,” Phys. Rev. D 101, 043015 (2020), arXiv:1908.09752 [astro-ph.CO] .
  * De Luca  _et al._ (2020b) V. De Luca, V. Desjacques, G. Franciolini, and A. Riotto, “The clustering evolution of primordial black holes,” JCAP 11, 028 (2020b), arXiv:2009.04731 [astro-ph.CO] .
  * Tkachev  _et al._ (2020) Maxim Tkachev, Sergey Pilipenko, and Gustavo Yepes, “Dark Matter Simulations with Primordial Black Holes in the Early Universe,” Mon. Not. Roy. Astron. Soc. 499, 4854–4862 (2020), arXiv:2009.07813 [astro-ph.CO] .
  * Gow  _et al._ (2020) Andrew D. Gow, Christian T. Byrnes, Alex Hall, and John A. Peacock, “Primordial black hole merger rates: distributions for multiple LIGO observables,” JCAP 01, 031 (2020), arXiv:1911.12685 [astro-ph.CO] .
  * De Luca  _et al._ (2020c) V. De Luca, G. Franciolini, P. Pani, and A. Riotto, “Primordial Black Holes Confront LIGO/Virgo data: Current situation,” JCAP 06, 044 (2020c), arXiv:2005.05641 [astro-ph.CO] .
  * Hall  _et al._ (2020) Alex Hall, Andrew D. Gow, and Christian T. Byrnes, “Bayesian analysis of LIGO-Virgo mergers: Primordial vs. astrophysical black hole populations,” Phys. Rev. D 102, 123524 (2020), arXiv:2008.13704 [astro-ph.CO] .
  * García-Bellido  _et al._ (2021) Juan García-Bellido, José Francisco Nuño Siles, and Ester Ruiz Morales, “Bayesian analysis of the spin distribution of LIGO/Virgo black holes,” Phys. Dark Univ. 31, 100791 (2021), arXiv:2010.13811 [astro-ph.CO] .
  * Wong  _et al._ (2021) Kaze W. K. Wong, Gabriele Franciolini, Valerio De Luca, Vishal Baibhav, Emanuele Berti, Paolo Pani, and Antonio Riotto, “Constraining the primordial black hole scenario with Bayesian inference and machine learning: the GWTC-2 gravitational wave catalog,” Phys. Rev. D 103, 023026 (2021), arXiv:2011.01865 [gr-qc] .
  * De Luca  _et al._ (2021) V. De Luca, G. Franciolini, P. Pani, and A. Riotto, “Bayesian Evidence for Both Astrophysical and Primordial Black Holes: Mapping the GWTC-2 Catalog to Third-Generation Detectors,” JCAP 05, 003 (2021), arXiv:2102.03809 [astro-ph.CO] .
  * Franciolini  _et al._ (2022b) Gabriele Franciolini, Vishal Baibhav, Valerio De Luca, Ken K. Y. Ng, Kaze W. K. Wong, Emanuele Berti, Paolo Pani, Antonio Riotto, and Salvatore Vitale, “Searching for a subpopulation of primordial black holes in LIGO-Virgo gravitational-wave data,” Phys. Rev. D 105, 083526 (2022b), arXiv:2105.03349 [gr-qc] .
  * Chen  _et al._ (2022) Zu-Cheng Chen, Chen Yuan, and Qing-Guo Huang, “Confronting the primordial black hole scenario with the gravitational-wave events detected by LIGO-Virgo,” Phys. Lett. B 829, 137040 (2022), arXiv:2108.11740 [astro-ph.CO] .
  * Kavanagh (2022) Bradley J. Kavanagh, “PBHbounds [Code v1.0, accessed 03/03/2022],” https://github.com/bradkav/PBHbounds, DOI:10.5281/zenodo.3538998 (2022).
  * Phinney (2001) E. S. Phinney, “A Practical theorem on gravitational wave backgrounds,” (2001), arXiv:astro-ph/0108028 .
  * Zhu  _et al._ (2011) Xing-Jiang Zhu, E. Howell, T. Regimbau, D. Blair, and Zong-Hong Zhu, “Stochastic Gravitational Wave Background from Coalescing Binary Black Holes,” Astrophys. J. 739, 86 (2011), arXiv:1104.3565 [gr-qc] .
  * Allsman  _et al._ (2001) R. A. Allsman _et al._ (Macho), “MACHO project limits on black hole dark matter in the 1-30 solar mass range,” Astrophys. J. Lett. 550, L169 (2001), arXiv:astro-ph/0011506 .
  * Tisserand  _et al._ (2007) P. Tisserand _et al._ (EROS-2), “Limits on the Macho Content of the Galactic Halo from the EROS-2 Survey of the Magellanic Clouds,” Astron. Astrophys. 469, 387–404 (2007), arXiv:astro-ph/0607207 .
  * Niikura  _et al._ (2019a) Hiroko Niikura _et al._ , “Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations,” Nature Astron. 3, 524–534 (2019a), arXiv:1701.02151 [astro-ph.CO] .
  * Niikura  _et al._ (2019b) Hiroko Niikura, Masahiro Takada, Shuichiro Yokoyama, Takahiro Sumi, and Shogo Masaki, “Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events,” Phys. Rev. D 99, 083503 (2019b), arXiv:1901.07120 [astro-ph.CO] .
  * Carr  _et al._ (2017) Bernard Carr, Martti Raidal, Tommi Tenkanen, Ville Vaskonen, and Hardi Veermäe, “Primordial black hole constraints for extended mass functions,” Phys. Rev. D 96, 023514 (2017), arXiv:1705.05567 [astro-ph.CO] .
  * Green and Kavanagh (2021) Anne M. Green and Bradley J. Kavanagh, “Primordial black holes as a dark matter candidate,” Journal of Physics G Nuclear Physics 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Nitz and Wang (2021a) Alexander H. Nitz and Yi-Fan Wang, “Search for Gravitational Waves from High-Mass-Ratio Compact-Binary Mergers of Stellar Mass and Subsolar Mass Black Holes,” Phys. Rev. Lett. 126, 021103 (2021a), arXiv:2007.03583 [astro-ph.HE] .
  * Nitz and Wang (2021b) Alexander H. Nitz and Yi-Fan Wang, “Search for Gravitational Waves from the Coalescence of Subsolar-Mass Binaries in the First Half of Advanced LIGO and Virgo’s Third Observing Run,” Phys. Rev. Lett. 127, 151101 (2021b), arXiv:2106.08979 [astro-ph.HE] .
  * Barsanti  _et al._ (2021) Susanna Barsanti, Valerio De Luca, Andrea Maselli, and Paolo Pani, “Detecting Subsolar-Mass Primordial Black Holes in Extreme Mass-Ratio Inspirals with LISA and Einstein Telescope,” (2021), arXiv:2109.02170 [gr-qc] .
  * Kavanagh  _et al._ (2018b) Bradley J. Kavanagh, Daniele Gaggero, and Gianfranco Bertone, “Merger rate of a subdominant population of primordial black holes,” Physical Review D 98 (2018b), 10.1103/physrevd.98.023536.
  * van den Bosch  _et al._ (2017) Frank C van den Bosch, Go Ogiya, Oliver Hahn, and Andreas Burkert, “Disruption of dark matter substructure: fact or fiction?” Monthly Notices of the Royal Astronomical Society 474, 3043–3066 (2017).
  * Adamek  _et al._ (2019) Julian Adamek, Christian T. Byrnes, Mateja Gosenca, and Shaun Hotchkiss, “WIMPs and stellar-mass primordial black holes are incompatible,” Phys. Rev. D 100, 023506 (2019), arXiv:1901.08528 [astro-ph.CO] .
  * Sadeghian  _et al._ (2013) Laleh Sadeghian, Francesc Ferrer, and Clifford M. Will, “Dark matter distributions around massive black holes: A general relativistic analysis,” Phys. Rev. D 88, 063522 (2013), arXiv:1305.2619 [astro-ph.GA] .
  * Ferrer  _et al._ (2017) Francesc Ferrer, Augusto Medeiros da Rosa, and Clifford M. Will, “Dark matter spikes in the vicinity of Kerr black holes,” Phys. Rev. D 96, 083014 (2017), arXiv:1707.06302 [astro-ph.CO] .
  * Lewis  _et al._ (2000) Antony Lewis, Anthony Challinor, and Anthony Lasenby, “Efficient computation of CMB anisotropies in closed FRW models,” Astrophys. J. 538, 473–476 (2000), arXiv:astro-ph/9911177 .
  * Howlett  _et al._ (2012) Cullan Howlett, Antony Lewis, Alex Hall, and Anthony Challinor, “CMB power spectrum parameter degeneracies in the era of precision cosmology,” JCAP 2012, 027 (2012), arXiv:1201.3654 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Kavanagh  _et al._ (2018) Bradley J. Kavanagh, Daniele Gaggero, and Gianfranco Bertone, “Merger rate of a subdominant population of primordial black holes,” Phys. Rev. D 98, 023536 (2018), arXiv:1805.09034 [astro-ph.CO] .
  * Chandrasekhar (1943a) S. Chandrasekhar, “Dynamical friction. i. general considerations: the coefficient of dynamical friction.” The Astrophysical Journal 97, 255 (1943a).
  * Chandrasekhar (1943b) S. Chandrasekhar, “Dynamical friction. II. the rate of escape of stars from clusters and the evidence for the operation of dynamical friction.” The Astrophysical Journal 97, 263 (1943b).
  * Chandrasekhar (1943c) S. Chandrasekhar, “Dynamical friction. III. a more exact theory of the rate of escape of stars from clusters.” The Astrophysical Journal 98, 54 (1943c).
  * Binney and Tremaine (2008) J. Binney and S. Tremaine,  _Galactic Dynamics: Second Edition_ (Princeton University Press, 2008).
  * Kavanagh (2020) Bradley J. Kavanagh, “HaloFeedback [Code, v0.9],” https://github.com/bradkav/HaloFeedback, DOI:10.5281/zenodo.3688813 (2020).
  * Hild  _et al._ (2011) S. Hild _et al._ , “Sensitivity Studies for Third-Generation Gravitational Wave Observatories,” Class. Quant. Grav. 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Evans  _et al._ (2016) Matthew Evans, Jan Harms, and Salvatore Vitale, “Exploring the sensitivity of next generation gravitational wave detectors,” https://dcc.ligo.org/LIGO-P1600143/public (2016).
  * Skilling (2004) John Skilling, “Nested sampling,” in  _AIP Conference Proceedings_ (AIP, 2004).
  * Higson  _et al._ (2019) Edward Higson, Will Handley, Mike Hobson, and Anthony Lasenby, “Dynamic nested sampling: an improved algorithm for parameter estimation and evidence calculation,” Statistics and Computing 29, 891–913 (2019), arXiv:1704.03459 [stat.CO] .
  * Speagle (2020) Joshua S. Speagle, “DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences,” Mon. Not. Roy. Astron. Soc. 493, 3132–3158 (2020), arXiv:1904.02180 [astro-ph.IM] .
  * Jeffreys (1998) H. Jeffreys,  _The Theory of Probability_, Oxford Classic Texts in the Physical Sciences (OUP Oxford, 1998).
  * Kass and Raftery (1995) Robert E. Kass and Adrian E. Raftery, “Bayes factors,” Journal of the American Statistical Association 90, 773–795 (1995).
  * Cutler  _et al._ (2019) Curt Cutler _et al._ , “What we can learn from multi-band observations of black hole binaries,” (2019), arXiv:1903.04069 [astro-ph.HE] .
  * Aasi (2014) J. et al. Aasi, “First all-sky search for continuous gravitational waves from unknown sources in binary systems,” Physical Review D 90 (2014), 10.1103/physrevd.90.062010.
  * Guo and Miller (2022) Huai-Ke Guo and Andrew Miller, “Searching for Mini Extreme Mass Ratio Inspirals with Gravitational-Wave Detectors,” (2022), arXiv:2205.10359 [astro-ph.IM] .
  * Graham  _et al._ (2020) M. J. Graham, K. E. S. Ford, B. McKernan, N. P. Ross, D. Stern, K. Burdge, M. Coughlin, S. G. Djorgovski, A. J. Drake, D. Duev, M. Kasliwal, A. A. Mahabal, S. van Velzen, J. Belecki, E. C. Bellm, R. Burruss, S. B. Cenko, V. Cunningham, G. Helou, S. R. Kulkarni, F. J. Masci, T. Prince, D. Reiley, H. Rodriguez, B. Rusholme, R. M. Smith, and M. T. Soumagnac, “Candidate electromagnetic counterpart to the binary black hole merger gravitational-wave event s190521g,” Phys. Rev. Lett. 124, 251102 (2020).
  * Baumann  _et al._ (2019) Daniel Baumann, Horng Sheng Chia, and Rafael A. Porto, “Probing Ultralight Bosons with Binary Black Holes,” Phys. Rev. D 99, 044001 (2019), arXiv:1804.03208 [gr-qc] .
  * Baumann  _et al._ (2022a) Daniel Baumann, Gianfranco Bertone, John Stout, and Giovanni Maria Tomaselli, “Ionization of gravitational atoms,” Phys. Rev. D 105, 115036 (2022a), arXiv:2112.14777 [gr-qc] .
  * Baumann  _et al._ (2022b) Daniel Baumann, Gianfranco Bertone, John Stout, and Giovanni Maria Tomaselli, “Sharp Signals of Boson Clouds in Black Hole Binary Inspirals,” Phys. Rev. Lett. 128, 221102 (2022b), arXiv:2206.01212 [gr-qc] .
  * Raidal  _et al._ (2017) Martti Raidal, Ville Vaskonen, and Hardi Veermäe, “Gravitational waves from primordial black hole mergers,” Journal of Cosmology and Astroparticle Physics 2017, 037–037 (2017).
  * Chernoff and Finn (1993) David F. Chernoff and Lee Samuel Finn, “Gravitational radiation, inspiraling binaries, and cosmology,” Astrophys. J. Lett. 411, L5–L8 (1993), arXiv:gr-qc/9304020 .
  * Ajith  _et al._ (2008) P. Ajith _et al._ , “A Template bank for gravitational waveforms from coalescing binary black holes. I. Non-spinning binaries,” Phys. Rev. D 77, 104017 (2008), [Erratum: Phys.Rev.D 79, 129901 (2009)], arXiv:0710.2335 [gr-qc] .
  * Coogan (2021) Adam Coogan, “pydd [Code],” https://github.com/adam-coogan/pydd (2021).