# Identifying heavy stellar black holes at cosmological distances
with next generation gravitational-wave observatories

Stephen Fairhurst  Gravity Exploration Institute, School of Physics and
Astronomy, Cardiff University, Cardiff, UK, CF24 3AA, UK  Cameron Mills
Albert-Einstein-Institut, Max-Planck-Institut for Gravitationsphysik, D-30167
Hannover, Germany  Leibniz Universitat Hannover, D-30167 Hannover, Germany
Monica Colpi  Department of Physics, University of Milano - Bicocca, Piazza
della Scienza 3,I20126 Milano, Italy  National Institute of Nuclear Physics
INFN, Milano - Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  Raffaella
Schneider  Dipartimento di Fisica, Universitá di Roma ‘La Sapienza’, P.le Aldo
Moro 2, I-00185 Roma, Italy  INAF-Osservatorio Astronomico di Roma, via di
Frascati 33, I-00078 Monteporzio Catone, Italy  INFN, Sezione di Roma I, P.le
Aldo Moro 2, I-00185 Roma, Italy  Alberto Sesana  Department of Physics,
University of Milano - Bicocca, Piazza della Scienza 3,I20126 Milano, Italy
National Institute of Nuclear Physics INFN, Milano - Bicocca, Piazza della
Scienza 3, 20126 Milano, Italy  Alessandro Trinca  INAF-Osservatorio
Astronomico di Roma, via di Frascati 33, I-00078 Monteporzio Catone, Italy
INFN, Sezione di Roma I, P.le Aldo Moro 2, I-00185 Roma, Italy  Rosa Valiante
INAF-Osservatorio Astronomico di Roma, via di Frascati 33, I-00078 Monteporzio
Catone, Italy  INFN, Sezione di Roma I, P.le Aldo Moro 2, I-00185 Roma, Italy

###### Abstract

We investigate the detectability of single-event coalescing black hole
binaries with total mass of $100-600\,\rm M_{\odot}$ at cosmological distances
($5\lesssim z\lesssim 20$) with the next generation of terrestrial
gravitational wave observatories, specifically Einstein Telescope and Cosmic
Explorer. Our ability to observe these binaries is limited by the low-
frequency performance of the detectors. Higher-order Multipoles of the
gravitational wave signal are observable in these systems, and detection of
such multipoles serves to both b the mass range over which black hole binaries
are observable and improve the recovery of their individual masses and
redshift. For high redshift systems of $\sim 200\,\rm M_{\odot}$ we will be
able to confidently infer that the redshift is at least $z=12$, and for
systems of $\sim 400\,\rm M_{\odot}$ we can infer a minimum redshift of at
least $z=8$. We discuss the impact that these observations will have in
narrowing uncertainties on the existence of the pair-instability mass-gap, and
their implications on the formation of the first stellar black holes that
could be seeds for the growth of supermassive black holes powering high-$z$
quasars.

##  1 Introduction

In the three observing runs by the LIGO-Virgo-KAGRA (LVK) Collaboration, many
tens of gravitational-wave (GW) transient signals consistent with the merger
of binary black holes have been detected (Abbott et al., 2019; The LIGO
Scientific Collaboration et al., 2021a, b). The majority of them, observed at
$z\lesssim 1$, are found to have one or both component masses between $20\,\rm
M_{\odot}$ and $50\,\rm M_{\odot}$ and total mass $\lesssim 80\,\rm
M_{\odot}$.111All masses in this paper are referred to as measured in the
source-frame. A variety of channels have been proposed for their origin:
formation as field binaries primarily in low-metallicity galaxies at high
redshifts, formation in dense stellar systems, in active galactic nucleus
(AGN) discs, or after generations of repeated mergers (see Graziani et al.,
2020; Mapelli, 2021; Gerosa & Fishbach, 2021; Mandel & Farmer, 2022, and
references therein).

The next-generation of ground-based GW observatories, specifically the
Einstein Telescope (ET, Punturo et al., 2014; Branchesi et al., 2023) and
Cosmic Explorer (CE, Evans et al., 2021), will open the prospect of detecting
the GW signatures of merging BBHs over a wider mass range and deeper
redshifts, extending the realm of observations to BBHs out to $z\sim 30$, when
the first stars began to shine, and into the intermediate-mass range ${\cal
O}(100-1000)\,\rm M_{\odot}$ (Kalogera et al., 2019; Maggiore et al., 2020).
Beyond redshift $z\sim 30-40$, merging primordial black holes of ${\cal
O}(10)\,\rm M_{\odot}$, formed by quantum processes in the early Universe
(Carr et al., 2021; Ng et al., 2022), may also be detected and studied (Ng et
al., 2022, 2023).

In this paper we study systems with individual masses extending from 60 to
480$\,\rm M_{\odot}$, and total masses of $180\,\rm M_{\odot}$ to $600\,\rm
M_{\odot}$, covering a mass range that is relevant for several reasons, as we
explore below. Measuring such systems is most interesting at cosmological
distances, which is only possible due to the enhanced sensitivity of ET and CE
at frequencies below $10$ Hz.

The formation of such heavy stellar-mass black holes requires the presence of
star forming regions of zero or extremely low metallicity, where fragmentation
and cooling of the parent gas cloud, and mass loss from stellar winds are
strongly suppressed (see Bromm, 2013; Spera & Mapelli, 2017, and references
therein). These are conditions that occur in the high-redshift Universe, and
are expected to result in a top-heavy mass function where stars heavier than
$150-300\,\rm M_{\odot}$ are more common than in the conventionally adopted,
Kroupa-like, stellar initial mass function (IMF) (see Klessen & Glover 2023
for a recent review on the first stars). At the highest redshifts, these heavy
BHs may represent systems not yet affected by accretion of surrounding gas
(Roupas & Kazanas, 2019; Safarzadeh & Haiman, 2020; van Son et al., 2020), and
hence their masses reflect their birth conditions. Detecting GWs from these
heavy stellar-mass binaries will let us constrain their merger rate which is
intimately related to the rate of formation of massive stellar binary systems
in pristine star forming galaxies (Stacy & Bromm, 2013; Chon & Hosokawa, 2019;
Sugimura et al., 2023).

Some of the BBH masses we investigate reside within the so called pair-
instability mass-gap (often referred to as upper-mass gap or Pair Instability
Supernova (PISN) gap. This gap is between about $65\,\rm M_{\odot}$ and
$135\,\rm M_{\odot}$ where no BH is expected to form in evolution models of
isolated stars. This mass-gap is attributed to a pair-instability, arising in
metal poor, massive stars between about $130\,\rm M_{\odot}$ and $250\,\rm
M_{\odot},$ which leads to a supernova explosion due to uncontrolled
${}^{12}{\rm C}(\alpha,\gamma)^{16}{\rm O}$ nuclear burning, leaving no
remnant (Farmer et al., 2020; Woosley & Heger, 2021). During the third
observing run of LVK a short duration signal, GW190521, was detected and
estimated to be consistent with the merger of two BHs with component masses of
about $85\,\rm M_{\odot}$ and $66\,\rm M_{\odot}$ and total mass of $142\,\rm
M_{\odot}$ (Abbott et al., 2020a, b). This is the heaviest BH observed in GWs
to date, with an intermediate-mass remnant and a primary component residing
within the pair-instability mass-gap.222Alternative analyses, (e.g. Nitz &
Capano, 2021) find that GW190521 could instead be $\sim 170\,\rm M_{\odot}$
black hole with a companion of $\sim 20\,\rm M_{\odot}$, suggesting that the
primary is already an intermediate-mass BH, with a mass beyond the mass gap
(Fishbach & Holz, 2020).

Detecting the GW signal from high-redshift heavy stellar BBH mergers, where
one or both components are in the upper-mass gap or straddling it, would be
highly informative. Various mechanisms could lead to the formation and
coalescence of such binaries, and among them, evolution in isolated binaries,
dynamical encounters in star clusters or a chain of Nth generation
mergers.333There are several proposed channels for the origin of the
components of GW190521-like systems observed at low redshift. For instance:
mergers from the relics of the first stars (known as Population III stars)
(Liu & Bromm, 2020; Tanikawa et al., 2021), isolated binary evolution with no-
gap (Belczynski, 2020), stellar collisions in star clusters (see Costa et al.,
2022; Arca-Sedda et al., 2021; Di Carlo et al., 2020; Renzo et al., 2020, and
references therein), and hierarchical mergers (Gerosa & Fishbach, 2021). But
at redshifts as high as $\sim 10$ the contributions from the two dynamical
channels appear to be negligible (Mapelli et al., 2022a). Consequently,
observing these systems at high $z$ would allow us to better probe the physics
of the isolated binary channel, the potential existence of an upper-mass gap
and its imprint on the mass function of the earliest stellar BHs.

Estimates of the location and width of the upper-mass gap are at best
approximate. Current uncertainties on the reaction rate, on rotation, and on
the presence of rich hydrogen envelopes may shift the instability interval for
the explosion to higher masses, and narrow further the gap (Farmer et al.,
2020; Belczynski, 2020; Woosley & Heger, 2021; Marchant & Moriya, 2020; Vink
et al., 2021) or even fill it entirely (Costa et al., 2021). Testing the
existence of this upper-mass gap and inferring it properties from GW
observations depends critically upon the accuracy with which the masses of the
individual BHs are measured from the merger signal. Here and for this purpose
we carry out a parameter estimation on BBHs with component masses which touch
the edges of the upper-mass gap, recognizing that all of the above arguments
become compelling if the redshift of the observed systems is $z\gtrsim 10.$

Determining the lowest redshift one can claim the source to be beyond, and
inferring posteriors for the distribution of the component masses is of
paramount importance (Mancarella et al., 2023). A key challenge, then, is to
accurately infer both the masses and redshift of the binary. There is a well-
known degeneracy between measurements of the distance to and inclination of a
binary from GW observations (Usman et al., 2019). At high redshifts, this
degeneracy further impacts our ability to infer the masses of the binary. In
GW observations, masses and redshift are degenerate, and only the redshifted
masses, $m_{1,2}(1+z)$, can be inferred from the signal. Given a cosmological
model, the measured distance can be used to obtain the redshift and hence the
source masses. However, if the distance is poorly constrained, this leads to
significant uncertainties on the mass. For example, it is not unusual to have
an uncertainty of $\sim 50\%$ in the distance measurement of BBH signals (The
LIGO Scientific Collaboration et al., 2021b). At a $z=10$ this translates to a
redshift uncertainty of $\pm 4$ and consequently an uncertainty in the masses
of $40\%$ due to redshift effects alone. The ability to accurately infer
redshifts and masses is improved by a detector network, which can provide more
accurate localization and distance measurements (Fairhurst, 2011; Mills et
al., 2018; Singer et al., 2014), as well as the observation of Higher order
Multipoles in the GW signal which help break the distance-inclination
degeneracy (Mills & Fairhurst, 2021; Fairhurst et al., 2023).

The paper is organized as follows. In Section 2, we discuss the observability
of high-mass, high-redshift binaries with a focus on the HoMs. In Section 3,
we provide detailed parameter estimation results for a number of
astrophysically interesting simulated BBH merger signals and in Section 4 we
summarize our results. We include two appendices. Appendix A provides
additional figures showing detector sensitivity for binaries of varying mass
ratio and Appendix B gives parameter estimation accuracy for low signal-to-
noise ratio (SNR) systems.

##  2 The Importance of Higher order Multipoles

High-mass, high-$z$ BBH coalesences are intrinsically low-frequency GW
sources. This is illustrated in Fig.1, where we show the frequency evolution
of the GW strain amplitude for a BBH of (120-60)$\,\rm M_{\odot}$ placed at
redshift $z=14$, with an inclination $\iota$ of $60\degree$ between the
orbital angular momentum and theof line of sight. The gravitational waveform
for this signal only extends to $15$ Hz and is therefore outside the sensitive
frequency range of current GW observatories. The leading GW emission from the
source, emitted in the (2, 2) multipole at twice the orbital frequency,
extends to only $7$ Hz in the detector, making discovery challenging.444We
recall that for all BBH observed to date, the (2, 2) multipole, which is
emitted at twice the orbital frequency, has been the dominant multipole
detected in the GW signal. Additional multipoles of the GW signal have been
observed for a handful of events (Abbott et al., 2020a, b) but, as their
amplitudes are lower, they are generally not identified for the majority of
sources. Although the (3, 3) and (4, 4) multipoles are intrinsically lower
amplitude, they extend to higher frequencies ($\sim$1.5 and 2 times the
frequency of the (2, 2)) and can therefore contribute significantly to the
observed SNR. This improves the prospects of detecting such a system.
Furthermore, the identification of these higher-order (higher-frequency)
multipoles in the signal can significantly improve the ability to infer the
parameters of the system, as they enable us to break measurement degeneracies
that exist with observation of only a single multipole.

Figure 1: Waveform amplitude from a (120,60) $\,\rm M_{\odot}$ binary at
$z=14$, inclined at an angle of $60\degree$ to the line of sight. The figure
shows the frequency evolution of the waveform, and also the individual
contributions from the three leading multipoles (2, 2), (3, 3) and (4, 4)
using the IMRPhenomXPHM waveform (Pratten et al., 2020). The oscillations in
the full waveform are due to constructive/destructive interference between the
multipoles. The projected noise curves for three future detectors: Cosmic
Explorer, Einstein Telescope and A# are also shown. The signal is represented
in such a way that area between the wave amplitude and noise is indicative of
the SNR (Abbott et al., 2016).

There are several well-known degeneracies in the emitted gravitational-
waveform, leading to some parameters being very well measured while others
being not. For our purposes, we are most concerned with a degeneracy between
the observed distance to and inclination of a binary, as discussed in Usman et
al. (2019). When only the (2, 2) multipole is observed, the amplitude gives a
good measurement of $\cos\iota/d_{L}$ where $\iota$ is the binary inclination
and $d_{L}$ is the luminosity distance. However, in many cases, the binary
inclination is only restricted to the range $\iota\in[0\degree,60\degree]$,
leading to a factor of two uncertainty in distance due to this degeneracy
alone. When the binary is observed at a high redshift, the measurement of the
masses also becomes degenerate with distance and inclination, and a factor of
two uncertainty in distance can lead to a similar uncertainty on the masses.
The observation of a second GW multipole can serve to break this degeneracy
(Abbott et al., 2020a; Mills & Fairhurst, 2021) as the relative amplitude of
the different multipoles depends upon the orientation of the binary555The
ratio of the amplitude of the (3, 3) multipole to the (2, 2) scales as
$\sin\iota$ while the (4, 4) multipole scales as $\sin^{2}\iota$ relative to
the (2, 2)..

Figure 2: The total SNR and SNR in each multipole in ET for a binary of total
mass of $180\,\rm M_{\odot}$, at redshift $z=14$ and overhead the detector.
Left: Variation of SNR with inclination of the binary, for a BBH of
($120,60)\,\rm M_{\odot}$ and the fixed redshift. Right: Variation of SNR with
mass ratio, for a $180\,\rm M_{\odot}$ binary with inclination of $60\degree$
at the fixed ratio. The total SNR for the system is computed using the (2, 2),
(2, 1), (3, 3), (3, 2) and (4, 4) multipoles and individual SNRs for the (2,
2), (3, 3) and (4, 4) multipoles are shown. At $\iota=0$ in the left panel the
full waveform SNR is higher than the (2,2) mode due to the contribution of the
(2, 1) and (3, 2) modes, which are not shown. The total SNR depends upon the
phase of the signal at merger, due to non-zero overlap between multipoles, and
is therefore shown as a shaded blue area. In both panels, horizontal lines
indicate SNR = 8, an approximate threshold for detection, and SNR = 3, an
approximate threshold for observation of a second multipole, given that
another multipole has been observed (Mills & Fairhurst, 2021).

In Fig. 2, we show the variation of the SNR with binary mass ratio
$q=m_{1}/m_{2}$ (assuming an inclination $\iota=60^{\circ}$) and inclination
(assuming $q=2$) in each of the multipoles for a binary of total mass of
$180\,\rm M_{\odot}$ at $z=14$ observed by ET. The SNR of the (2, 2) multipole
is greatest for face-on signals ($\iota=0^{\circ}$) with equal mass components
($q=1$). For a face-on signal, the (2, 2) multipole is circularly polarized
and, as the inclination increases, the amplitude of both polarizations
decreases to a minimum for edge-on systems, $\iota=90^{\circ}$, whose emission
is linearly polarized. For the other multipoles considered, the SNR vanishes
at face-on and peaks at $\sim 50\degree$ for the (3, 3) multipole and $\sim
60\degree$ for the (4, 4) multipole. The binary would be observable in ET at
any orientation. For inclinations $\iota\gtrsim 10\degree$ or $30\degree$ the
(3, 3) and (4, 4) multipoles would be identifiable, respectively. Since this
waveform lasts only a few cycles in the detector band, the contributions from
the different multipoles are not orthogonal. Consequently, the total SNR
varies with the merger phase of the binary.

The SNR of each different multipole, and the full signal, also varies with
mass ratio. The (2, 2) multipole is largest for equal mass systems and
decreases by a factor of two by mass ratio $q=5$, while the (3, 3) vanishes
for equal mass and peaks around $q=3$. For this signal, the SNR in the (4, 4)
multipole does not vary significantly with mass ratio. The (2, 2) and (4, 4)
multipoles would be identifiable at any mass ratio, and the (3, 3) for
binaries with mass ratio above $\sim$1.2. Identification of more than one
multipole enables an improved measurement of mass ratio, as well as binary
orientation.

Figure 3: The same as Fig. 2 but for CE.

In Fig. 3, we show the same dependence of SNR with inclination and mass ratio
for the CE detector. Since CE has sensitivity to the signal above 5 Hz, rather
than 3 Hz for ET, the overall SNR is lower and the signal would be marginally
observable. Furthermore, a broad range both in inclination and mass ratio,
where the (3,3) multiple gives the dominant contribution to the SNR, becomes
accessible. This provides a clear example of a signal where the HoMs enable
detection as well improved parameter recovery.

Figure 4: Maximum observable redshift of binaries with mass ratio $q=2$ with
next-generation observatories. Left: ET, Right: CE. Shaded regions show
redshifts at which 10%, 50% and 90% of sources are observable. The maximum
observable redshift of the (2, 2) multipole at SNR = 8 is shown in orange, the
observability of the full waveform at SNR = 8 is shown in blue, and the
observability of the second most significant multipole above SNR = 3 is shown
in black. The red asterisks represent the masses and redshifts of the $q=2$
systems we study in detail in Section 3.

Given the above discussion, we are interested in identifying the regions of
the mass space where HoMs can contribute to either the observability or
parameter measurement accuracy of high-mass, high-redshift binaries. In Fig. 4
we show the sensitivity of the proposed ET and CE observatories to BBH mergers
with mass ratio of 2 as a function of redshift. We show the maximum redshift
at which a binary can be observed, at an SNR of 8, and also the redshifts at
which 10%, 50% and 90% of binaries, averaged over sky location and
orientation, will be observed. The detector sensitivity is shown for both the
(2, 2) multipole, in orange, and the full waveform, in blue. At low masses,
the (2, 2) multipole dominates the observable signal and therefore the
distance to which the full waveform can be observed is essentially equal to
that of the (2, 2) multipole. However, at high masses, the (3, 3) and (4, 4)
multipoles contribute more significantly and incorporating them increases the
sensitivity of the detectors to these systems. When a system has been
observed, the identification of a second multipole, at SNR above 3, can
greatly improve parameter recovery by breaking degeneracies between distance
and inclination and improving mass ratio measurement. The range of masses and
redshifts for which the binary would be observed with SNR above 8, and with
SNR above 3 in at least two multipoles, is shown in black in Fig. 4. For
example, in ET a 4,000$\,\rm M_{\odot}$ system is visible at $z\approx 1$ with
the (2, 2) multipole but up to $z\approx 2$ with the full waveform.

Remarkably, for the majority of binaries with $M\gtrsim 100\,\rm M_{\odot}$
observed by ET, and $M\gtrsim 30\,\rm M_{\odot}$ observed by CE, a second
multipole will be observable. At lower masses, it is the (2, 2) and (3, 3)
multipoles which contain most power, while at high masses it is the (3, 3) and
(4, 4) multipoles that are observed, with the (2, 2) multipole power occurring
at frequencies below the instrumental sensitivity. The picture is similar at
different mass ratios, and figures showing the sensitivity to binaries with
$q=1,2,4$ and $10$ are provided in Appendix A, for ET in Fig. 11 and CE in
Fig. 12. The most significant difference occurs for equal mass binaries, where
the (3, 3) multipole vanishes and we therefore require both (2, 2) and (4, 4)
multipoles to be observable. This limits the range for which two multipoles
can be seen and increases the minimum mass at which we expect to observe two
multipoles to $\sim 200\,\rm M_{\odot}$ in ET and $\sim 50\,\rm M_{\odot}$ for
CE. Nonetheless, for the majority of high-mass, high-redshift binaries, we
expect to observe multipole multipoles, and therefore obtain good estimates of
both the masses and redshift of the system. In the next section, we
investigate those expectations in detail through parameter recovery of a
series of systems.

##  3 Parameter recovery for high-mass, high-redshift binaries

Observation in a single GW observatory leads to large uncertainties in the sky
location of the binary (Singer et al., 2014), and this is again degenerate
with the inferred distance and redshift. A network of detectors with
comparable sensitivity can provide accurate localization (Mills et al., 2018)
and therefore improved redshift and mass accuracy. Binaries with black hole
spins misaligned with the orbital angular momentum will precess. In principle,
the observation of precession can further improve parameter estimates.
However, given that so few cycles of the waveform are visible in the
detectors, the prospects for observing precession are slim (Green et al.,
2020). Therefore, in what follows we neglect precession effects.

To illustrate the expected performance of a next-generation GW network in
observing and measuring these binaries, we perform a number of simulations and
obtain parameter estimates with the LALInference (Veitch et al., 2015) package
and a uniform in comoving volume distance prior. We simulate four different
binary mass combinations, denoting (in the source frame) with $m_{1}$
($m_{2}$) the primary (secondary) mass and with $M$ the total mass. We
consider $(120,60)\,\rm M_{\odot}$ and $(90,90)\,\rm M_{\odot}$ binaries,
chosen so that component BHs lie in, or close to the upper mass-gap, and
$(240,120)\,\rm M_{\odot}$ and $(480,120)\,\rm M_{\odot}$ binaries chosen to
probe observability of high-redshift intermediate mass black hole (IMBH) in
binaries. In all cases, we simulate quasi-circular non-spinning BBHs, but
allow for non-zero, aligned spins when performing parameter estimation666The
restriction to non-spinning BBH is solely to simplify presentation — all
results presented here could be easily extended to aligned-spin BHs.. This is
important as the degeneracy between the binary mass ratio and BH spins (Hannam
et al., 2013; Fairhurst et al., 2023) greatly impacts the accuracy with which
mass ratio can be measured.

The simulated signals are added to data from a three-detector network of
observatories with sensitivity matching ET (Punturo et al., 2014) and CE
(Evans et al., 2021). Specifically, we use a single, triangular ET detector
located in Europe and two 40 km CE observatories, one in the US and one in
India. The simulations are performed at the optimal sky location for the
network. Given the greater low-frequency sensitivity of ET, this leads to the
binaries being essentially overhead ET. The signals are generated at varying
inclination angle, to enable us to investigate the importance of HoMs. We
choose the redshift of the sources to ensure a fixed SNR for all signals. In
the main text, we use an SNR of 30, while in Appendix B we investigate quieter
signals with an SNR of 15.

###  3.1 Observing mass-gap objects

A mass gap in the BH mass distribution is expected due to the presence of the
PISN. Farmer et al. (2020) investigated the location of this pair-instability
region as a function of the temperature-dependent uncertainty in the
${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction rate. Determining the value of
the ${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction rate is extremely important
for tracing the evolution of massive stars. Thus, restricting this rate
through GW observations would be of considerable astrophysical interest.
According to Farmer et al. (2020), the width of the mass-gap remains roughly
constant as a function of the (unknown) reaction rate, but the mass-range
where no black hole can form varies. At the lowest rate relative to the
median, the mass-gap extends from $\sim 90\,\rm M_{\odot}$ to $\sim 175\,\rm
M_{\odot}$. At the highest rate, the location of the mass-gap is between $\sim
60\,\rm M_{\odot}$ and $\sim 120\,\rm M_{\odot}$. Interestingly, there exists
a region of BH masses between $90\,\rm M_{\odot}$ and $120\,\rm M_{\odot}$
where we should not expect any black hole to form, for any rate. We refer to
this region as the forbidden strip. Consequently, we choose to investigate
systems which host at least one member with mass touching this narrow strip.
Then, if their masses were to be determined with sufficient accuracy, their
detection could constrain the ${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction
rate to be at the extreme of the allowed range (see Fig. 5 of Farmer et al.,
2020). In particular, we focus on $(120,60)\,\rm M_{\odot}$ and $(90,90)\,\rm
M_{\odot}$ binaries which have components at the lower and upper range of the
forbidden strip.

As seen in Fig. 4, BBHs with masses $(120,60)\,\rm M_{\odot}$ will be
detectable at a maximum redshift of $z\sim 25$, for an optimally located and
oriented system, with 50% of mergers at $z\sim 17$ and the vast majority of
events at $z\sim 10$ being detectable. The sensitivity to $(90,90)\,\rm
M_{\odot}$ systems is comparable. If $(90,90)\,\rm M_{\odot}$ systems were to
be observed, this will allow us to constrain the strength of the uncertain
${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction rate. In particular, black holes
with such masses would imply the rate to be at the lower end of the explored
range. A binary with mass $(120,60)\,\rm M_{\odot}$ would be challenging to
form through stellar evolution. Specifically, allowing for variation in the
reaction rate, the mass of the primary would require a very high reaction rate
for ${}^{12}\rm C(\alpha,\gamma)^{16}O$, while the mass of the secondary would
be compatible with a value below the median. Therefore this would be a system
where only one of the two black holes could originate from stellar evolution
and the other would require a different formation channel.

Figure 5:  Inferred posteriors of component masses for binaries of total mass
of $180\,\rm M_{\odot}$ which straddle, or sit within the upper mass gap.
Results are shown for inclinations of $\iota=0\degree,30\degree,60\degree$ and
$90\degree$, with the redshift of the system varied, as indicated in the
labels, so that the signals are all simulated with an SNR of 30 in the ET–CE
network, at the most sensitive sky location for the network. The grey region
denotes the pair-instability mass-gap and the darker grey denotes the
forbidden strip where no black hole is expected to form for any value of the
${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction rate (Farmer et al., 2020).
Simulated values are denoted by a black cross and contours show the 90%
credible region. Left: Binary with masses just below the mass gap:
$m_{1}=m_{2}=90\,\rm M_{\odot}$. Right: Binary with masses that straddle the
mass gap: $m_{1}=120\,\rm M_{\odot}$ and $m_{2}=60\,\rm M_{\odot}$. Greek
letter $\rho$ indicates the SNR, throughout the paper.

Figure 6: Inferred posteriors for the redshift for binaries of total mass of
$180\,\rm M_{\odot}$ which straddle, or sit within the upper mass gap. Signals
are simulated with varying inclination (in degrees) and redshift ($z$), as
indicated in the labels, with a fixed SNR of 30 in the ET–CE network, at the
most sensitive sky location for the network. Symbols $\rho_{\ell m}$ indicate
the SNR in the $(\ell,m)$ multipole HoMs. Black vertical lines indicate the
90% credible intervals, violins show the 99.5% range and yellow vertical lines
show simulated values. Left: Binary with masses just below the mass gap:
$m_{1}=m_{2}=90\,\rm M_{\odot}$. Right: Binary with masses that straddle the
mass gap: $m_{1}=120\,\rm M_{\odot}$ and $m_{2}=60\,\rm M_{\odot}$.

In Fig. 5 and 6, we show the recovered parameter accuracies for both mass and
redshift for $(90,90)\,\rm M_{\odot}$ and $(120,60)\,\rm M_{\odot}$ binaries
observed with SNR = 30. The first thing to note is that these high-mass, high-
redshift systems could be identified with good accuracy by the next-generation
GW network, as would be expected due to the relatively large SNRs. For all
events, there is, at most, a factor of two uncertainty in the mass of the
systems and a $50\%$ uncertainty in the redshift, with both numbers quoted at
the 90% confidence interval. However, we also notice a substantial variation
in the accuracy of parameter measurement between the systems. The parameters
of systems close to face-on ($\iota=0\degree$ or $30\degree$) are recovered
with significantly larger uncertainties than those which are more highly
inclined ($\iota=60\degree$ or $90\degree$). When the binary is close to face-
on, the uncertainty in the component mass posterior is much greater and, for
$\iota=0\degree$, the 90% mass region for the $(120,60)\,\rm M_{\odot}$ binary
includes equal masses. As the inclination of the binary is increased,
parameter accuracy improves significantly: already at $\iota=30\degree$ the
posterior for the $(120,60)\,\rm M_{\odot}$ binary is inconsistent with equal
masses, although large uncertainty in the mass ratio remains. For binaries
inclined at $60\degree$ or $90\degree$ the parameter accuracy is excellent. In
both cases, the mass ratio is very well constrained and uncertainty in total
mass and redshift is $\pm(10-20)\%$.

Figure 7: Inferred posteriors for redshift, mass ratio ($1/q=m_{2}/m_{1}$) and
inclination for $(120,60)\,\rm M_{\odot}$ binary observed at SNR=30. Top:
Parameter recovery for a system at $z=14$ and $\iota=60\degree$. Bottom:
Parameter recovery for a system at $z=17$ and $\iota=0\degree$.

We can explain why the results of our parameter accuracy display the features
they do, using our understanding of the importance of HoMs discussed in
Section 2. Let us begin by considering the $(120,60)\,\rm M_{\odot}$ system at
$z=14$, inclined at $\iota=60\degree$, that we have plotted in Figs. 1, 2 and
3. The signal is observed with high SNR in ET and the (3, 3) multipole is also
clearly observed. In CE, the total SNR of the system is sufficient for it to
be observed, with the (3, 3) multipole providing the dominant contribution.
Since the event is observed in three detectors, it is relatively well
localized,777Since these are very low-frequency systems, the localization is
poorer than for events in GWTC-3 (The LIGO Scientific Collaboration et al.,
2021b) as the localization depends upon the frequency bandwidth of the signal
(Fairhurst, 2009). with a 90% localization area of $300\deg^{2}$. The
observation of HoMs in the waveform enables the accurate inference of both the
binary inclination and mass ratio. Since the sky location, mass ratio and
inclination are well measured, this enables accurate inference of the distance
to the binary and consequently the redshift of the source. In the top panel of
Fig. 7 we show the recovered values for the redshift, binary inclination and
mass ratio. All three are recovered with an accuracy better than 10%.

Next, we consider a comparable system observed face-on ($\iota=0\degree$) at
the same SNR, and at a redshift of $17$. In that case, the power in both the
(3, 3) and (4, 4) multipoles vanishes. Consequently, the binary is no longer
observable in CE as it has an SNR of 1.888The waveform shown in Fig. 3
corresponds to a binary at $z=14$ and we are now considering the same mass
binary at $z=17$. The larger redshift reduces the SNR primarily through
redshifting the signal which lowers the frequency by 20%. This leads to an SNR
which is lower than that shown in the figure. At this SNR, the CE
observatories are unable to provide localization of the source, with a 90%
localization area of $10,000\deg^{2}$. Furthermore, the vanishing HoMs means
that only the (2, 2) multipole is observed in ET. Therefore, it is not
possible to break the degeneracy between binary orientation and distance, nor
to place a tight constraint upon the mass ratio. The bottom panel of Fig. 7
shows a the recovered redshift, mass ratio and inclination for this system.
The mass ratio is not accurately recovered and, indeed, the binary is inferred
most likely be (close-to) equal mass, although the distribution does extend to
$1/q=0.5$. In addition, the binary orientation is not accurately recovered,
with a broad distribution of $\iota\lesssim 25\degree$ — more inclined systems
are excluded as they would have observable power in the (4, 4) multipole. The
mass ratio–inclination distribution does show a secondary peak close to the
simulated value ($1/q\approx 0.5$ and $\iota<10\degree$), however, the
preference is for an equal mass system. Despite both mass ratio and
inclination being offset from the true values, the inferred redshift matches
well with the simulated value. However, due to the uncertainties in other
parameters, the redshift uncertainty is now close to $25\%$. The comparison of
parameter accuracy for these two systems highlights the importance of both a
network of detectors and also observability of the HoMs in accurate inference
of binary properties.

It is worth noting that our intuition from current GW observations that the
majority of sources are close to face-on (or face-off) no longer holds in the
next-generation network (Vitale, 2016). In the nearby Universe, where sources
are approximately uniform in volume, a signal observed with a given SNR is
most likely to originate from a distant binary which is close to face-on (or
face off) as the number of sources increases as $d_{L}^{2}$. For a high-
redshift source, whose redshift is past the peak of the redshift distribution
— likely around a redshift of $z\approx 2$ at which star formation peaks —
this is no longer the case. Now, the most likely origin is from a binary which
is at lower redshift, where the intrinsic rate is higher, and is either poorly
oriented or from a region of the sky where the detectors have lower
sensitivity. Thus, the results from sources inclined at $60\degree$ and
$90\degree$ are more typical of the observed population.

Let us return to the implications for probing the location of the Pair
Instability Supernova (PISN) mass-gap. For both mass pairs, binaries inclined
at $60\degree$ or $90\degree$ are those which provide the best mass
measurements. For the $(90,90)\,\rm M_{\odot}$ system, we have $m_{1}$ and
$m_{2}$ measured in the interval $\in[70,100]\,\rm M_{\odot}$. So, this system
is consistent with both components lying below the mass-gap provided the
${}^{12}\rm C(\alpha,\gamma)^{16}O$ rate is low. For the $(120,60)\,\rm
M_{\odot}$ system, we have $m_{1}\in[100,140]\,\rm M_{\odot}$ and
$m_{2}\in[40,80]\,\rm M_{\odot}$. The masses are consistent with one above and
one below the gap, provided the reaction rate is high. If both signals were
observed, this would be inconsistent with our current understanding of the
PISN mass gap.

To investigate the observability at even higher redshifts, we have simulated a
second set of signals, with the same masses and inclinations but with a lower
SNR fixed at 15. For these systems, the redshifts range from $z\gtrsim 20$ for
face-on systems to $z\approx 15$ for edge on systems. Broadly, the results are
consistent with those in Fig. 5 and 6, but with larger uncertainties due to
the lower signal amplitude. In particular, for all but the face-on systems, we
are able to clearly identify that the $(120,60)\,\rm M_{\odot}$ binary is of
unequal mass, due to the observed power in HoMs. For the inclined systems, the
uncertainty in total mass and redshift is around a factor of two (from 150 to
$300\,\rm M_{\odot}$ and $z=12$ to 25). Thus, while it is possible to identify
these systems as unambiguously high-mass and high-redshift sources, the
uncertainties in masses and redshifts make it difficult to perform precision
astrophysics. For the $(90,90)\,\rm M_{\odot}$ system, it is only at
$\iota=60\degree$ or $90\degree$ that the parameters are well recovered. These
results are shown in Fig. 13 and 14 in Appendix B.

Figure 8: Inferred posteriors for redshift and redshifted mass for a simulated
$(120,60)\,\rm M_{\odot}$ BBH merger with $\iota=0\degree$ observed at $z=21$.

For the face-on systems, we see an interesting feature whereby the binary can
be mistaken for a different system with very different properties. Fig. 8
shows the inferred redshift and redshifted mass –$M(1+z)$ – distributions for
a $(120,60)\,\rm M_{\odot}$ system at $z=21$. The primary peak is at $z=21$
and $M(1+z)\approx 4,000\,\rm M_{\odot}$ corresponding to the simulated value.
However, there is a secondary peak around $z=5$, with a redshifted mass around
$6,000\,\rm M_{\odot}$ corresponding to a binary with mass of $1,000\,\rm
M_{\odot}$. For this system, it is the (3, 3) multipole which is consistent
with the simulated signal. This provides another example of the challenges
which arise when identifying high-mass, high-redshift binaries. The signal
would be observed only in ET, and have only one observable multipole. Not only
does this lead to poor parameter recovery, but also the inability to
distinguish between a $180\,\rm M_{\odot}$ binary at $z=21$ and a $1,000\,\rm
M_{\odot}$ binary at $z=5$. Given the GW data alone, it would not be possible
to distinguish between the two scenarios. The relative significance of the two
will depend upon astrophysical knowledge of the mass and redshift
distributions of BBHs. Here, we have used priors which are uniform in comoving
volume and component masses. Other choices might lead to different conclusions
about the mass and redshift of the binary.

In summary, for these two representative sources, and for high SNRs we would
confidently identify the systems as high-mass BBHs at high redshift. These
could be potential seeds for the growth of supermassive black holes. The first
system of $(90,90)\,\rm M_{\odot}$ would be marginally consistent with being a
binary formed at the lower edge of the mass gap, in correspondence of the
lowest value of the ${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction rate.
Inconsistent otherwise. The second system of $(120,60)\,\rm M_{\odot}$ would
be consistent with one BH (the lightest) originating from the core-collapse of
a massive star (provided the ${}^{12}\rm C(\alpha,\gamma)^{16}O$ rate is low)
and the second, in the midst of the pair-instability gap, would have a
different origin. For a large fraction of the computed rates, the discovery of
the latter system would be inconsistent with the explosion scenario implied by
the pair instability that would predict no BHs, and therefore a different
channel has to be called for, for both components.

###  3.2 Observing intermediate mass black holes in binaries

Next, let us consider higher mass systems, containing an IMBH. For
concreteness, we consider binaries with masses $(240,120)\,\rm M_{\odot}$ and
$(480,120)\,\rm M_{\odot}$ which would be observable up to $z\sim 10$. While
the $120\,\rm M_{\odot}$ BH would be at the border of the pair-instability
gap, both BHs are well above the gap. Their observation would either require a
sufficiently high upper-mass end of a Kroupa-like IMF (extending to at least
$300\,\rm M_{\odot}$) or a top-heavy IMF (Klessen & Glover, 2023).
Alternatively the primary BH (particularly the most massive one of $480\,\rm
M_{\odot}$) could have had time to increase its original mass due to accretion
of gas from its surroundings. Thus for these systems, we are interested in
determining whether the mass of the primary and the redshift can be accurately
inferred in order to identify these early IMBHs.

Assuming the mass gap to be that predicted by the median of the ${}^{12}\rm
C(\alpha,\gamma)^{16}O$ reaction rate, between $\sim 50\,\rm M_{\odot}$ and
$130\,\rm M_{\odot}$, and a standard Kroupa-like IMF, in the interval between
$0.1\,\rm M_{\odot}$ and $150\,\rm M_{\odot}$, all our systems are expected to
have dynamical origin in dense star clusters (see, eg. Mapelli et al.
2021).999We warn that black holes with these masses could be of primordial
origin or outcome of post-formation accretion, as mentioned in the
Introduction. Among these, Nuclear Star Clusters (NSCs) could be the sites for
3rd- or 4th-generation BBHs with observed individual masses up to $\sim
600\,\rm M_{\odot}$, merging on timescales smaller than 500 Myr compatible
with the redshift of observation. Formation in Globular Clusters (GCs) would
be marginally compatible with our lightest systems. According to Mapelli et
al. (2022b), above $z\sim 10$ there are two available BBH formation channels,
the isolated and the dynamical formation in young star clusters. However, the
maximum masses for these channels are $\sim 50\,\rm M_{\odot}$ and $\sim
100\,\rm M_{\odot}$ respectively. This implies that detecting BBHs with
individual masses larger than $\sim 100\,\rm M_{\odot}$ at $z>10$ could point
to a top-heavy IMF, as predicted for the first stellar generation.

Figure 9:  Inferred posteriors for the component masses for heavy binaries.
Results are shown for inclinations of $\iota=0\degree,30\degree,60\degree$ and
$90\degree$, with the redshift of the system varied, as indicated in the
labels, so that the signals are all simulated with an SNR of 30 in the ET-CE
network, at the most sensitive sky location for the network. The grey region
denotes the pair-instability mass-gap and the darker grey denotes the
forbidden strip where no black hole is expected to form for any value of the
${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction rate (Farmer et al., 2020).
Simulated values are denoted by a black cross and contours show the 90%
credible region. Left: $m_{1}=240\,\rm M_{\odot}$ and $m_{2}=120\,\rm
M_{\odot}$. Right: $m_{1}=480\,\rm M_{\odot}$ and $m_{2}=120\,\rm M_{\odot}$.

Figure 10: Inferred posteriors for the redshift for heavy binaries. Signals
are simulated with varying inclination and redshift, as indicated in the
labels, with a fixed SNR of 30 in the ET–CE network, at the most sensitive sky
location for the network. Black vertical lines indicate the 90% credible
intervals, violins show the 99.5% range and yellow vertical lines show
simulated values. Symbols $\rho_{\ell m}$ indicate the SNR in the $(\ell,m)$
multipole HoMs. Left: $m_{1}=240\,\rm M_{\odot}$ and $m_{2}=120\,\rm
M_{\odot}$. Right: $m_{1}=480\,\rm M_{\odot}$ and $m_{2}=120\,\rm M_{\odot}$.

Figures 9 and 10 show the accuracy with which we can measure the masses and
redshifts of the events. The broad features are similar to what we have
already observed for the lower mass systems, namely that the parameter
recovery is significantly worse for face-on systems, due to the vanishing
HoMs. Even though the mass ratio of the systems is 2 or 4, both are inferred
to be consistent with equal mass (or nearly equal mass) binaries when viewed
face-on. Furthermore, the uncertainty in redshift and total mass is about a
factor of two. For the inclined systems, the recovery of masses and redshifts
improves significantly, particularly for $\iota\geq 60\degree$. In that case,
component masses and redshifts are recovered with a $\sim 20\%$ accuracy. In
particular, the mass of the $120\,\rm M_{\odot}$ BH will be constrained to be
between $90\,\rm M_{\odot}$ and $\sim 150\,\rm M_{\odot}$ for all except for
the face-on system. This is consistent with a black hole above the mass-gap
and, due to uncertainties in the mass measurement, will not significantly
restrict the ${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction-rate. In Appendix B,
we also show results for events simulated at higher redshifts and at a lower
SNR of 15. The results are comparable to those discussed above, with the
masses and redshifts for inclined systems better measured, and masses
constrained to be unequal. For face-on systems, the parameter recovery is
significantly worse and we again see multiple peaks in the mass-redshift
distributions corresponding to different multipoles matching with the signal.

Remarkably, the next-generation GW observatories have the capability to detect
and accurately identify mergers involving $240\,\rm M_{\odot}$ BH at a
redshift of 10, and confidently infer a minimum redshift of 7, and mergers
involving a $480\,\rm M_{\odot}$ BH at a redshift of 6, and infer a redshift
of at least 4. These systems will be interesting to observe because we do not
know if BHs of those masses exist, and we can hope to shed light on their
formation routes, either by accretion from lower-mass BHs or by direct
collapse of very massive stars.

##  4 Discussion

The next-generation of GW detectors provide a unique way to probe the
existence of heavy stellar black holes in the high-redshift Universe. Future
GW observations of BHs with masses above $\sim 50\,\rm M_{\odot}$ at redshift
$z\sim 10-15$ will enable us to probe the properties of the first stars
forming in the Universe and their initial mass function. If BHs in the mass
range explored here exist, they can contribute, as seeds 101010the so called
“light seeds” explored in the literature (Valiante et al., 2017; Inayoshi et
al., 2020; Volonteri et al., 2021)., to the rapid growth of the population of
quasars observed close to the recombination epoch, at $z\approx 7.5$, and
housing accreting BHs of ${\cal O}(10^{8}-10^{9})\,\rm M_{\odot}$ (e.g.
Volonteri & Rees, 2005; Madau et al., 2014; Pezzulli et al., 2016; Volonteri
et al., 2015; Ricarte & Natarajan, 2018; Valiante et al., 2021; Trinca et al.,
2022). Whether and how the bridge between stellar and supermassive black holes
was established when the first galaxies were forming is currently unknown
(Volonteri et al., 2021). The revolutionary data coming from JWST, with the
recent discovery of more than 40 new faint accreting supermassive BHs of
${\cal O}(10^{5}-10^{7}\,\rm M_{\odot})$ at $4<z\leq 10.6$ (see Maiolino et
al., 2023, and references therein), is an outstanding confirmation of the rich
BH landscape at cosmic dawn predicted by theoretical models (Schneider et al.,
2023). In the future, with the Laser Interferometer Space Antenna (LISA) in
operation (Amaro-Seoane et al., 2017), we will detect low-frequency GWs from
merging massive BHs of ${\cal O}(10^{4}-10^{6}\,\rm M_{\odot})$ out to $z\sim
10-15$. By combining and confronting statistically all observations of both
merging and accreting BHs, we will be able to shed light into the origin and
evolution of the BH populations, from the stellar to the supermassive through
the intermediate-mass ones, across the cosmic epochs (Fragione & Loeb, 2023;
Valiante et al., 2021).

In this paper, we focused on the observability of high-redshift stellar BBHs
with high masses, and, equally importantly, on the accuracy with which their
masses and redshifts can be inferred. We have shown that both the observation
of systems and the accurate measurement of their parameters depend critically
on the inclusion of HoMs in the GW waveform. At the highest masses and
redshifts, HoMs, which extend the signal to higher frequencies than the (2, 2)
multipole, can significantly increase the sensitive range of the detectors.
Across a broad range of masses and redshifts, we expect to see multiple GW
multipoles in signals observed by CE and ET. Observation of more than one
multipole, typically a HoM in addition to the (2, 2) multipole, enables the
breaking of the degeneracy between binary inclination and distance, as well as
a more accurate determination of the mass ratio. Additionally, a network of
observatories is required for source’s localizations. When a signal is seen in
only a single detector, the sky location is poorly measured and, since the
detector response varies significantly over the sky, this leads to large
uncertainties in the distance. For very high-redshift sources, accurate
distance/redshift measurement is vital for the measurement of the BH
individual masses, as the observed signal depends upon $M(1+z)$. By performing
full parameter estimation on a set of representative systems, we demonstrated
that it will be possible to measure masses and redshifts with an accuracy of
$10-20\%$, for signals at redshifts up to at least 15. Those systems which can
be observed and accurately measured are typically seen in both CE and ET
detectors, so they are well localized and also tend to be viewed away from
face-on (or face-off) so that more than one GW multipole is observed.

We examined systems with masses $(120,60)\,\rm M_{\odot}$ and $(90,90)\,\rm
M_{\odot}$, which lie in, or around the pair-instability mass-gap. For the
best-measured of the examples we investigated, at a redshift of 10, we could
measure redshift and component masses with 10% uncertainly. This would,
unambiguously, place these sources in the high-redshift Universe and serve to
constrain the, currently unknown, location of the pair-instability mass-gap.
We also investigated mergers of $(240,120)\,\rm M_{\odot}$ and $(480,120)\,\rm
M_{\odot}$ binaries, which enable us to probe the observability of early IMBH.
It will be possible to observe these IMBH at redshifts up to $z=10$ and
constrain the redshift to be at least $z=7$.

The results in this paper complement those in Ng et al. (2022, 2023) which
investigate lower-mass BH mergers in the next-generation GW network, and
Mancarella et al. (2023) who introduced the concept of an “inference horizon”
which maps the redshift a source can confidently be placed beyond. In all
cases, it is shown that next-generation GW network provides a unique
capability to probe high-redshift black hole formation.

The most critical feature of detector sensitivity for observing these systems
is the low-frequency sensitivity of the detectors. In our study, we used a
low-frequency limit of $3$ Hz for ET and $5$ Hz for CE. Even the relatively
small change from $3$ Hz to $5$ Hz can have a profound impact on sensitivity
to high-mass, high-redshift sources as shown in Fig. 2 and 3. Achieving the
instrumental design sensitivity at low frequencies has been challenging in the
current LIGO and Virgo observatories. As the detailed technical designs of the
next-generation observatories are finalised, the desire to probe the remnants
of high-mass stars in the early Universe should be considered as a motivation
to optimize sensitivity at low frequencies.

## Acknowledgements

We thank Riccardo Buscicchio for his careful reading and valuable comments. SF
acknowledges the support of STFC grant ST/V005618/1 and a Leverhulme Trust
International Fellowship. RS acknowledges support from the Amaldi Research
Centre funded by the MIUR programme “Dipartimento di Eccellenza”
(CUP:B81I18001170001). MC, RS, AT, RV acknowledge the INFN TEONGRAV specific
initiative. MC acknowledges support by the 2017-NAZ- 0418/PER grant, and by
the Italian Ministry for Universities and Research (MUR) program “Dipartimenti
di Eccellenza 2023-2027”, within the framework of the activities of the
“Centro Bicocca di Cosmologia Quantitativa (BiCoQ)”. AS acknowledges the
financial support provided under the European Union’s H2020 ERC Consolidator
Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement:
818691). MC and RS thank the Institut d’Astrophysique de Paris for kind
hospitality.

##  Appendix A Detector Sensitivity

Figure 11: Maximum observable redshift of BBH mergers with ET for mass ratios
$q=1,2,4$ and $10$. Shaded regions show redshifts at which 10%, 50% and 90% of
sources are observable. The maximum observable redshift of the (2, 2)
multipole at SNR = 8 is shown in orange, the observability of the full
waveform at SNR = 8 is shown in blue, and the observability of the second most
significant mode above SNR = 3 is shown in black. The red asterisks represent
the masses and redshifts of the ($q=2$) systems we study in detail in Section
3.

Figure 12: Same as 11 but for CE.

In this Appendix, we show the sensitivity of ET and CE to binary mergers as a
function of mass and redshift. Fig. 4 shows the sensitivity to mergers with
mass ratio of two. In Fig. 11 we show the ET sensitivity for binaries with
mass ratios 1, 2, 4 and 10. In Fig. 12, we show the same for CE. The maximum
reach of the detectors is for equal mass binaries. However, at equal mass the
(3, 3) multipole vanishes so there is a larger range for which only one
multipole is visible. As we increase the mass ratio, the maximum sensitive
redshift decreases, as the amplitude of the emitted GW also decreases.
However, the relative significance of the HoMs increases so that an increasing
fraction of sources will be observed with at least two multipoles.

##  Appendix B Low SNR signals

In this Appendix we present parameter estimation results for systems that
would have a network SNR of 15 in the ET-CE network described in Section 3, to
complement the results for SNR 30 presented in that section of the paper.
Since these events are at a lower SNR, they are also at a higher redshift,
with the $(120,60)\,\rm M_{\odot}$ and $(90,90)\,\rm M_{\odot}$ binaries at
redshifts between 23 and 15 (depending upon inclination). At these redshifts,
the signal is shifted to such low frequencies that it lies essentially outside
of the sensitive band of CE — the SNR of these events in CE is less than 2 in
all cases. Consequently, all sources are poorly localized on the sky, with a
typical 90% localization of thousands of square degrees.

Figure 13:  Inferred posteriors of component masses for binaries of total mass
of $180\,\rm M_{\odot}$ which straddle, or sit within the upper mass gap.
Results are shown for inclinations of $\iota=0\degree,30\degree,60\degree$ and
$90\degree$, with the redshift of the system varied, as indicated in the
labels, so that the signals are all simulated with an SNR of 15 in the ET-CE
network, at the most sensitive sky location for the network. The grey region
denotes the mass gap, whose location varies with the ${}^{12}\rm
C(\alpha,\gamma)^{16}O$ reaction rate, and the darker grey denotes the
forbidden strip (Farmer et al., 2020). Simulated values are denoted by a black
cross and contours show the 90% credible region. Left: Binary with masses just
below the mass gap: $m_{1}=m_{2}=90\,\rm M_{\odot}$. Right: Binary with masses
that straddle the mass gap: $m_{1}=120\,\rm M_{\odot}$ and $m_{2}=60\,\rm
M_{\odot}$.

Figure 14: Inferred posteriors for the redshift for binaries of total mass of
$180\,\rm M_{\odot}$ which straddle, or sit within the upper mass gap. Signals
are simulated with varying inclination and redshift, as indicated in the
labels, with a fixed SNR of 15 in the ET–CE network, at the most sensitive sky
location for the network. Black vertical lines indicate the 90% credible
intervals, violins show the 99.5% range and yellow vertical lines show
simulated values. Symbols $\rho_{\ell m}$ indicate the SNR in the $(\ell,m)$
multipole HoMs. Left: Binary with masses just below the mass gap:
$m_{1}=m_{2}=90\,\rm M_{\odot}$. Right: Binary with masses that straddle the
mass gap: $m_{1}=120\,\rm M_{\odot}$ and $m_{2}=60\,\rm M_{\odot}$.

In Fig. 13 and 14, we show the accuracy with which the masses and redshift are
recovered for the $(120,60)\,\rm M_{\odot}$ and $(90,90)\,\rm M_{\odot}$
binaries. The qualitative results are similar to those for the SNR 30 signals
presented in Section 3, with broader posteriors as expected due to the lower
SNR. Specifically, the masses and redshifts are poorly measured for face-on
systems, and measurement accuracy improves for inclined systems (particularly
$\iota=60\degree,90\degree$) where there is observable power in the HoMs. In
the best-case scenarios, masses and redshifts are measured with $\sim 25\%$
accuracy. For all systems other than $\iota=0\degree$ the $(120,60)\,\rm
M_{\odot}$ system is clearly identified as having unequal masses. However, the
mass distributions are broad enough that limited information about the
location of the pair-insatability mass-gap can be extracted.

For the $\iota=0\degree$ systems, and $\iota=30\degree$ for the $(90,90)\,\rm
M_{\odot}$ binary, there is a bimodality in the recovered redshift. In
addition, the inferred mass distribution is broader than that shown in Fig.13
and extends to $\sim 1000\,\rm M_{\odot}$. For these events, there is zero (or
limited) power in the HoMs so only a single GW multipole is observable. The
secondary peak at high masses and $z\approx 5$ corresponds to a binary
configuration where the (3, 3) multipole has the correct amplitude and
frequency content to match the simulated signal. This is discussed in more
detail in Section 3, around Fig. 8.

Figure 15:  Inferred posteriors for the component masses for heavy binaries.
Results are shown for inclinations of $\iota=0\degree,30\degree,60\degree$ and
$90\degree$, with the redshift of the system varied, as indicated in the
labels, so that the signals are all simulated with an SNR of 15 in the ET-CE
network, at the most sensitive sky location for the network. The grey region
denotes the pair-instability mass-gap and the darker grey denotes the
forbidden strip where no black hole is expected to form for any value of the
${}^{12}\rm C(\alpha,\gamma)^{16}O$ reaction rate (Farmer et al., 2020).
Simulated values are denoted by a black cross and contours show the 90%
credible region. Left: $m_{1}=240\,\rm M_{\odot}$ and $m_{2}=120\,\rm
M_{\odot}$. Right: $m_{1}=480\,\rm M_{\odot}$ and $m_{2}=120\,\rm M_{\odot}$.

Figure 16:  Inferred posteriors for the redshift for heavy binaries. Signals
are simulated with varying inclination and redshift, as indicated in the
labels, with a fixed SNR of 15 in the ET–CE network, at the most sensitive sky
location for the network. Black vertical lines indicate the 90% credible
intervals, violins show the 99.5% range and yellow vertical lines show
simulated values. Symbols $\rho_{\ell m}$ indicate the SNR in the $(\ell,m)$
multipole HoMs. Left: $m_{1}=240\,\rm M_{\odot}$ and $m_{2}=120\,\rm
M_{\odot}$. Right: $m_{1}=480\,\rm M_{\odot}$ and $m_{2}=120\,\rm M_{\odot}$.

In Fig. 15 and 16, we show the accuracy with which the masses and redshift are
recovered for the $(240,120)\,\rm M_{\odot}$ and $(480,120)\,\rm M_{\odot}$
binaries. As for the lower-mass systems, the qualitative results are similar
to those for the SNR 30 signals presented in Section 3, with broader
posteriors as expected due to the lower SNR. Nonetheless, other than the face-
on ($\iota=0\degree$) systems, the binaries are clearly identified as unequal
mass systems containing an IMBH with minimum mass $200\,\rm
M_{\odot}$/$400\,\rm M_{\odot}$ for the two system. Redshifts are generally
underestimated, likely due to the poor sky localization, and lower bounds on
the redshift are no better than for the higher SNR systems shown in Fig. 10.

Again, the face-on systems show significant bimodality with a second peak at
much lower redshifts and higher masses. As before, this corresponds to a
system where the HoMs, rather that the (2, 2) multipole, are associated with
the observed waveform.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, , 6, 041015, doi: 10.1103/PhysRevX.6.041015
  * Abbott et al. (2019) —. 2019, Physical Review X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, Phys. Rev. Lett., 125, 101102, doi: 10.1103/PhysRevLett.125.101102
  * Abbott et al. (2020b) —. 2020b, ApJ, 900, L13, doi: 10.3847/2041-8213/aba493
  * Abbott et al. (2020a) Abbott, R., et al. 2020a, Phys. Rev. D, 102, 043015, doi: 10.1103/PhysRevD.102.043015
  * Abbott et al. (2020b) —. 2020b, Astrophys. J. Lett., 896, L44, doi: 10.3847/2041-8213/ab960f
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., et al. 2017, ArXiv e-prints.  https://arxiv.org/abs/1702.00786
  * Arca-Sedda et al. (2021) Arca-Sedda, M., Rizzuto, F. P., Naab, T., et al. 2021, ApJ, 920, 128, doi: 10.3847/1538-4357/ac1419
  * Belczynski (2020) Belczynski, K. 2020, ApJ, 905, L15, doi: 10.3847/2041-8213/abcbf1
  * Branchesi et al. (2023) Branchesi, M., Maggiore, M., Alonso, D., et al. 2023, J. Cosmology Astropart. Phys, 2023, 068, doi: 10.1088/1475-7516/2023/07/068
  * Bromm (2013) Bromm, V. 2013, Reports on Progress in Physics, 76, 112901, doi: 10.1088/0034-4885/76/11/112901
  * Carr et al. (2021) Carr, B., Kohri, K., Sendouda, Y., & Yokoyama, J. 2021, Reports on Progress in Physics, 84, 116902, doi: 10.1088/1361-6633/ac1e31
  * Chon & Hosokawa (2019) Chon, S., & Hosokawa, T. 2019, MNRAS, 488, 2658, doi: 10.1093/mnras/stz1824
  * Costa et al. (2022) Costa, G., Ballone, A., Mapelli, M., & Bressan, A. 2022, MNRAS, 516, 1072, doi: 10.1093/mnras/stac2222
  * Costa et al. (2021) Costa, G., Bressan, A., Mapelli, M., et al. 2021, MNRAS, 501, 4514, doi: 10.1093/mnras/staa3916
  * Di Carlo et al. (2020) Di Carlo, U. N., Mapelli, M., Giacobbo, N., et al. 2020, MNRAS, 498, 495, doi: 10.1093/mnras/staa2286
  * Evans et al. (2021) Evans, M., et al. 2021.  https://arxiv.org/abs/2109.09882
  * Fairhurst (2009) Fairhurst, S. 2009, New J. Phys., 11, 123006, doi: 10.1088/1367-2630/11/12/123006
  * Fairhurst (2011) —. 2011, Class.Quant.Grav., 28, 105021, doi: 10.1088/0264-9381/28/10/105021
  * Fairhurst et al. (2023) Fairhurst, S., Hoy, C., Green, R., Mills, C., & Usman, S. A. 2023.  https://arxiv.org/abs/2304.03731
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36, doi: 10.3847/2041-8213/abbadd
  * Fishbach & Holz (2020) Fishbach, M., & Holz, D. E. 2020, ApJ, 904, L26, doi: 10.3847/2041-8213/abc827
  * Fragione & Loeb (2023) Fragione, G., & Loeb, A. 2023, ApJ, 944, 81, doi: 10.3847/1538-4357/acb34e
  * Gerosa & Fishbach (2021) Gerosa, D., & Fishbach, M. 2021, Nature Astronomy, 5, 749, doi: 10.1038/s41550-021-01398-w
  * Graziani et al. (2020) Graziani, L., Schneider, R., Marassi, S., et al. 2020, MNRAS, 495, L81, doi: 10.1093/mnrasl/slaa063
  * Green et al. (2020) Green, R., Hoy, C., Fairhurst, S., et al. 2020.  https://arxiv.org/abs/2010.04131
  * Hannam et al. (2013) Hannam, M., Brown, D. A., Fairhurst, S., Fryer, C. L., & Harry, I. W. 2013, ApJ, 766, L14, doi: 10.1088/2041-8205/766/1/L14
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27, doi: 10.1146/annurev-astro-120419-014455
  * Kalogera et al. (2019) Kalogera, V., Berry, C. P. L., Colpi, M., et al. 2019, BAAS, 51, 242.  https://arxiv.org/abs/1903.09220
  * Klessen & Glover (2023) Klessen, R. S., & Glover, S. C. O. 2023, arXiv e-prints, arXiv:2303.12500, doi: 10.48550/arXiv.2303.12500
  * Liu & Bromm (2020) Liu, B., & Bromm, V. 2020, ApJ, 903, L40, doi: 10.3847/2041-8213/abc552
  * Madau et al. (2014) Madau, P., Haardt, F., & Dotti, M. 2014, ApJ, 784, L38, doi: 10.1088/2041-8205/784/2/L38
  * Maggiore et al. (2020) Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, J. Cosmology Astropart. Phys, 2020, 050, doi: 10.1088/1475-7516/2020/03/050
  * Maiolino et al. (2023) Maiolino, R., Scholtz, J., Curtis-Lake, E., et al. 2023, arXiv e-prints, arXiv:2308.01230, doi: 10.48550/arXiv.2308.01230
  * Mancarella et al. (2023) Mancarella, M., Iacovelli, F., & Gerosa, D. 2023, Phys. Rev. D, 107, L101302, doi: 10.1103/PhysRevD.107.L101302
  * Mandel & Farmer (2022) Mandel, I., & Farmer, A. 2022, Phys. Rep., 955, 1, doi: 10.1016/j.physrep.2022.01.003
  * Mapelli (2021) Mapelli, M. 2021, in Handbook of Gravitational Wave Astronomy, 16, doi: 10.1007/978-981-15-4702-7_16-1
  * Mapelli et al. (2022a) Mapelli, M., Bouffanais, Y., Santoliquido, F., Arca Sedda, M., & Artale, M. C. 2022a, MNRAS, 511, 5797, doi: 10.1093/mnras/stac422
  * Mapelli et al. (2022b) —. 2022b, MNRAS, 511, 5797, doi: 10.1093/mnras/stac422
  * Mapelli et al. (2021) Mapelli, M., Dall’Amico, M., Bouffanais, Y., et al. 2021, MNRAS, 505, 339, doi: 10.1093/mnras/stab1334
  * Marchant & Moriya (2020) Marchant, P., & Moriya, T. J. 2020, A&A, 640, L18, doi: 10.1051/0004-6361/202038902
  * Mills & Fairhurst (2021) Mills, C., & Fairhurst, S. 2021, Phys. Rev. D, 103, 024042, doi: 10.1103/PhysRevD.103.024042
  * Mills et al. (2018) Mills, C., Tiwari, V., & Fairhurst, S. 2018, Physical Review D, 97, 104064, doi: 10.1103/PhysRevD.97.104064
  * Ng et al. (2022) Ng, K. K. Y., Chen, S., Goncharov, B., et al. 2022, ApJ, 931, L12, doi: 10.3847/2041-8213/ac6bea
  * Ng et al. (2023) Ng, K. K. Y., Goncharov, B., Chen, S., et al. 2023, Phys. Rev. D, 107, 024041, doi: 10.1103/PhysRevD.107.024041
  * Nitz & Capano (2021) Nitz, A. H., & Capano, C. D. 2021, ApJ, 907, L9, doi: 10.3847/2041-8213/abccc5
  * Pezzulli et al. (2016) Pezzulli, E., Valiante, R., & Schneider, R. 2016, MNRAS, 458, 3047, doi: 10.1093/mnras/stw505
  * Pratten et al. (2020) Pratten, G., et al. 2020.  https://arxiv.org/abs/2004.06503
  * Punturo et al. (2014) Punturo, M., Lück, H., & Beker, M. 2014, in Astrophysics and Space Science Library, Vol. 404, Advanced Interferometers and the Search for Gravitational Waves, ed. M. Bassan, 333, doi: 10.1007/978-3-319-03792-9_13
  * Renzo et al. (2020) Renzo, M., Cantiello, M., Metzger, B. D., & Jiang, Y. F. 2020, ApJ, 904, L13, doi: 10.3847/2041-8213/abc6a6
  * Ricarte & Natarajan (2018) Ricarte, A., & Natarajan, P. 2018, MNRAS, 481, 3278, doi: 10.1093/mnras/sty2448
  * Roupas & Kazanas (2019) Roupas, Z., & Kazanas, D. 2019, A&A, 632, L8, doi: 10.1051/0004-6361/201937002
  * Safarzadeh & Haiman (2020) Safarzadeh, M., & Haiman, Z. 2020, ApJ, 903, L21, doi: 10.3847/2041-8213/abc253
  * Schneider et al. (2023) Schneider, R., Valiante, R., Trinca, A., et al. 2023, arXiv e-prints, arXiv:2305.12504, doi: 10.48550/arXiv.2305.12504
  * Singer et al. (2014) Singer, L. P., Price, L. R., Farr, B., et al. 2014.  https://arxiv.org/abs/1404.5623
  * Spera & Mapelli (2017) Spera, M., & Mapelli, M. 2017, ArXiv e-prints.  https://arxiv.org/abs/1706.06109
  * Stacy & Bromm (2013) Stacy, A., & Bromm, V. 2013, MNRAS, 433, 1094, doi: 10.1093/mnras/stt789
  * Sugimura et al. (2023) Sugimura, K., Matsumoto, T., Hosokawa, T., Hirano, S., & Omukai, K. 2023, arXiv e-prints, arXiv:2307.15108, doi: 10.48550/arXiv.2307.15108
  * Tanikawa et al. (2021) Tanikawa, A., Kinugawa, T., Yoshida, T., Hijikawa, K., & Umeda, H. 2021, MNRAS, 505, 2170, doi: 10.1093/mnras/stab1421
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2021a, arXiv e-prints, arXiv:2108.01045.  https://arxiv.org/abs/2108.01045
  * The LIGO Scientific Collaboration et al. (2021b) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021b, arXiv e-prints, arXiv:2111.03606.  https://arxiv.org/abs/2111.03606
  * Trinca et al. (2022) Trinca, A., Schneider, R., Maiolino, R., et al. 2022, arXiv e-prints, arXiv:2211.01389.  https://arxiv.org/abs/2211.01389
  * Usman et al. (2019) Usman, S. A., Mills, J. C., & Fairhurst, S. 2019, Astrophys. J., 877, 82, doi: 10.3847/1538-4357/ab0b3e
  * Valiante et al. (2017) Valiante, R., Agarwal, B., Habouzit, M., & Pezzulli, E. 2017, PASA, 34, e031, doi: 10.1017/pasa.2017.25
  * Valiante et al. (2021) Valiante, R., Colpi, M., Schneider, R., et al. 2021, MNRAS, 500, 4095, doi: 10.1093/mnras/staa3395
  * van Son et al. (2020) van Son, L. A. C., De Mink, S. E., Broekgaarden, F. S., et al. 2020, ApJ, 897, 100, doi: 10.3847/1538-4357/ab9809
  * Veitch et al. (2015) Veitch, J., Raymond, V., Farr, B., et al. 2015, Phys. Rev. D, D91, 042003, doi: 10.1103/PhysRevD.91.042003
  * Vink et al. (2021) Vink, J. S., Higgins, E. R., Sander, A. A. C., & Sabhahit, G. N. 2021, MNRAS, 504, 146, doi: 10.1093/mnras/stab842
  * Vitale (2016) Vitale, S. 2016, Phys. Rev. D, 94, 121501, doi: 10.1103/PhysRevD.94.121501
  * Volonteri et al. (2021) Volonteri, M., Habouzit, M., & Colpi, M. 2021, Nature Reviews Physics, 3, 732, doi: 10.1038/s42254-021-00364-9
  * Volonteri & Rees (2005) Volonteri, M., & Rees, M. J. 2005, ApJ, 633, 624, doi: 10.1086/466521
  * Volonteri et al. (2015) Volonteri, M., Silk, J., & Dubus, G. 2015, ApJ, 804, 148, doi: 10.1088/0004-637X/804/2/148
  * Woosley & Heger (2021) Woosley, S. E., & Heger, A. 2021, ApJ, 912, L31, doi: 10.3847/2041-8213/abf2c4

  *[LVK]: LIGO-Virgo-KAGRA
  *[GW]: gravitational-wave
  *[AGN]: active galactic nucleus
  *[BBHs]: binary black hole
  *[IMF]: initial mass function
  *[GWs]: gravitational-wave
  *[BBH]: binary black hole
  *[BHs]: black hole
  *[HoMs]: Higher order Multipole
  *[SNR]: signal-to-noise ratio
  *[ET]: The Einstein Telescope
  *[SNRs]: signal-to-noise ratio
  *[CE]: Cosmic Explorer
  *[IMBH]: intermediate mass black hole
  *[BH]: black hole
  *[PISN]: Pair Instability Supernova
  *[IMBHs]: intermediate mass black hole
  *[HoM]: Higher order Multipole