# Measuring properties of primordial black hole mergers at cosmological
distances: effect of higher order modes in gravitational waves

Ken K. Y. Ng<EMAIL_ADDRESS>William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, USA  LIGO,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, Massachusetts 02139, USA  Boris Goncharov  Gran
Sasso Science Institute (GSSI), I-67100 L’Aquila, Italy  INFN, Laboratori
Nazionali del Gran Sasso, I-67100 Assergi, Italy  Shiqi Chen  LIGO,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, Massachusetts 02139, USA  Ssohrab Borhanian
Institute for Gravitation and the Cosmos, Department of Physics, Pennsylvania
State University, University Park, PA, 16802, USA  Theoretisch-Physikalisches
Institut, Friedrich-Schiller-Universität Jena, 07743, Jena, Germany  Ulyana
Dupletsa  Gran Sasso Science Institute (GSSI), I-67100 L’Aquila, Italy  INFN,
Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy  Gabriele
Franciolini  Dipartimento di Fisica, Sapienza Università di Roma, Piazzale
Aldo Moro 5, 00185, Roma, Italy  INFN, Sezione di Roma, Piazzale Aldo Moro 2,
00185, Roma, Italy  Marica Branchesi  Gran Sasso Science Institute (GSSI),
I-67100 L’Aquila, Italy  INFN, Laboratori Nazionali del Gran Sasso, I-67100
Assergi, Italy  Jan Harms  Gran Sasso Science Institute (GSSI), I-67100
L’Aquila, Italy  INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi,
Italy  Michele Maggiore  Département de Physique Théorique, Université de
Genève, 24 quai E. Ansermet, CH-1211 Geneva, Switzerland  Gravitational Wave
Science Center (GWSC), Université de Genève, CH-1211 Geneva, Switzerland
Antonio Riotto  Département de Physique Théorique, Université de Genève, 24
quai E. Ansermet, CH-1211 Geneva, Switzerland  Gravitational Wave Science
Center (GWSC), Université de Genève, CH-1211 Geneva, Switzerland  B. S.
Sathyaprakash  Institute for Gravitation and the Cosmos, Department of
Physics, Pennsylvania State University, University Park, PA, 16802, USA
Department of Astronomy & Astrophysics, Pennsylvania State University,
University Park, PA, 16802, USA  School of Physics and Astronomy, Cardiff
University, Cardiff, UK, CF24 3AA  Salvatore Vitale  LIGO, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA  Kavli Institute
for Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA

###### Abstract

Primordial black holes (PBHs) may form from the collapse of matter
overdensities shortly after the Big Bang. One may identify their existence by
observing gravitational wave (GW) emissions from merging PBH binaries at high
redshifts $z\gtrsim 30$, where astrophysical binary black holes (BBHs) are
unlikely to merge. The next-generation ground-based GW detectors, Cosmic
Explorer and Einstein Telescope, will be able to observe BBHs with total
masses of $\mathcal{O}(10-100)~{}M_{\odot}$ at such redshifts. This paper
serves as a companion paper of Ref. [1], focusing on the effect of higher-
order modes (HoMs) in the waveform modeling, which may be detectable for these
high redshift BBHs, on the estimation of source parameters. We perform
Bayesian parameter estimation to obtain the measurement uncertainties with and
without HoM modeling in the waveform for sources with different total masses,
mass ratios, orbital inclinations and redshifts observed by a network of next-
generation GW detectors. We show that including HoMs in the waveform model
reduces the uncertainties of redshifts and masses by up to a factor of two,
depending on the exact source parameters. We then discuss the implications for
identifying PBHs with the improved single-event measurements, and expand the
investigation of the model dependence of the relative abundance between the
BBH mergers originating from the first stars and the primordial BBH mergers as
shown in Ref. [1].

††preprint: ET-0231A-22, CE-P2200006

##  I Introduction

An interesting possibility is that a fraction of the merger events detected by
the LIGO-Virgo-KAGRA (LVK) Collaboration may be due to primordial BHs (PBHs)
[2, 3, 4, 5] formed from the collapse of sizable overdensities in the
radiation-dominated early universe [6, 7, 8, 9]. In this scenario, PBHs are
not clustered at formation [10, 11, 12, 13, 14, 15], they are born spinless
[16, 17] and may assemble in binaries via gravitational decoupling from the
Hubble flow before the matter-radiation equality [18, 19] (see [20, 21, 22,
23, 24] for reviews). After their formation, PBH binaries may be affected by a
phase of baryonic mass accretion at redshifts smaller than $z\sim 30$, which
would modify the PBH masses, spins and merger rate [25, 26].

Analysing the population properties of masses, spins, and redshifts of binary
black holes (BBHs) in the LVK’s second catalog [27], several studies
constrained the potential contribution from PBH binaries to current data [28,
29, 26, 30, 31, 32, 33, 34, 35]. However, these analyses require precise
knowledge of the astrophysical BBH “foreground” in order to verify if there is
a PBH subpopulation within the BBHs observed at low redshifts [33, 36]. Such
analyses are limited by the horizon of current GW detectors, $z\lesssim 3$ at
their design sensitivity [37], and are subject to significant uncertainties on
the mechanisms of BBH formation in different astrophysical environments, such
as galactic fields [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48], dense star
clusters [49, 50, 51, 52, 53, 54, 55, 56, 52, 53, 54, 50], active galactic
nuclei [57, 58, 59, 60, 61, 62, 63, 64, 65], or from the collapse of
Population III (Pop III) stars [66, 67, 68, 69].

Instead, searching for PBHs at high redshifts where astrophysical BHs have not
merged yet may mitigate most of the issues caused by the astrophysical
foreground. The PBH merger rate increases with redshift [70], while the
astrophysical contribution decreases at $z\lesssim 30$ [69, 71, 72]. The
proposed next-generation detectors, such as the Cosmic Explorer (CE) [73, 74,
75] and the Einstein Telescope (ET) [76, 77], whose horizons are up to $z\sim
100$ for stellar-mass BBHs [37, 78], may provide a unique opportunity to test
and shed light on the primordial origin of BH mergers at high redshifts. A key
question is therefore to understand the uncertainties related to the
measurements of the source parameters, such as the redshift, masses and spins.

In Ref. [1], we established the possibility of identifying the PBH mergers
with masses of 20 and 40 $M_{\odot}$ at $z\geq 40$ using single-event redshift
measurements. We also discussed how the prior knowledge of relative abundance
between Pop III and PBH mergers affects the statistical significance, assuming
that there is a critical redshift, $z_{\mathrm{crit}}=30$, above which no
astrophysical BBHs are expected to merge. The results were based on full
Bayesian parameter estimation with a waveform model, IMRPhenomXPHM, which
includes the effects of spin precession and higher-order modes (HoMs) [79, 80,
81]. In this paper, we show the importance of HoMs to the parameter estimation
of the high redshift BBHs at $z\geq 10$ in the context of PBH detections. We
compare the Bayesian posteriors of the relevant parameters obtained by
IMRPhenomXPHM and the similar waveform family without HoMs, IMRPhenomPv2 [82,
83, 84] to systematically study the improvement on measurements due to the HoM
modeling in the waveform.

We first recap the details of our simulations and the settings of the
parameter estimation in Sec. II. Then, we show whether and how IMRPhenomXPHM
performs better when measuring redshift (Sec. III), as well as masses and
spins (Sec. IV), for BBHs with different sets of the source-frame total mass,
mass ratio, orbital inclination, and redshift. Finally, in Sec. V, we re-
examine the estimation of the probability that a single source originated from
PBHs using redshift measurements under different choices of
$z_{\mathrm{crit}}$, and discuss the possible implications of the mass and
spin measurements for PBH detections.

##  II Simulation details

As in Ref. [1], we simulate BBHs at five different redshifts, $\hat{z}=10,$
20, 30, 40 and 50. The hat symbol denotes the true value of a parameter here
and throughout the paper. To encompass the detectable mass range, we choose
the total masses in the source frame to be $\hat{M}_{\rm tot}=5$, 10, 20, 40,
and 80 $M_{\odot}$, with mass ratios $\hat{q}=1$, 2, 3, 4 and 5. Here, we
define $q\equiv m_{1}/m_{2}$ for $m_{1}>m_{2}$, where $m_{1}$ and $m_{2}$ are
the primary and secondary mass, respectively. For each mass pair, we further
choose four orbital inclination angles, $\hat{\iota}=0$ (face-on), $\pi/6$,
$\pi/3$, and $\pi/2$ (edge-on). All simulated BBHs are non-spinning, as we
expect that PBHs are born with negligible spins [85, 17, 16] and may be spun-
up by accreting materials at later times [85, 25, 86, 26]. However, we do
_not_ assume zero spins when performing parameter estimation of the source
parameters and instead allow for generic spin-precession. For each of these
500 sources, the sky location and polarization angle are chosen to maximize
the signal-to-noise ratio (SNR) for each source. The reference orbital phase
and GPS time are fixed at $0$ and $1577491218$, respectively. The baseline
detector network is a 40-km CE in the United States, 20-km CE in Australia,
and ET in Europe. We only analyze simulated sources whose network SNRs are
larger than 12. We use Planck 2018 Cosmology when calculating the luminosity
distance $d_{L}$ at a given redshift [87].

We employ a nested sampling algorithm [88, 89] packaged in Bilby [90] to
obtain posterior probability densities. As we are only interested in the
uncertainty caused by the loudness of the signal, we use a zero-noise
realization [91] for the Bayesian inference and mitigate the offsets
potentially caused by Gaussian fluctuations [92]. To ensure our results are
free from the systematics due to the difference in the two waveform families,
we use the same waveform family for both simulating the waveforms and
calculating the likelihood. That is, we use the IMRPhenomXPHM (IMRPhenomPv2)
waveform template to analyze the IMRPhenomXPHM (IMRPhenomPv2) simulated
waveforms 111See Ref. [93] for the analysis of waveform systematics for the
next-generation GW detectors.. The low-frequency cut off in the likelihood
calculations is 5 Hz for all sources.

As in Ref. [1], we first sample the parameter space with uniform priors on the
detector-frame total mass, $M_{\rm tot}^{D}=M_{\rm tot}(1+z)$, between
$[0.5,1.5]\hat{M}_{\rm tot}^{D}$, and $q$ between $[1,10]$. The prior on
redshift is uniform in the comoving rate density,
$\propto\frac{dV_{c}}{dz}\frac{1}{1+z}$, between
$[z(\hat{d}_{L}/10),z(5\hat{d}_{L})]$. In Sec. V, we will revisit the
physically motivated prior on redshift. We use uniform priors for other
parameters: the sky position, the polarization angle, the orbital inclination,
the spin orientations, the spin magnitudes, the arrival time and the phase of
the signal at the time of arrival.

Then, we reweigh the posteriors into uniform prior on the source-frame primary
mass, $m_{1}$, and the inverse mass ratio $1/q$ (which is between $[0.1,1]$).
Strictly speaking, the marginalized one-dimensional priors on $m_{1}$ and
$1/q$ are not exactly uniform after the reweighing because the boundary of the
square domain of $(M_{\rm tot}^{D},q)$ transforms into a different shape
according to the Jacobian. For example, the marginalized prior on the redshift
and that on the inverse mass ratio have additional factors of $1/(1+z)$ and
$q/(q+1)$, respectively, upon the coordinate transformation. However, we find
that such boundary effect has negligible effect on the posteriors. As we will
discuss below, the degeneracy among different parameters and the scaling in
$p_{0}(z)$ is more significant.

##  III Redshift measurement in presence of higher-order modes

Since each HoM has a different angular emission spectrum, including HoMs in
the waveform model breaks the distance-inclination degeneracy characteristic
of the dominant $(2,2)$ harmonic mode [94, 95]. The interference of additional
HoMs can result in amplitude modulation, similar to what can be induced by
spin precession [96, 81, 80]. For example, in the top panel of Fig. 1 we show
the Fourier amplitude of a BBH with $(\hat{M}_{\rm
tot},\hat{z},\hat{q})=(80~{}M_{\odot},30,1)$ and
$\hat{\iota}=0^{\circ},30^{\circ},60^{\circ}$ and $90^{\circ}$ (blue, orange,
green and red, respectively). To reduce the systematics between waveform
families due to differences in precessing frame mapping, we compare
IMRPhenomXPHM (solid lines) and IMRPhenomXP [79] (dotted lines) instead. The
amplitude modulation of the waveforms with HoMs – which is stronger for
inclination angles close to $90^{\circ}$ – is apparent and helps improving the
estimation of the distance and the inclination. By contrast, for the waveforms
without HoMs the main effect of increasing the inclination angle is to reduce
the Fourier amplitude, which qualitatively shows why the two parameters are
partially degenerate when only the (2,2) mode is used222The degeneracy is
worst at small inclination angles, see Sec. III.1.. The other contribution is
the phase modulation in the later part of the waveform due to HoMs. To
visualize this effect, we show the phase difference between IMRPhenomXP and
IMRPhenomXPHM, $\Delta\Phi(f)$, in the bottom panel of Fig. 1. Whereas the
(2,2) mode of the inspiral defines the waveform up to $\approx 8$ Hz, after
which the ringdown takes over, HoMs of the inspiral extend to higher
frequencies. The interference of the HoMs and the (2,2) ringdown piles up a
significant phase modulation, and hence improves the measurement of
inclination.

(a)

(b)

Figure 1: Comparison between waveforms without HoMs (IMRPhenomXP) and with
HoMs (IMRPhenomXP) for BBHs with $(\hat{M}_{\rm
tot},\hat{z},\hat{q})=(80~{}M_{\odot},30,1)$ and
$\hat{\iota}=0^{\circ},30^{\circ},60^{\circ}$ and $90^{\circ}$ (blue, orange,
green and red, respectively). Top panel: Strain amplitudes for IMRPhenomXPHM
(solid lines) and IMRPhenomXP (dotted lines) projected on ET’s detector frame.
Bottom panel: Phase difference between IMRPhenomXPHM and IMRPhenomXP at each
frequency. In all systems, the right ascension angle, declination angle and
polarization angle are $110^{\circ}$, $45^{\circ}$ and $93^{\circ}$,
respectively.

Moreover, the parameters, $q$, $\iota$ and $M_{\rm tot}$, determine the
amplitude of each mode. The uncertainties of distance (and hence redshift) are
thus sensitive to the values of $(q,\iota,M_{\rm tot})$ with other parameters
fixed. In this section, we will quantify the variation of the redshift
uncertainty due to each intrinsic parameter one at a time. We will also show
which region of redshift gain the most from the presence of HoMs in the
waveform model. In the following figures, blue (red) violins represent the
posteriors obtained by IMRPhenomXPHM (IMRPhenomPv2).

###  III.1 Orbital inclination

We first discuss the role of orbital inclination in the redshift measurements
using HoM waveforms. Waveform models which only contain (2, 2) mode suffer
from the distance-inclination degeneracy of the mode, especially for nearly
face-on $(\iota\simeq 0)$ systems whose amplitude scales as
$\sim(1-\iota^{2}/2)/d_{L}$. On the other hand, each HoM corresponds to
spherical harmonics with a different angular response as a function of
$\iota$. If the waveform models are sensitive to HoMs, measuring the relative
amplitudes of the HoMs provides better constraints on the orbital inclination
angle, and thus reduces the distance-inclination degeneracy.

Figure 2: Posteriors of redshift for sources with $(\hat{M}_{\rm
tot},\hat{z},\hat{q})=(40M_{\odot},30,1)$ at
$\hat{\iota}=0^{\circ},30^{\circ},60^{\circ}$ and $90^{\circ}$, obtained with
HoM (blue, IMRPhenomXPHM) and without HoM (red, IMRPhenomPv2). The solid
horizontal lines show the 95% credible intervals, whereas the dashed lines
mark $\hat{z}$. Top axis shows the optimal SNR of IMRPhenomXPHM waveform,
$\rho_{\rm PHM}$, and that of IMRPhenomPv2 waveform, $\rho_{\rm P}$. The
detector network is CE-CES20-ET. We also note that the SNR does not
necessarily increase with the addition of HoM because they may interfere
destructively at large values of $\hat{\iota}$.

We now quantify the improvement on the redshift measurements due to the
presence of HoMs with varying inclination angles. In Fig. 2, we show the
redshift posteriors obtained by the two waveform models for sources with
$(\hat{M}_{\rm tot},\hat{z},\hat{q})=(40M_{\odot},30,1)$ at
$\hat{\iota}=0^{\circ},30^{\circ},60^{\circ}$ and $90^{\circ}$. Indeed, the
redshift uncertainties are generally smaller in the cases of IMRPhenomXPHM
than those of IMRPhenomPv2. The decrease in the uncertainties is about $\sim
30\%-50\%$.

Notably, the lower bound of redshift uncertainties increases from $z\sim 10$
in the cases of IMRPhenomPv2 to $z\sim 20$ in the cases of IMRPhenomXPHM. This
improvement due to HoMs is particularly interesting for determining the
astrophysical or primordial origin of BBHs. If the redshift measurement of a
system is precise enough to rule out the epoch of the astrophysical BBHs, one
may even use a single measurement to identify the existence of primordial BBH,
as discussed in Ref. [1]. On the other hand, even if a single measurement is
not conclusive enough, the improvement on redshift measurements reduces the
required number of events to conduct a statistical analysis over a population
of high-redshift BBHs [97, 98].

As $\hat{\iota}$ increases, the SNR decreases with the amplitude of the
dominant $(2,2)$ harmonic. One might expect the redshift uncertainty would
then increase with $\hat{\iota}$. Indeed, as shown in Fig. 2, the uncertainty
raises from $\hat{\iota}=0^{\circ}$ to $\hat{\iota}=60^{\circ}$ but shrinks
when the system is edge on $(\hat{\iota}=90^{\circ})$. For the edge-on system,
the $\times$-polarization content is very sensitive to small changes of
$\iota$, which makes it possible to obtain better estimates of $\iota$ and to
break the distance-inclination degeneracy. The other feature worth-mentioning
is the apparent bias in face-on systems. As discussed in Ref. [1], this is
caused by the physical cut-off of the parameter space: an overestimation of
the redshift cannot be compensated by the increase of $\cos\iota$ beyond 1.

We also note that the HoM-improvements and features of the redshift
uncertainties discussed above are similar for systems with other values of
$\hat{q}$ and $\hat{M}_{\rm tot}$.

###  III.2 Mass ratio

It is well known that the excitation of each HoM is sensitive to the mass
ratio of a binary system. In particular, Ref. [99] identified the mass ratio
dependence of each mode amplitude in Post-Newtonian (PN) theory, which was
later found to be in a good agreement with the full numerical relativity
simulations [100]. Broadly speaking, more HoMs are excited as the system is
more asymmetric in component masses. When the waveform model is sensitive to
identify the presence (systems with more asymmetric masses) or the absence
(systems with nearly equal masses), it can provide extra constraints on the
parameter space along the degeneracy between inclination and distance. Thus,
the redshift uncertainties are generally smaller when one uses the HoM-
waveform template in the inference, even if the true waveform does not contain
substantial HoM contents.

Figure 3: Redshift measurements for sources with $(\hat{M}_{\rm
tot},\hat{z},\hat{\iota})=(40M_{\odot},30,30^{\circ})$ at $\hat{q}=1,2,3$ and
$4$. The format is the same as in Fig. 2.

In Fig. 3, we show the redshift posteriors obtained by the two waveform models
for sources with $(\hat{M}_{\rm
tot},\hat{z},\hat{\iota})=(40M_{\odot},30,30^{\circ})$ at $\hat{q}=1,2,3$ and
4. First, for all mass ratios, the redshift uncertainties obtained by
IMRPhenomXPHM decrease by $\sim 20\%$ to $\sim 40\%$ when compared to those
obtained by IMRPhenomPv2. Second, the scaling between the uncertainties and
the SNRs in the IMRPhenomXPHM is different from that in IMRPhenomPv2. While
the increase in $q$ leads to an increase in the redshift uncertainty due to
the decrease in the SNR, it also excites more the HoMs and enrich the
structure in the IMRPhenomXPHM waveform. The additional HoM contents may
provide more constraining power to break the distance-inclination degeneracy.
Hence, the redshift uncertainties obtained by IMRPhenomXPHM do not increase
with SNR as fast as those obtained by IMRPhenomPv2.

###  III.3 Total mass

Unlike mass ratio and inclination, the total mass does not change the HoM
contents, i.e., it simply scales the overall amplitudes and frequencies of all
modes. However, from the perspective of detection, it is a crucial parameter
that determines the detectability of HoMs, which then affects the ability of
using HoMs to break the distance-inclination degeneracy. HoMs have higher
frequency, but weaker amplitudes, than the dominating $(2,2)$ mode. If the
power of all HoMs is much weaker than the given noise floor, only the $(2,2)$
mode is detectable and the distance-inclination degeneracy still remains,
regardless of the presence of HoMs in the waveform model.

Figure 4: Redshift measurements for sources with
$(\hat{q},\hat{z},\hat{\iota})=(1,30,30^{\circ})$ at $\hat{M}_{\rm
tot}=10,20,40$ and $80~{}M_{\odot}$. The format is the same as in Fig. 2.

As an example in Fig. 4, we show the posterior for the redshift obtained with
the two waveform models for sources with
$(\hat{q},\hat{z},\hat{\iota})=(1,30,30^{\circ})$ at $\hat{M}_{\rm
tot}=10,20,30$ and $40M_{\odot}$. Indeed, the redshift uncertainty for
$\hat{M}_{\rm tot}=10M_{\odot}$ obtained by IMRPhenomXPHM is comparable to
those obtained by IMRPhenomPv2. On the other hand, as $\hat{M}_{\rm tot}$
increases, the contribution of HoMs becomes significant and breaks the
distance-inclination degeneracy. The redshift uncertainties in the presence of
HoMs can be reduced by $\sim 30-40\%$ when $\hat{M}_{\rm tot}$ increases from
$40$ to $80M_{\odot}$. As $\hat{M}_{\rm tot}$ increases further, the waveform
drifts out of the sensitive frequency band, and the signal is not detectable.

###  III.4 Redshift

Similar to the total mass, the redshift only scales the overall amplitudes of
all modes. As the detector-frame total mass increases with the redshift, we
expect the improvement on the redshift measurement to be larger at higher
redshift. On the other hand, the SNR decreases as the redshift increases, and
the overall uncertainties increase in both waveform models. These two trends
can be seen in Fig. 5, in which we compare the posteriors of redshifts
obtained by the two waveform models for sources with $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota})=(40M_{\odot},1,30^{\circ})$ at
$\hat{z}=10,20,30,40$. As the redshift increases, the posteriors of redshift
in the cases of IMRPhenomXPHM are centered at the true values, while those of
IMRPhenomPv2 shift towards lower redshift. This is a combined effect of the
redshift prior and the distance-inclination degeneracy. In the matter-
dominated regime, $1\lesssim z\lesssim 1000$, this prior scales as
$p_{0}(z)\sim(1+z)^{-5/2}$ and favors smaller redshift. Due to the distance-
inclination degeneracy in IMRPhenomPv2, the posterior leans towards the region
of larger inclination angle but smaller distance. On the other hand, owing to
the presence of HoMs, the posteriors obtained with IMRPhenomXPHM do not suffer
from this degeneracy and hence are less influenced by the inverse power-law
scaling in $p_{0}(z)$.

Figure 5: Redshift measurements for sources with $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota})=(40M_{\odot},1,30^{\circ})$ at $\hat{z}=10,20,30,40$
and $50$. We offset the posterior to the true redshift $\hat{z}$. The format
is the same as in Fig. 2.

###  III.5 Correlation between distance, inclination and mass ratio

We end this section by discussing how the presence of HoM aids the measurement
in $(q,\iota,z)$. We illustrate the two-dimensional (2D) posteriors among the
pairs of $(q,\cos\iota,z)$, as shown in Fig. 6. In this example for the source
with $(\hat{M}_{\rm
tot},\hat{z},\hat{q},\hat{\iota})=(40M_{\odot},40,2,60^{\circ})$, the 2D
posteriors obtained by IMRPhenomXPHM (blue) show drastically different
behaviors when compared to those obtained by IMRPhenomPv2 (red).

First, the 2D contours, and hence the marginalized posteriors, are more
localized in IMRPhenomXPHM. In particular, the marginalized posteriors of
$1/q$ and $\cos\iota$ do not rail towards the prior edge at $q=\cos\iota=1$.
Second, the contours in IMRPhenomPv2 posteriors have different structures than
those in IMRPhenomXPHM posteriors. In the pair $(1/q,\cos\iota)$, the
IMRPhenomXPHM posterior shows an anti-correlation, but the IMRPhenomPv2
posterior is uncorrelated. In the pair $(\cos\iota,z)$, the IMRPhenomXPHM
posterior partially follows the correlation as in IMRPhenomPv2 posterior for
$\cos\iota\lesssim 0.6$, and becomes anti-correlated $\cos\iota\gtrsim 0.8$.
Similar situation happens in the pair $(1/q,z)$, in which the IMRPhenomXPHM
posterior becomes uncorrelated as $1/q$ increases to $\sim 0.4$. These
features emerge from the different dependence of $q$ and $\iota$ in each HoM.

Figure 6: Comparison of the posteriors among $(1/q,\cos{\iota},z)$ obtained by
IMRPhenomPv2 (red) and IMRPhenomXPHM (blue) at $(\hat{M}_{\rm
tot},\hat{z},\hat{q},\hat{\iota})=(40M_{\odot},40,2,60^{\circ})$. The contours
represent the boundaries of 68%, 95%, and 99.7% credible regions. The true
values are indicated by the black solid lines.

##  IV Mass and spin measurement

In this section, we explore the impact of HoM on the improvement of
measurements of masses and spins. These intrinsic parameters can provide
additional evidence for the (non)existence of PBHs, if they can be well
measured. For example, the spins of PBHs formed in the standard scenario from
the collapse of density perturbations generated during inflation [22] are
expected to be below the percent level, see Refs. [16, 17]. The spins of PBHs
are expected to remain negligible at $z\gtrsim 30$ because there is not enough
time for the BHs to gain angular momentum through accretion of surrounding
materials. This prediction remain valid for $z\gtrsim 10$ in the case of weak
accretion [25], which we will take as a benchmark scenario in the following.

The mass distribution of a population of PBHs depends on the properties of the
collapsing perturbations. A large class of models predicts a distribution that
can be approximated with a log-normal shape, whose central mass scale and
width are model dependent and not currently observationally constrained. If
the masses of high-redshift BBHs are measured to be contradictory to the
astrophysical predictions, these outliers can be smoking gun evidence for the
existence of PBHs as well. The next section would be dedicated to the
comparison with astrophysical populations.

We will first discuss the improvement on the measurement of the mass ratio
$1/q$, then the source-frame primary mass $m_{1}$, and finally the effective
spin $\chi_{\rm eff}$.

###  IV.1 Mass ratio

From the example of the $q-\iota-z$ joint measurement in Sec. III.5, we expect
that the presence of HoMs in the waveform model reduces the uncertainty in the
mass ratio measurement. Such improvement will also affect the measurement of
the primary mass, which will be explored in the next subsection.

(a)  (b)

Figure 7: Posteriors of mass ratio for sources with $(\hat{M}_{\rm
tot},\hat{\iota})=(40M_{\odot},30^{\circ})$ and $\hat{q}=1$ (top panel) or 2
(bottom panel) at $\hat{z}=10,20,30,40$ and 50, obtained by IMRPhenomPv2 (red)
and IMRPhenomXPHM (blue). Other plot settings are the same as in Fig. 2.

In Fig. 7, we compare the posteriors of mass ratio obtained by the two
waveform models for sources with $(\hat{M}_{\rm
tot},\hat{\iota})=(40M_{\odot},30^{\circ})$ and $\hat{q}=1$ (top panel) or 2
(bottom panel) at $\hat{z}=10,20,30,40$ and 50. Indeed, the uncertainties are
generally smaller in the cases of IMRPhenomXPHM. The reduction of the
uncertainties is more prominent when the redshift increases, because the
detector-frame masses are larger and more HoMs are detectable. In particular,
at $\hat{z}\geq 40$, the mass ratio is unconstrained in the cases of
IMRPhenomPv2, or even slightly biased away from $\hat{q}$. Similar to the
discussion in Sec. III.4, this is caused by the combined effect of the scaling
in the redshift prior and the correlation between $(1/q,z)$ as shown in the
IMRPhenomPv2 contours of Fig. 6. On the other hand, the posteriors in the
cases of IMRPhenomXPHM clearly show a peak around $\hat{q}$, i.e., consistent
with the true value.

###  IV.2 Primary mass

In GW astronomy, the mass parameters measured directly are the detector-frame
chirp mass (total mass) if the waveform is dominated by the inspiral (merger-
ringdown) phase within the sensitive frequency band. Converting from the
detector-frame chirp mass or total mass to the source-frame component masses,
there are a factor of $(1+z)$ and a Jacobian term as a function of $q$.
Therefore, the improvement on both measurements of $q$ and $z$ propagate to
the measurement of the source-frame primary mass.

In the upper panel of Fig. 8, we compare the posteriors of the source-frame
primary mass obtained by the two waveform models for sources with
$(\hat{m}_{1},\hat{m}_{2},\hat{\iota})=(20M_{\odot},20M_{\odot},30^{\circ})$
at $\hat{z}=10,20,30,40$ and 50. The uncertainty improves by $\sim 40\%$ at
$\hat{z}=10$ to $\sim 60\%$ at $\hat{z}=50$, owing to the presence of HoMs in
the waveform which reduces the uncertainty in both the redshift and the mass
ratio.

Figure 8: Posteriors of the source-frame primary mass for different sources
obtained by IMRPhenomPv2 (red) and IMRPhenomXPHM (blue). The parameters of the
sources are
$(\hat{m}_{1},\hat{m}_{2},\hat{\iota})=(20M_{\odot},20M_{\odot},30^{\circ})$
at $\hat{z}=10,20,30,40$ and 50 (upper panel),
$(\hat{m}_{1},\hat{m}_{2},\hat{\iota})=(80M_{\odot},80M_{\odot},30^{\circ})$
at $\hat{z}=10$ and 20 (lower left panel), and
$(\hat{m}_{1},\hat{m}_{2},\hat{\iota})=(125M_{\odot},125M_{\odot},30^{\circ})$
at $\hat{z}=10$ (lower right panel). Plot settings are the same as in Fig. 2.

As the possible PBH masses are not restricted within the stellar mass region,
we also show systems with a larger mass range, $\hat{m}_{1}=80$ and
$125M_{\odot}$, but at smaller redshifts in the lower panel of Fig. 8. The
posterior of $m_{1}$ for the systems with $\hat{m}_{1}=80M_{\odot}$ at
$\hat{z}=10$ (20) is well-constrained within a relative uncertainty of $\sim
10\%$ $(\sim 50\%)$ in the cases of IMRPhenomXPHM, while those in the cases of
IMRPhenomPv2 have $\geq 5$ times larger uncertainties. Similarly, for the
system with $m_{1}=125M_{\odot}$, the relative uncertainty of $m_{1}$ is
greatly reduced from $\sim 130\%$ (IMRPhenomPv2) to $\sim 15\%$
(IMRPhenomXPHM).

###  IV.3 Effective spin

In this subsection, we focus on the weakly-accreting scenario [25] and
demonstrate if the zero spin can be well measured, and whether the measurement
can be improved by the presence of HoMs in the waveform model. In Fig. 9, we
compare the posterior of $\chi_{\rm eff}$ obtained by the two waveform models
for zero-spin sources with $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota})=(40M_{\odot},1,30^{\circ})$ at $\hat{z}=10,20,30,40$
and 50. The uncertainties are comparable in both waveform models. We note that
the posteriors in the cases of IMRPhenomXPHM are more asymmetric around 0,
with larger probability masses in the negative $\chi_{\rm eff}$ region, unlike
the trends shown in Ref. [101]. Hence, we conclude that the presence of HoMs
in the waveform model only changes the morphology of the likelihood function
along $\chi_{\rm eff}$, and does not improve the measurement of the spin
parameters.

Figure 9: Posteriors of effective spin for zero-spin sources with
$(\hat{M}_{\rm tot},\hat{q},\hat{\iota})=(40M_{\odot},1,30^{\circ})$ at
$\hat{z}=10,20,30,40$ and 50, obtained by IMRPhenomPv2 (red) and IMRPhenomXPHM
(blue). Grey solid lines indicate the prior on $\chi_{\rm eff}$ sampled from
the uniform prior on spin magnitudes, spin orientations, and $1/q$. For the
sake of visualization, the areas of the violins are not normalized in common.
Thus, the grey lines are only normalized to their corresponding violins and
have different widths. Plot settings are the same as in Fig. 2.

We also note that it is generally harder to measure the spins of the systems
when the redshift increases. This is because the spin effect first goes into
the 1.5 PN inspiral phase, but the signal drifts to the lower frequency and
most of the SNR is dominated by the merger-ringdown phase instead.

##  V Implications for the PBH detection

In the previous section, we have established that including HoMs in the
waveform model improves the single measurement of the redshift and masses,
while the uncertainty of effective spin remains unchanged. In this section, we
discuss how the single-event measurement of the redshift, primary mass and
effective spin obtained by IMRPhenomXPHM can be used to identify or infer the
properties of PBHs.

###  V.1 Redshift measurement

The redshift may be one of the most significant parameter in the
identification of PBHs. This is because PBHs are formed much earlier than the
astrophysical BHs. While the time of birth of the first BH is still uncertain,
theoretical and simulation studies suggest that the astrophysical epoch of
BBHs is about $z\sim 30$ [66, 67, 68, 69, 102, 103, 104, 105, 106].

In Refs. [32, 107, 1], the critical redshift to approximate when the merger
rate density of Pop III BBHs turns off has been chosen to be
$z_{\mathrm{crit}}=30$, i.e., the redshift region is fully primordial above
$z_{\mathrm{crit}}$. With such definition, we showed that the relative
abundance of PBH and Pop III mergers can affect the significance of confirming
the primordial origin of a single BBH detection. The heuristic prior $p_{\rm
tot}$ based on the expected merger rate densities is then

$\displaystyle p_{\rm tot}$
$\displaystyle\left(z|f^{\mathrm{PBH}}_{\mathrm{III}}\right)\propto$
$\displaystyle\left[f^{\mathrm{PBH}}_{\mathrm{III}}\frac{\dot{n}_{\rm
PBH}(z)}{\dot{n}_{\rm PBH}(z_{\mathrm{crit}})}+\frac{\dot{n}_{\rm
III}(z)}{\dot{n}_{\rm
III}(z_{\mathrm{crit}})}\right]\frac{dV_{c}}{dz}\frac{1}{1+z},$ (1)

where

$\displaystyle\dot{n}_{\rm III}(z)\propto\begin{cases}\frac{e^{a_{\rm
III}(z-z_{\rm III})}}{b_{\rm III}+a_{\rm III}e^{(a_{\rm III}+b_{\rm
III})(z-z_{\rm III})}}&\mathrm{if}\,z<z_{\mathrm{crit}}\\\
0&\mathrm{otherwise}\end{cases},$ (2)

is the phenomenological fit to the Pop III merger rate density based on the
simulation study in Ref. [69, 97], with $(a_{\rm III},b_{\rm III},z_{\rm
III})=(0.66,0.3,11.6)$ from [97],

$\displaystyle\dot{n}_{\rm
PBH}(z)\propto\left(\frac{t(z)}{t(0)}\right)^{-34/37},$ (3)

is the analytic PBH merger rate density obtained from the dynamics of early
PBH binary formation [108, 109, 70, 110, 26] (with $t(z)$ being the age of the
Universe at $z$), and

$\displaystyle
f^{\mathrm{PBH}}_{\mathrm{III}}(z_{\mathrm{crit}})\equiv\frac{\dot{n}_{\rm
PBH}(z_{\mathrm{crit}})}{\dot{n}_{\rm III}(z_{\mathrm{crit}})}$

is the ratio between the two merger rate densities at $z_{\mathrm{crit}}$.

We also employed two statistical indicator in Ref. [1]. The first indicator is
the probability of primordial origin, $P_{\text{p}}$ as the fraction of the
redshift posterior with $z\geq z_{\mathrm{crit}}$,

$\displaystyle
P_{\text{p}}(z_{\mathrm{crit}}|f^{\mathrm{PBH}}_{\mathrm{III}})=\frac{1}{Z(f^{\mathrm{PBH}}_{\mathrm{III}})}\int_{z_{\mathrm{crit}}}^{\infty}\frac{p\left(z|d\right)}{p_{0}(z)}p_{\rm
tot}\left(z|f^{\mathrm{PBH}}_{\mathrm{III}}\right)dz,$ (4)

where $p\left(z|d\right)$ is the redshift posterior obtained with the default
prior $p_{0}(z)\propto dV_{c}/dz/(1+z)$ (uniform source-frame rate density),
and $Z(f^{\mathrm{PBH}}_{\mathrm{III}})$ is the evidence of the merger rate
density model $p_{\rm tot}$ parameterized by
$f^{\mathrm{PBH}}_{\mathrm{III}}$. We found that the best system in our
simulation set, $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota},\hat{z})=(40M_{\odot},1,60^{\circ},40)$, i.e., the
posterior mass has the most support at $z\geq z_{\mathrm{crit}}$, can only
provides mild evidence for its primordial origin. The second indicator is the
Bayes factor between models differed by their priors evaluated at different
$f^{\mathrm{PBH}}_{\mathrm{III}}$. As the estimation of the Bayes factor is
sensitive to the sampling algorithm and computationally expensive, we do not
repeat the statistical analysis using the Bayes factor. In the following, we
relax the assumption of $z_{\mathrm{crit}}=30$, and explore how $P_{\text{p}}$
varies as a function of $z_{\mathrm{crit}}\in[15,40]$.

First, we revisit the system $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota},\hat{z})=(40M_{\odot},1,60^{\circ},40)$ and show
$P_{\text{p}}(z_{\mathrm{crit}})$ evaluated at
$f^{\mathrm{PBH}}_{\mathrm{III}}=100,10,1,0.1$ and 0.01 in Fig. 10. A smaller
$f^{\mathrm{PBH}}_{\mathrm{III}}$ requires a lower $z_{\mathrm{crit}}$ to
maintain $P_{\text{p}}=0.9$ (black dotted line). This is expected since the
decrease in the prior volume of $z\geq z_{\mathrm{crit}}$ due to a smaller
$f^{\mathrm{PBH}}_{\mathrm{III}}$ can be compensated by lowering
$z_{\mathrm{crit}}$ for a larger primordial region. In this example, even if
$f^{\mathrm{PBH}}_{\mathrm{III}}=0.01$ at $z_{\mathrm{crit}}\sim 23$, the
redshift measurement still allows for $P_{\text{p}}=0.9$.

Figure 10: Probabilities of a primordial origin $P_{\text{p}}$ as a function
of $z_{\mathrm{crit}}$ for the system $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota},\hat{z})=(40M_{\odot},1,60^{\circ},40)$ evaluated at
$f^{\mathrm{PBH}}_{\mathrm{III}}=$100 (blue), 10 (orange), 1 (green), 0.1
(red) and 0.01 (purple). The calculations are based on the IMRPhenomXPHM
posteriors.

Next, we choose the optimistic scenario,
$f^{\mathrm{PBH}}_{\mathrm{III}}=100$, and illustrate how
$P_{\text{p}}(z_{\mathrm{crit}})$ differs for various sets of source
parameters. We pivot the system $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota},\hat{z})=(40M_{\odot},1,30^{\circ},30)$, and vary the
values of these four parameters one by one, as shown in each panel of Fig. 11.
The intersection of the black dotted line and each colored line indicates the
maximum $z_{\mathrm{crit}}$ to reach $P_{\text{p}}=0.9$ for each system. Most
of the systems have the values of maximum $z_{\mathrm{crit}}$ above 20, except
for the system at $\hat{z}=20$ which has a maximum $z_{\mathrm{crit}}\sim 18$.
The shift of maximum $z_{\mathrm{crit}}$ is a direct result of the change of
the redshift uncertainties due to the change of
$\hat{\iota},\hat{q},\hat{M}_{\rm tot}$ and $\hat{z}$ as shown in Figs. 2, 3,
4 and 5, respectively.

(a) $(\hat{M}_{\rm tot},\hat{z},\hat{q})=(40M_{\odot},30,1)$ (b)
$(\hat{M}_{\rm tot},\hat{\iota},\hat{z})=(40M_{\odot},30^{\circ},30)$ (c)
$(\hat{q},\hat{\iota},\hat{z})=(1,30^{\circ},30)$ (d) $(\hat{M}_{\rm
tot},\hat{q},\hat{\iota})=(40M_{\odot},1,30^{\circ})$

Figure 11: Probabilities of a primordial origin $P_{\text{p}}$ as a function
of $z_{\mathrm{crit}}$ for the systems with (a) $(\hat{M}_{\rm
tot},\hat{q},\hat{z})=(40M_{\odot},1,30)$ and $\iota=0^{\circ}$ (blue),
$30^{\circ}$ (orange), $60^{\circ}$ (green) and $90^{\circ}$ (red); (b)
$(\hat{M}_{\rm tot},\hat{\iota},\hat{z})=(40M_{\odot},30^{\circ},30)$ and
$\hat{q}=$1 (blue), 2 (orange), 3 (green) and 4 (red); (c)
$(\hat{q},\hat{\iota},\hat{z})=(1,30^{\circ},30)$ and $\hat{M}_{\rm tot}=$10
(blue), 20 (orange), 40 (green) and 80$M_{\odot}$ (red); as well as (d)
$(\hat{M}_{\rm tot},\hat{q},\hat{\iota})=(40M_{\odot},1,30^{\circ})$ and
$z=20$ (blue), 30 (orange), 40 (green) and 50 (red), all evaluated at
$f^{\mathrm{PBH}}_{\mathrm{III}}=100$. The calculations are based on the
IMRPhenomXPHM posteriors. The choice of $z_{\mathrm{crit}}=30$ is marked by
the black solid line.

###  V.2 Mass measurement

The improvement on the primary mass measurement as discussed in Sec. IV.2
certainly helps the statistical inference of the PBH mass spectrum (if they
exist within the detectable mass range). Roughly speaking, the uncertainties
of the population properties (which can be the median mass and variance of a
log-normal mass spectrum in the context of PBHs) scale with the uncertainties
of the single measurement, and inversely with the square root of the number of
detections. A factor of two improvement on the mass measurement (as seen in
Fig. 8) translates to a factor of four reduction of the required number of
events to reach the same statistical uncertainties of the population
properties.

Moreover, the improvement on the mass measurement may also increase the
statistical power for the identification of PBHs whose masses are outside the
astrophysically allowed values. Astrophysical BBHs originated from Pop I/II
stars are not expected to have masses lying in the pair-instability supernova
(PISN) mass gap $(\sim 50-130M_{\odot})$ [111, 112, 113], unless the mergers
consist of remnants of previous mergers in dense stellar clusters [114] or
gas-rich environment [59]. Recent studies of BBHs originated from Pop III
stars suggest that the mass gap may be narrowed to $\sim 100-130M_{\odot}$
[71]. In both cases, a good measurement of $m_{1}$ is valuable in both
measuring the mass spectra of Pop III and PBH mergers, or serves as a possible
indicator for the existence of PBHs provided that the Pop III mergers cannot
fill up the mass gap efficiently at high redshifts. Given the model
uncertainties in Pop III mergers, let us stick with the “standard” PISN mass
gap scenario: $\sim 50-130M_{\odot}$. In Fig. 8, the mass uncertainties of the
in-gap systems obtained by IMRPhenomXPHM lie within the mass gap, while those
obtained by IMRPhenomPv2 do not. This shows the importance of HoMs for
identifying outliers that are inconsistent with the astrophysical predictions,
such as the existence of PBHs.

###  V.3 Spin measurement

The spin spectrum of the Pop III BBHs is still highly uncertain, since both
the initial conditions of the Pop III stars and the binary evolution under
very low metallicity are yet to be understood. Pop III BBHs may have non-
negligible aligned spins, i.e., the majority of the BBHs have $\chi_{\rm
eff}\geq 0$ [71]. On the other hand, PBH binaries are dynamically formed. Even
if each component PBH acquires a significant spin through accretion, the
distribution of $\chi_{\rm eff}$ is expected to be symmetric around 0 [25].
This may be another signature to distinguish between Pop III and primordial
BBH population at high redshift. In fact, such method has been proposed to
measure the branching ratio between the formation channels of the dynamical
capture and the isolated binary evolution using the LIGO/Virgo BBHs within
$z<1$ [115, 116, 117].

As seen in Fig. 9, including HoMs does not improve the effective spin
measurement. Also, the typical uncertainties of these high-redshift BBHs are
comparable to those of the current detections [118, 27, 119]. Therefore, the
required number of detections needed for distingushing the formation channels
may be $\mathcal{O}(100)$, similar to that predicted in previous studies [116,
115].

##  VI Discussions

In this work, we have simulated BBHs merging at $z\geq 10$ with different mass
ratios, total masses, and inclinations. Performing full Bayesian parameter
estimation with a waveform model that includes HoMs and spin precession, we
quantified the uncertainties in redshifts, mass ratios, source-frame primary
masses, and effective spins, as measured by a network of next-generation
ground-based GW detectors. We compared the uncertainties obtained by a
waveform model without HoMs, and found that adding HoMs to the waveform model
can improve the measurements of redshifts and component masses by up to a
factor of 5. Such improvements originate from the different dependence of mass
ratio and inclination in the amplitude of each HoM, which helps breaking the
$(2,2)$ distance-inclination degeneracy and constraining the mass ratio.
Generally speaking, BBHs with more asymmetric mass ratios, more inclined
orbit, or larger masses in the detector frame are more benefited by the HoMs.
We note that the improvements are sensitive to the exact configuration of the
BBH, which needs to generate large enough amplitudes of HoMs in order to break
the degeneracy [120]. We also showed that including HoMs in the waveform
models has no effect on the spin measurements for zero-spin sources.

We note that waveform systematics due to different descriptions of the spin
precession dynamics are still present in the waveform families without HoMs,
namely IMRPhenomPv2 and IMRPhenomXP [80, 81]. However, the improvements shown
in this work are driven by the inclusion of HoMs, which breaks the distance-
inclination degeneracy (see Fig. 1), not by the detailed description of the
spin dynamics. Therefore, we expect our results would remain unchanged if the
systems have non-negligible spins. Given the improved uncertainties in
redshifts, we revisit the evaluation of primordial probabilities, which depend
on both the critical redshift that defines the epoch of astrophysical BBHs,
and the relative abundance characterized by the rate density ratio between Pop
III and PBH mergers at the critical redshift. With the best redshift
measurements in our simulation set, one may ensure the primordial probability
larger than $90\%$ for a relative abundance
$f^{\mathrm{PBH}}_{\mathrm{III}}(z_{\mathrm{crit}})\geq 0.01$, if the critical
redshift is $\gtrsim 20$. We also explored how BBHs with different parameters
react to the calculation of primordial probabilities. Assuming a large
relative abundance $f^{\mathrm{PBH}}_{\mathrm{III}}(z_{\mathrm{crit}})=100$,
the minimum $z_{\mathrm{crit}}$ to achieve a primordial probability of 90% is
$\gtrsim 20$, depending on the true redshifts of the systems mostly.

The improvements on the measurement of masses enable a better inference of the
mass spectra of both the Pop III and PBH mergers. While the measurements of
effective spins remain mostly unchanged, the uncertainties are comparable to
those of the current measurements. One may borrow the method of using the
distribution of effective spins to distinguish these two high-redshift BBHs
based on different formation scenarios as suggested by the existing
literature. In any case, we expect the mass and spin measurements to provide
additional evidence for identifying PBHs using high-redshift GW observations.
We will investigate this avenue as a future work.

##  VII Acknowledgments

We would like to thank Valerio De Luca for suggestions and comments. KKYN and
SV are supported by the NSF through the award PHY-1836814. SC is supported by
the Undergraduate Research Opportunities Program of Massachusetts Institute of
Technology. The work of MM is supported by the Swiss National Science
Foundation and by the SwissMap National Center for Competence in Research. MB
acknowledges support from the European Union’s Horizon 2020 Programme under
the AHEAD2020 project (grant agreement n. 871158). BSS is supported in part by
NSF Grant No. PHY-2207638, AST-2006384 and PHY- 2012083. SB is supported by
NSF Grant No.PHY-1836779. BG is supported by the Italian Ministry of
Education, University and Research within the PRIN 2017 Research Program
Framework, n. 2017SYRTCN. A.R. are supported by the Swiss National Science
Foundation (SNSF), project The Non-Gaussian Universe and Cosmological
Symmetries, project number: 200020-178787. G.F. acknowledges financial support
provided under the European Union’s H2020 ERC, Starting Grant agreement no.
DarkGRA–757480, and under the MIUR PRIN and FARE programmes (GW-NEXT, CUP:
B84I20000100001), and support from the Amaldi Research Center funded by the
MIUR program “Dipartimento di Eccellenza” (CUP: B81I18001170001) and
H2020-MSCA-RISE-2020 GRU.

## References

  * Ng  _et al._ [2022a] K. K. Y. Ng, S. Chen, B. Goncharov, U. Dupletsa, S. Borhanian, M. Branchesi, J. Harms, M. Maggiore, B. S. Sathyaprakash, and S. Vitale, On the Single-event-based Identification of Primordial Black Hole Mergers at Cosmological Distances, Astrophys. J. Lett. 931, L12 (2022a), arXiv:2108.07276 [astro-ph.CO] .
  * Zel’dovich and Novikov [1967] Y. B. Zel’dovich and I. D. Novikov, The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model, Soviet Astron. AJ (Engl. Transl. ), 10, 602 (1967).
  * Hawking [1974] S. W. Hawking, Black hole explosions, Nature 248, 30 (1974).
  * Chapline [1975] G. F. Chapline, Cosmological effects of primordial black holes, Nature 253, 251 (1975).
  * Carr [1975] B. J. Carr, The Primordial black hole mass spectrum, Astrophys. J. 201, 1 (1975).
  * Ivanov  _et al._ [1994] P. Ivanov, P. Naselsky, and I. Novikov, Inflation and primordial black holes as dark matter, Phys. Rev. D 50, 7173 (1994).
  * Garcia-Bellido  _et al._ [1996] J. Garcia-Bellido, A. D. Linde, and D. Wands, Density perturbations and black hole formation in hybrid inflation, Phys. Rev. D 54, 6040 (1996), arXiv:astro-ph/9605094 .
  * Ivanov [1998] P. Ivanov, Nonlinear metric perturbations and production of primordial black holes, Phys. Rev. D 57, 7145 (1998), arXiv:astro-ph/9708224 .
  * Blinnikov  _et al._ [2016] S. Blinnikov, A. Dolgov, N. K. Porayko, and K. Postnov, Solving puzzles of GW150914 by primordial black holes, JCAP 1611, 036, arXiv:1611.00541 [astro-ph.HE] .
  * Ali-Haïmoud [2018] Y. Ali-Haïmoud, Correlation Function of High-Threshold Regions and Application to the Initial Small-Scale Clustering of Primordial Black Holes, Phys. Rev. Lett. 121, 081304 (2018), arXiv:1805.05912 [astro-ph.CO] .
  * Desjacques and Riotto [2018] V. Desjacques and A. Riotto, Spatial clustering of primordial black holes, Phys. Rev. D 98, 123533 (2018), arXiv:1806.10414 [astro-ph.CO] .
  * Ballesteros  _et al._ [2018] G. Ballesteros, P. D. Serpico, and M. Taoso, On the merger rate of primordial black holes: effects of nearest neighbours distribution and clustering, JCAP 10, 043, arXiv:1807.02084 [astro-ph.CO] .
  * Moradinezhad Dizgah  _et al._ [2019] A. Moradinezhad Dizgah, G. Franciolini, and A. Riotto, Primordial Black Holes from Broad Spectra: Abundance and Clustering, JCAP 11, 001, arXiv:1906.08978 [astro-ph.CO] .
  * Inman and Ali-Haïmoud [2019] D. Inman and Y. Ali-Haïmoud, Early structure formation in primordial black hole cosmologies, Phys. Rev. D 100, 083528 (2019), arXiv:1907.08129 [astro-ph.CO] .
  * De Luca  _et al._ [2020a] V. De Luca, V. Desjacques, G. Franciolini, and A. Riotto, The clustering evolution of primordial black holes, JCAP 11, 028, arXiv:2009.04731 [astro-ph.CO] .
  * De Luca  _et al._ [2019] V. De Luca, V. Desjacques, G. Franciolini, A. Malhotra, and A. Riotto, The initial spin probability distribution of primordial black holes, JCAP 05, 018, arXiv:1903.01179 [astro-ph.CO] .
  * Mirbabayi  _et al._ [2020] M. Mirbabayi, A. Gruzinov, and J. Noreña, Spin of Primordial Black Holes, JCAP 03, 017, arXiv:1901.05963 [astro-ph.CO] .
  * Nakamura  _et al._ [1997] T. Nakamura, M. Sasaki, T. Tanaka, and K. S. Thorne, Gravitational waves from coalescing black hole MACHO binaries, Astrophys. J. Lett. 487, L139 (1997), arXiv:astro-ph/9708060 .
  * Ioka  _et al._ [1998] K. Ioka, T. Chiba, T. Tanaka, and T. Nakamura, Black hole binary formation in the expanding universe: Three body problem approximation, Phys. Rev. D 58, 063003 (1998), arXiv:astro-ph/9807018 .
  * Polnarev and Khlopov [1985] A. G. Polnarev and M. Y. Khlopov, COSMOLOGY, PRIMORDIAL BLACK HOLES, AND SUPERMASSIVE PARTICLES, Sov. Phys. Usp. 28, 213 (1985).
  * Khlopov [2010] M. Y. Khlopov, Primordial Black Holes, Res. Astron. Astrophys. 10, 495 (2010), arXiv:0801.0116 [astro-ph] .
  * Sasaki  _et al._ [2018] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Primordial black holes—perspectives in gravitational wave astronomy, Class. Quant. Grav. 35, 063001 (2018), arXiv:1801.05235 [astro-ph.CO] .
  * Green and Kavanagh [2021] A. M. Green and B. J. Kavanagh, Primordial Black Holes as a dark matter candidate, J. Phys. G 48, 043001 (2021), arXiv:2007.10722 [astro-ph.CO] .
  * Franciolini [2021] G. Franciolini,  _Primordial Black Holes: from Theory to Gravitational Wave Observations_ , Ph.D. thesis, Geneva U., Dept. Theor. Phys. (2021), arXiv:2110.06815 [astro-ph.CO] .
  * De Luca  _et al._ [2020b] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, The evolution of primordial black holes and their final observable spins, JCAP 04, 052, arXiv:2003.02778 [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, arXiv:2005.05641 [astro-ph.CO] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * De Luca  _et al._ [2021a] V. De Luca, V. Desjacques, G. Franciolini, P. Pani, and A. Riotto, GW190521 Mass Gap Event and the Primordial Black Hole Scenario, Phys. Rev. Lett. 126, 051101 (2021a), arXiv:2009.01728 [astro-ph.CO] .
  * Wong  _et al._ [2021] K. W. K. Wong, G. Franciolini, V. De Luca, V. Baibhav, E. Berti, P. Pani, and A. 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] .
  * Hütsi  _et al._ [2021] G. Hütsi, M. Raidal, V. Vaskonen, and H. Veermäe, Two populations of LIGO-Virgo black holes, JCAP 03, 068, arXiv:2012.02786 [astro-ph.CO] .
  * Hall  _et al._ [2020] A. Hall, A. D. Gow, and C. 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] .
  * De Luca  _et al._ [2021b] 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, arXiv:2102.03809 [astro-ph.CO] .
  * Franciolini  _et al._ [2022] G. Franciolini, V. Baibhav, V. De Luca, K. K. Y. Ng, K. W. K. Wong, E. Berti, P. Pani, A. Riotto, and S. Vitale, Searching for a subpopulation of primordial black holes in LIGO-Virgo gravitational-wave data, Phys. Rev. D 105, 083526 (2022), arXiv:2105.03349 [gr-qc] .
  * Mukherjee and Silk [2021] S. Mukherjee and J. Silk, Can we distinguish astrophysical from primordial black holes via the stochastic gravitational wave background?, Mon. Not. Roy. Astron. Soc. 506, 3977 (2021), arXiv:2105.11139 [gr-qc] .
  * Franciolini  _et al._ [2022] G. Franciolini, I. Musco, P. Pani, and A. Urbano, From inflation to black hole mergers and back again: Gravitational-wave data-driven constraints on inflationary scenarios with a first-principle model of primordial black holes across the QCD epoch, arXiv e-prints , arXiv:2209.05959 (2022), arXiv:2209.05959 [astro-ph.CO] .
  * Franciolini  _et al._ [2022] G. Franciolini, R. Cotesta, N. Loutrel, E. Berti, P. Pani, and A. Riotto, How to assess the primordial origin of single gravitational-wave events with mass, spin, eccentricity, and deformability measurements, Phys. Rev. D 105, 063510 (2022), arXiv:2112.10660 [astro-ph.CO] .
  * Hall and Evans [2019] E. D. Hall and M. Evans, Metrics for next-generation gravitational-wave detectors, Class. Quant. Grav. 36, 225002 (2019), arXiv:1902.09485 [astro-ph.IM] .
  * O’Shaughnessy  _et al._ [2017] R. O’Shaughnessy, J. Bellovary, A. Brooks, S. Shen, F. Governato, and C. Christensen, The effects of host galaxy properties on merging compact binaries detectable by LIGO, Mon. Not. Roy. Astron. Soc. 464, 2831 (2017), arXiv:1609.06715 [astro-ph.GA] .
  * Dominik  _et al._ [2012] M. Dominik, K. Belczynski, C. Fryer, D. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Double Compact Objects I: The Significance of the Common Envelope on Merger Rates, Astrophys. J. 759, 52 (2012), arXiv:1202.4901 [astro-ph.HE] .
  * Dominik  _et al._ [2013] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Double Compact Objects II: Cosmological Merger Rates, Astrophys. J. 779, 72 (2013), arXiv:1308.1546 [astro-ph.HE] .
  * Dominik  _et al._ [2015] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik, and F. Pannarale, Double Compact Objects III: Gravitational Wave Detection Rates, Astrophys. J. 806, 263 (2015), arXiv:1405.7016 [astro-ph.HE] .
  * de Mink and Belczynski [2015] S. de Mink and K. Belczynski, Merger rates of double neutron stars and stellar origin black holes: The Impact of Initial Conditions on Binary Evolution Predictions, Astrophys. J. 814, 58 (2015), arXiv:1506.03573 [astro-ph.HE] .
  * Belczynski  _et al._ [2016a] K. Belczynski, D. E. Holz, T. Bulik, and R. O’Shaughnessy, The first gravitational-wave source from the isolated evolution of two 40-100 Msun stars, Nature 534, 512 (2016a), arXiv:1602.04531 [astro-ph.HE] .
  * Stevenson  _et al._ [2017] S. Stevenson, A. Vigna-Gómez, I. Mandel, J. W. Barrett, C. J. Neijssel, D. Perkins, and S. E. de Mink, Formation of the first three gravitational-wave observations through isolated binary evolution, Nature Commun. 8, 14906 (2017), arXiv:1704.01352 [astro-ph.HE] .
  * Mapelli  _et al._ [2019] M. Mapelli, N. Giacobbo, F. Santoliquido, and M. C. Artale, The properties of merging black holes and neutron stars across cosmic time, Mon. Not. Roy. Astron. Soc. 487, 2 (2019), arXiv:1902.01419 [astro-ph.HE] .
  * Breivik  _et al._ [2020] K. Breivik _et al._ , COSMIC Variance in Binary Population Synthesis, Astrophys. J. 898, 71 (2020), arXiv:1911.00903 [astro-ph.HE] .
  * Bavera  _et al._ [2020] S. S. Bavera, T. Fragos, Y. Qin, E. Zapartas, C. J. Neijssel, I. Mandel, A. Batta, S. M. Gaebel, C. Kimball, and S. Stevenson, The origin of spin in binary black holes: Predicting the distributions of the main observables of Advanced LIGO, Astron. Astrophys. 635, A97 (2020), arXiv:1906.12257 [astro-ph.HE] .
  * Broekgaarden  _et al._ [2019] F. S. Broekgaarden, S. Justham, S. E. de Mink, J. Gair, I. Mandel, S. Stevenson, J. W. Barrett, A. Vigna-Gómez, and C. J. Neijssel, STROOPWAFEL: Simulating rare outcomes from astrophysical populations, with application to gravitational-wave sources, Mon. Not. Roy. Astron. Soc. 490, 5228 (2019), arXiv:1905.00910 [astro-ph.HE] .
  * Portegies Zwart and McMillan [2000] S. F. Portegies Zwart and S. L. W. McMillan, Black Hole Mergers in the Universe, ApJ 528, L17 (2000), arXiv:astro-ph/9910061 [astro-ph] .
  * Antonini and Gieles [2020] F. Antonini and M. Gieles, Merger rate of black hole binaries from globular clusters: theoretical error bars and comparison to gravitational wave data from GWTC-2, Phys. Rev. D 102, 123016 (2020), arXiv:2009.01861 [astro-ph.HE] .
  * Santoliquido  _et al._ [2021] F. Santoliquido, M. Mapelli, N. Giacobbo, Y. Bouffanais, and M. C. Artale, The cosmic merger rate density of compact objects: impact of star formation, metallicity, initial mass function and binary evolution, Mon. Not. Roy. Astron. Soc. 502, 4877 (2021), arXiv:2009.03911 [astro-ph.HE] .
  * Rodriguez  _et al._ [2015] C. L. Rodriguez, M. Morscher, B. Pattabiraman, S. Chatterjee, C.-J. Haster, and F. A. Rasio, Binary Black Hole Mergers from Globular Clusters: Implications for Advanced LIGO, Phys. Rev. Lett. 115, 051101 (2015), [Erratum: Phys.Rev.Lett. 116, 029901 (2016)], arXiv:1505.00792 [astro-ph.HE] .
  * Rodriguez  _et al._ [2016] C. L. Rodriguez, S. Chatterjee, and F. A. Rasio, Binary Black Hole Mergers from Globular Clusters: Masses, Merger Rates, and the Impact of Stellar Evolution, Phys. Rev. D 93, 084029 (2016), arXiv:1602.02444 [astro-ph.HE] .
  * Rodriguez and Loeb [2018] C. L. Rodriguez and A. Loeb, Redshift Evolution of the Black Hole Merger Rate from Globular Clusters, Astrophys. J. Lett. 866, L5 (2018), arXiv:1809.01152 [astro-ph.HE] .
  * Di Carlo  _et al._ [2019] U. N. Di Carlo, N. Giacobbo, M. Mapelli, M. Pasquato, M. Spera, L. Wang, and F. Haardt, Merging black holes in young star clusters, Mon. Not. Roy. Astron. Soc. 487, 2947 (2019), arXiv:1901.00863 [astro-ph.HE] .
  * Kremer  _et al._ [2020] K. Kremer, M. Spera, D. Becker, S. Chatterjee, U. N. Di Carlo, G. Fragione, C. L. Rodriguez, C. S. Ye, and F. A. Rasio, Populating the upper black hole mass gap through stellar collisions in young star clusters, Astrophys. J. 903, 45 (2020), arXiv:2006.10771 [astro-ph.HE] .
  * Bartos  _et al._ [2017] I. Bartos, B. Kocsis, Z. Haiman, and S. Márka, Rapid and Bright Stellar-mass Binary Black Hole Mergers in Active Galactic Nuclei, Astrophys. J. 835, 165 (2017), arXiv:1602.03831 [astro-ph.HE] .
  * Yi and Cheng [2019] S.-X. Yi and K. Cheng, Where Are the Electromagnetic-wave Counterparts of Stellar-mass Binary Black Hole Mergers?, Astrophys. J. Lett. 884, L12 (2019), arXiv:1909.08384 [astro-ph.HE] .
  * Yang  _et al._ [2019] Y. Yang _et al._ , Hierarchical Black Hole Mergers in Active Galactic Nuclei, Phys. Rev. Lett. 123, 181101 (2019), arXiv:1906.09281 [astro-ph.HE] .
  * Yang  _et al._ [2020] Y. Yang, I. Bartos, Z. Haiman, B. Kocsis, S. Márka, and H. Tagawa, Cosmic Evolution of Stellar-mass Black Hole Merger Rate in Active Galactic Nuclei, Astrophys. J. 896, 138 (2020), arXiv:2003.08564 [astro-ph.HE] .
  * Gröbner  _et al._ [2020] M. Gröbner, W. Ishibashi, S. Tiwari, M. Haney, and P. Jetzer, Binary black hole mergers in AGN accretion discs: gravitational wave rate density estimates, Astron. Astrophys. 638, A119 (2020), arXiv:2005.03571 [astro-ph.GA] .
  * Tagawa  _et al._ [2020a] H. Tagawa, Z. Haiman, and B. Kocsis, Formation and Evolution of Compact Object Binaries in AGN Disks, Astrophys. J. 898, 25 (2020a), arXiv:1912.08218 [astro-ph.GA] .
  * Tagawa  _et al._ [2021] H. Tagawa, B. Kocsis, Z. Haiman, I. Bartos, K. Omukai, and J. Samsing, Mass-gap Mergers in Active Galactic Nuclei, Astrophys. J. 908, 194 (2021), arXiv:2012.00011 [astro-ph.HE] .
  * Tagawa  _et al._ [2020b] H. Tagawa, Z. Haiman, I. Bartos, and B. Kocsis, Spin Evolution of Stellar-mass Black Hole Binaries in Active Galactic Nuclei, Astrophys. J. 899, 26 (2020b), arXiv:2004.11914 [astro-ph.HE] .
  * Samsing  _et al._ [2020] J. Samsing, I. Bartos, D. D’Orazio, Z. Haiman, B. Kocsis, N. Leigh, B. Liu, M. Pessah, and H. Tagawa, Active Galactic Nuclei as Factories for Eccentric Black Hole Mergers, Arxiv e-print (2020), arXiv:2010.09765 [astro-ph.HE] .
  * Kinugawa  _et al._ [2014] T. Kinugawa, K. Inayoshi, K. Hotokezaka, D. Nakauchi, and T. Nakamura, Possible Indirect Confirmation of the Existence of Pop III Massive Stars by Gravitational Wave, Mon. Not. Roy. Astron. Soc. 442, 2963 (2014), arXiv:1402.6672 [astro-ph.HE] .
  * Kinugawa  _et al._ [2016] T. Kinugawa, A. Miyamoto, N. Kanda, and T. Nakamura, The detection rate of inspiral and quasi-normal modes of Population III binary black holes which can confirm or refute the general relativity in the strong gravity region, Mon. Not. Roy. Astron. Soc. 456, 1093 (2016), arXiv:1505.06962 [astro-ph.SR] .
  * Hartwig  _et al._ [2016] T. Hartwig, M. Volonteri, V. Bromm, R. S. Klessen, E. Barausse, M. Magg, and A. Stacy, Gravitational Waves from the Remnants of the First Stars, Mon. Not. Roy. Astron. Soc. 460, L74 (2016), arXiv:1603.05655 [astro-ph.GA] .
  * Belczynski  _et al._ [2017] K. Belczynski, T. Ryu, R. Perna, E. Berti, T. Tanaka, and T. Bulik, On the likelihood of detecting gravitational waves from Population III compact object binaries, Mon. Not. Roy. Astron. Soc. 471, 4702 (2017), arXiv:1612.01524 [astro-ph.HE] .
  * Raidal  _et al._ [2019] M. Raidal, C. Spethmann, V. Vaskonen, and H. Veermäe, Formation and Evolution of Primordial Black Hole Binaries in the Early Universe, JCAP 02, 018, arXiv:1812.01930 [astro-ph.CO] .
  * Tanikawa  _et al._ [2022] A. Tanikawa, T. Yoshida, T. Kinugawa, A. A. Trani, T. Hosokawa, H. Susa, and K. Omukai, Merger Rate Density of Binary Black Holes through Isolated Population I, II, III and Extremely Metal-poor Binary Star Evolution, Astrophys. J. 926, 83 (2022), arXiv:2110.10846 [astro-ph.HE] .
  * Hijikawa  _et al._ [2021] K. Hijikawa, A. Tanikawa, T. Kinugawa, T. Yoshida, and H. Umeda, On the population III binary black hole mergers beyond the pair-instability mass gap, Mon. Not. Roy. Astron. Soc. 505, L69 (2021), arXiv:2104.13384 [astro-ph.HE] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific), Exploring the Sensitivity of Next Generation Gravitational Wave Detectors, Class. Quant. Grav. 34, 044001 (2017), arXiv:1607.08697 [astro-ph.IM] .
  * Reitze  _et al._ [2019] D. Reitze _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Evans  _et al._ [2021] M. Evans _et al._ , A horizon study for cosmic explorer science, observatories, and community, dcc.cosmicexplorer.org/CE-P2100003/public (2021).
  * 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] M. Maggiore _et al._ , Science Case for the Einstein Telescope, JCAP 03, 050, arXiv:1912.02622 [astro-ph.CO] .
  * Iacovelli  _et al._ [2022] F. Iacovelli, M. Mancarella, S. Foffa, and M. Maggiore, Forecasting the detection capabilities of third-generation gravitational-wave detectors using GWFAST, arXiv e-prints (2022), arXiv:2207.02771 [gr-qc] .
  * Pratten  _et al._ [2021] G. Pratten _et al._ , Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes, Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Pratten  _et al._ [2020] G. Pratten, S. Husa, C. Garcia-Quiros, M. Colleoni, A. Ramos-Buades, H. Estelles, and R. Jaume, Setting the cornerstone for a family of models for gravitational waves from compact binaries: The dominant harmonic for nonprecessing quasicircular black holes, Phys. Rev. D 102, 064001 (2020), arXiv:2001.11412 [gr-qc] .
  * García-Quirós  _et al._ [2020] C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, Multimode frequency-domain model for the gravitational wave signal from nonprecessing black-hole binaries, Phys. Rev. D 102, 064002 (2020), arXiv:2001.10914 [gr-qc] .
  * Hannam  _et al._ [2014] M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Husa  _et al._ [2016] S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. I. New numerical waveforms and anatomy of the signal, Phys. Rev. D 93, 044006 (2016), arXiv:1508.07250 [gr-qc] .
  * Khan  _et al._ [2016] S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Frequency-domain gravitational waves from nonprecessing black-hole binaries. II. A phenomenological model for the advanced detector era, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Bianchi  _et al._ [2018] E. Bianchi, A. Gupta, H. M. Haggard, and B. S. Sathyaprakash, Small Spins of Primordial Black Holes from Random Geometries: Bekenstein-Hawking Entropy and Gravitational Wave Observations, arXiv (2018), arXiv:1812.05127 [gr-qc] .
  * De Luca  _et al._ [2020d] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, Constraints on Primordial Black Holes: the Importance of Accretion, Phys. Rev. D 102, 043505 (2020d), arXiv:2003.12589 [astro-ph.CO] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Skilling [2006] J. Skilling, Nested sampling for general Bayesian computation, Bayesian Analysis 1, 833 (2006).
  * Speagle [2020] J. S. Speagle, DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences, MNRAS 493, 3132 (2020), arXiv:1904.02180 [astro-ph.IM] .
  * Ashton  _et al._ [2019] G. Ashton _et al._ , BILBY: A user-friendly Bayesian inference library for gravitational-wave astronomy, Astrophys. J. Suppl. 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Vallisneri [2008] M. Vallisneri, Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects, Phys. Rev. D 77, 042001 (2008), arXiv:gr-qc/0703086 .
  * Rodriguez  _et al._ [2014] C. L. Rodriguez, B. Farr, V. Raymond, W. M. Farr, T. B. Littenberg, D. Fazi, and V. Kalogera, Basic Parameter Estimation of Binary Neutron Star Systems by the Advanced LIGO/Virgo Network, Astrophys. J. 784, 119 (2014), arXiv:1309.3273 [astro-ph.HE] .
  * Pürrer and Haster [2020] M. Pürrer and C.-J. Haster, Gravitational waveform accuracy requirements for future ground-based detectors, Phys. Rev. Res. 2, 023151 (2020), arXiv:1912.10055 [gr-qc] .
  * Usman  _et al._ [2019] S. A. Usman, J. C. Mills, and S. Fairhurst, Constraining the Inclinations of Binary Mergers from Gravitational-wave Observations, Astrophys. J. 877, 82 (2019), arXiv:1809.10727 [gr-qc] .
  * Chen  _et al._ [2019] H.-Y. Chen, S. Vitale, and R. Narayan, Viewing angle of binary neutron star mergers, Phys. Rev. X 9, 031028 (2019), arXiv:1807.05226 [astro-ph.HE] .
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * Ng  _et al._ [2021] K. K. Y. Ng, S. Vitale, W. M. Farr, and C. L. Rodriguez, Probing multiple populations of compact binaries with third-generation gravitational-wave detectors, Astrophys. J. Lett. 913, L5 (2021), arXiv:2012.09876 [astro-ph.CO] .
  * Ng  _et al._ [2022b] K. K. Y. Ng, G. Franciolini, E. Berti, P. Pani, A. Riotto, and S. Vitale, Constraining High-redshift Stellar-mass Primordial Black Holes with Next-generation Ground-based Gravitational-wave Detectors, Astrophys. J. Lett. 933, L41 (2022b), arXiv:2204.11864 [astro-ph.CO] .
  * Blanchet  _et al._ [2008] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, The Third post-Newtonian gravitational wave polarisations and associated spherical harmonic modes for inspiralling compact binaries in quasi-circular orbits, Class. Quant. Grav. 25, 165003 (2008), [Erratum: Class.Quant.Grav. 29, 239501 (2012)], arXiv:0802.1249 [gr-qc] .
  * Borhanian  _et al._ [2020] S. Borhanian, K. G. Arun, H. P. Pfeiffer, and B. S. Sathyaprakash, Comparison of post-Newtonian mode amplitudes with numerical relativity simulations of binary black holes, Class. Quant. Grav. 37, 065006 (2020), arXiv:1901.08516 [gr-qc] .
  * Ng  _et al._ [2018] K. K. Y. Ng, S. Vitale, A. Zimmerman, K. Chatziioannou, D. Gerosa, and C.-J. Haster, Gravitational-wave astrophysics with effective-spin measurements: asymmetries and selection biases, Phys. Rev. D98, 083007 (2018), arXiv:1805.03046 [gr-qc] .
  * Inayoshi  _et al._ [2017] K. Inayoshi, R. Hirai, T. Kinugawa, and K. Hotokezaka, Formation pathway of Population III coalescing binary black holes through stable mass transfer, Mon. Not. Roy. Astron. Soc. 468, 5020 (2017), arXiv:1701.04823 [astro-ph.HE] .
  * Liu and Bromm [2020a] B. Liu and V. Bromm, The Population III origin of GW190521, Astrophys. J. Lett. 903, L40 (2020a), arXiv:2009.11447 [astro-ph.GA] .
  * Liu and Bromm [2020b] B. Liu and V. Bromm, Gravitational waves from Population III binary black holes formed by dynamical capture, Mon. Not. Roy. Astron. Soc. 495, 2475 (2020b), arXiv:2003.00065 [astro-ph.CO] .
  * Kinugawa  _et al._ [2020] T. Kinugawa, T. Nakamura, and H. Nakano, Chirp Mass and Spin of Binary Black Holes from First Star Remnants, Mon. Not. Roy. Astron. Soc. 498, 3946 (2020), arXiv:2005.09795 [astro-ph.HE] .
  * Tanikawa  _et al._ [2021] A. Tanikawa, H. Susa, T. Yoshida, A. A. Trani, and T. Kinugawa, Merger rate density of Population III binary black holes below, above, and in the pair-instability mass gap, Astrophys. J. 910, 30 (2021), arXiv:2008.01890 [astro-ph.HE] .
  * De Luca  _et al._ [2021c] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, The minimum testable abundance of primordial black holes at future gravitational-wave detectors, JCAP 11, 039, arXiv:2106.13769 [astro-ph.CO] .
  * Raidal  _et al._ [2017] M. Raidal, V. Vaskonen, and H. Veermäe, Gravitational Waves from Primordial Black Hole Mergers, JCAP 09, 037, arXiv:1707.01480 [astro-ph.CO] .
  * Chen and Huang [2018] Z.-C. Chen and Q.-G. Huang, Merger Rate Distribution of Primordial-Black-Hole Binaries, Astrophys. J. 864, 61 (2018), arXiv:1801.10327 [astro-ph.CO] .
  * Chen and Huang [2020] Z.-C. Chen and Q.-G. Huang, Distinguishing Primordial Black Holes from Astrophysical Black Holes by Einstein Telescope and Cosmic Explorer, JCAP 08, 039, arXiv:1904.02396 [astro-ph.CO] .
  * Belczynski  _et al._ [2016b] K. Belczynski _et al._ , The Effect of Pair-Instability Mass Loss on Black Hole Mergers, Astron. Astrophys. 594, A97 (2016b), arXiv:1607.03116 [astro-ph.HE] .
  * Spera and Mapelli [2017] M. Spera and M. Mapelli, Very massive stars, pair-instability supernovae and intermediate-mass black holes with the SEVN code, Mon. Not. Roy. Astron. Soc. 470, 4739 (2017), arXiv:1706.06109 [astro-ph.SR] .
  * Woosley [2017] S. E. Woosley, Pulsational Pair-Instability Supernovae, Astrophys. J. 836, 244 (2017), arXiv:1608.08939 [astro-ph.HE] .
  * Rodriguez  _et al._ [2018] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, and F. A. Rasio, Post-Newtonian Dynamics in Dense Star Clusters: Highly-Eccentric, Highly-Spinning, and Repeated Binary Black Hole Mergers, Phys. Rev. Lett. 120, 151101 (2018), arXiv:1712.04937 [astro-ph.HE] .
  * Vitale  _et al._ [2017] S. Vitale, R. Lynch, R. Sturani, and P. Graff, Use of gravitational waves to probe the formation channels of compact binaries, Class. Quant. Grav. 34, 03LT01 (2017), arXiv:1503.04307 [gr-qc] .
  * Farr  _et al._ [2018] B. Farr, D. E. Holz, and W. M. Farr, Using Spin to Understand the Formation of LIGO and Virgo’s Black Holes, Astrophys. J. Lett. 854, L9 (2018), arXiv:1709.07896 [astro-ph.HE] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, Virgo), Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog, Astrophys. J. Lett. 913, L7 (2021b), arXiv:2010.14533 [astro-ph.HE] .
  * Abbott  _et al._ [2019] B. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, ArXiv e-print (2021c), arXiv:2111.03606 [gr-qc] .
  * Mills and Fairhurst [2021] C. Mills and S. Fairhurst, Measuring gravitational-wave higher-order multipoles, Phys. Rev. D 103, 024042 (2021), arXiv:2007.04313 [gr-qc] .