# The impact of a primordial gravitational wave background on LISA resolvable
sources

Matteo Braglia<EMAIL_ADDRESS>Center for Cosmology and Particle Physics, New
York University, 726 Broadway, New York, NY 10003, USA  Mauro Pieroni
<EMAIL_ADDRESS>Theoretical Physics Department, CERN, 1211 Geneva 23,
Switzerland  Sylvain Marsat<EMAIL_ADDRESS>Laboratoire des 2
Infinis - Toulouse (L2IT-IN2P3), Université de Toulouse, CNRS, UPS, F-31062
Toulouse Cedex 9, France

###### Abstract

The existence of a primordial stochastic gravitational wave background (SGWB)
is a common prediction in various models of the early Universe. Despite
constraints at different frequency ranges and claims of detection in the nHz
range by Pulsar Timing Arrays, the amplitude and spectral dependence of the
SGWB in the mHz range remain largely unknown. Plausible models of early
Universe Physics predict a wide range of SGWB amplitudes, from undetectable to
exceeding the constraints from Big Bang Nucleosynthesis. This paper explores
the potential impact of a prominent primordial SGWB on LISA’s main scientific
targets. Our main analyses focuses on Massive Black Hole Binaries (MBHBs). By
employing publicly available MBHB population models and state-of-the-art
LISA’s forecasting pipeline, we analyze the effects of the SGWB on MBHB
detections. We find that the decrease in the signal-to-noise ratio induced by
a strong primordial GWB can significantly reduce the number of detectable
events, compromise the precision of constraints, and even hinder sky
localization for some events. We also examine the impact of the SGWB on
Stellar Origin Black Hole Binaries (SOBHBs) and Galactic Binaries (GBs), which
are fainter sources than MBHBs. Depending on the spectral properties of the
SGWB, we conclude that these sources could be either marginally affected or
rendered completely undetectable. This largely unexplored aspect raises
critical questions about the potential challenges posed by a prominent SGWB to
LISA’s astrophysical objectives, including MBHBs, SOBHBs, and GBs.

Suggested keywords

###### pacs:

Valid PACS appear here

††preprint: CERN-TH-2024-073

##  I Introduction

The Laser Interferometer Space Antenna (LISA) [1] represents a pioneering
effort in the realm of gravitational wave (GW) astronomy, slated for launch in
the mid-2030s. Designed to operate within the mHz frequency band, LISA’s
primary mission is to detect and study the merging of Massive Black Hole
Binaries (MBHBs) in the distant Universe. However, its reach extends beyond
this focus, encompassing a diverse array of GW sources including stellar
remnant or compact object binary systems, Double White Dwarfs (DWDs) within
our galaxy and the inspiral stage of Stellar Origin Black Hole Binaries
(SOBHBs), ultimately merging into the sensitivity band of current and future
ground-based detectors. See, e.g., [2], for a review of the astrophysical
sources LISA will probe.

The significance of these detections cannot be overstated. They promise
profound insights into fundamental questions about our Universe, enabling
tests of gravity in uncharted regimes [3, 4, 5, 6], exploration of the
environments surrounding black holes [7], investigations into galaxy evolution
[8, 9, 10, 11], and even larger-scale inquiries into the expansion of our
Universe [12, 13, 14, 15, 16]. Drawing upon the expertise of researchers from
diverse fields within astrophysics and cosmology, these ambitious scientific
endeavors hinge upon the accurate detection and parameter estimation of
gravitational wave signals.

Current assessments of LISA’s ability to detect these sources rely on noise
curves derived from various factors, including the constellation’s arm length,
non-inertial perturbations of geodesic references, and interferometric
measurement noise. As documented in several definition study reports [17, 18],
ESA, in collaboration with the LISA consortium, continuously refines its
understanding of these elements, contributing to a relatively robust estimate
of LISA’s sensitivity. In this paper we investigate the impact on GW parameter
estimation of an additional, theoretically uncontrollable noise source: the
presence of a large stochastic gravitational wave background (SGWB).

Predictions for SGWBs of astrophysical origin, arising from the superposition
of many subthreshold galactic [19, 20, 21], or extragalactic binaries [22],
are reasonably well-understood and often considered in the estimation of
LISA’s scientific capabilities. However, the same level of understanding does
not extend to SGWBs of primordial origin, potentially leading to significant
theoretical uncertainty111Extreme mass ratio inspirals also produce another
(astrophysical) SGWB which is very often neglected [23, 24, 25]. However, even
within the theoretical uncertainties, the amplitude of this SGWB is expected
to be significantly smaller than that of the primordial SGWBs considered in
this work.. Primordial SGWBs can originate shortly after the Big Bang through
various types of early Universe physics, including inflation, phase
transitions, or cosmic strings [26, 27, 28, 29, 30, 31, 32, 33]. The amplitude
of the generated SGWBs can span from undetectable levels to extremely large,
and their spectral shape varies significantly depending on the underlying
mechanism. Currently, the only established constraints on the SGWB amplitude
in the mHz frequency band come from observations of the Cosmic Microwave
Background (CMB) and Big Bang Nucleosynthesis (BBN), limiting the (integrated
over frequency) amplitude of the SGWB’s energy density to be
$h^{2}\,\Omega_{\rm GW}\,<\,1.2\times 10^{-6}$ at $95$% CL [34]. Well-
motivated models of the early universe could potentially produce anything
below this upper bound. However, this constraint is not highly stringent. In
fact, when translated into power spectral density, this limit is nearly 2
orders of magnitude above the nominal sensitivity of LISA at certain
frequencies.

Motivated by these uncertainties, as well as the recent detection of an SGWB
by Pulsar Timing Array collaborations [35, 36, 37, 38, 39, 40, 41, 42, 43] at
smaller frequencies, which indicate that SGWB may indeed play a prominent role
in uncharted frequency bands of GW astronomy, for the first time, we examine
quantitatively the effects of such an SGWB on LISA sources. We begin in
Section II by extending the well-known match-filtered parameter estimation
methodology to account for the presence of an SGWB. Subsequently, we conduct
quantitative forecasts on parameter estimation for the most relevant sources
in the LISA band in the presence of an SGWB, comparing them to the ideal, and
well studied, scenario of an absent SGWB. To accomplish this, we utilize the
state-of-the-art forecasting pipeline for LISA, as implemented in the lisabeta
code. We provide an in-depth analysis of the science case for MBHBs in Section
III, presenting results for several mock populations simulated based on the
state-of-the-art catalogs of MBHBs. Our findings suggest that parameter
estimation could be significantly impacted, leading to fewer detected sources
if the SGWB is sufficiently loud. In Sections IV.1 and IV.2, we draw more
generalized conclusions about SOBHB and Galactic Binaries (GBs), demonstrating
a degradation in scientific capabilities in these scenarios as well,
potentially even worse than for MBHBs.

This work focuses on the most extreme scenario where the SGWB saturates the
CMB+BBN bound. As such our results represent the most pessimistic scenario for
what concerns the detection of individual GW signals with LISA and provide an
upper bound on the degradation due to the unexpected presence of an SGWB.
Nonetheless, given that LISA will target a completely unexplored frequency
range, we must remain open to unexpected phenomena. We believe that our
findings underscore the importance of considering the potential impact of an
SGWB seriously and warrant further study in this area.

##  II Methodology

In this section, we outline our forecast methodology, beginning with a review
of widely used data analysis techniques. We then discuss the incorporation of
the primordial SGWB as an additional source of noise for LISA, along with
illustrations of representative primordial signals consistent with existing
experimental constraints.

Data and likelihood. Let us start by considering the case where no SGWB is
present in the data. Moreover, let us work under the idealized assumption that
each data stream has been cleaned of glitches and other instrumental
artifacts. Denoting the time-domain data stream for each of the 3 LISA
channels as a combination of signal and noise, $d_{i}(t)=s_{i}(t)+n_{i}(t)$,
where $i$ ranges over the orthogonal combinations222Assuming the LISA
configuration to be a rigid equilateral triangle, with the same noise in all
links, it can be shown that the signal and noise are simultaneously diagonal
in such basis, which is defined through a frequency-independent rotation
starting from the Michelson XYZ basis [44]. Moreover, under these assumptions,
the T channel is nearly signal-insensitive at low frequencies, making it a
null channel that can be used to monitor the noise. $A$, $E$, and $T$ (see,
e.g., Ref. [45]), we can express the data in the frequency domain by taking a
Fourier transform of the time-domain data as:

$\tilde{d}_{i}(f)=\int_{-T/2}^{T/2}dt\,e^{2\pi ift}d_{i}(t).$ (1)

In the following, we assume the noise to be Gaussian and stationary with the
ensemble variance given by:

$\langle\tilde{n}_{i}(f)\tilde{n}_{j}^{*}(f^{\prime})\rangle=\frac{1}{2}\delta(f-f^{\prime})\delta_{ij}S_{n,\,i}(f),$
(2)

where $\tilde{n}_{i}(f)$ denotes the noise in channel $i$ in the frequency
domain, and $S_{n,\,i}(f)$ is the single-sided noise power spectral density
(PSD) for TDI channel $i$. Details regarding noise PSDs for each channel can
be found, e.g., in Ref. [46].

While the noise is a stochastic variable with zero mean, which can only be
characterized through its variance, the signal $s_{i}(t)$ is a deterministic
variable. Moreover, $s_{i}(t)$ is a combination of the “true” GW signal,
typically expressed in terms of a waveform 333It is important to note that our
template waveform, $h(t,\bm{\theta})$, serves as an approximation of the
physical signal. In the subsequent analysis, we make the simplifying
assumption that our template waveforms perfectly match the astrophysical
waveforms. See e.g. [47, 48] for studies on this topic. $h(t,\bm{\theta}_{\rm
GW})$ depending on some physical parameters $\bm{\theta}_{\rm GW}$
characterizing the system, modulated onto each of the data streams by the
detector response function $r_{i}(t,\bm{\theta}_{\rm r})$ with
$\bm{\theta}_{\rm r}$ collectively denoting the sky-localization and the
orientation of the source (for details on the computation of the LISA response
function, see, e.g., Ref [49, 50, 46, 51]). In the following, we will
collectively denote the signal parameters with $\bm{\theta}$, and the signal
with $s_{i}(\bm{\theta})$. In our analysis, we employ IMRPhenomXHM as a
waveform approximant [52, 53].

An optimal estimator for the signal can be defined by taking advantage of the
difference in the statistical properties of the signal and noise. Typically,
this is done using the standard matched-filter inner product [54]:

$(a|b)=2\int_{f_{1}}^{f_{2}}\frac{a(f)b^{*}(f)+a^{*}(f)b(f)}{S_{n}(f)}\,\textrm{d}f,$
(3)

through which the log-likelihood can be expressed as:

$\ln\mathcal{L}=-\sum_{i=A,\,E,\,T}\frac{1}{2}\left(d_{i}-s_{i}(\bm{\theta})|d_{i}-s_{i}(\bm{\theta})\right),$
(4)

i.e., as a Gaussian log-likelihood for the noise, where three matched-filter
inner products are evaluated with the PSD of the three TDI channels.

Parameter estimation: posterior sampling and Fisher formalism. We employ two
methodologies to forecast constraints on the parameters $\bm{\theta}$ starting
from the likelihood above.

The first approach involves calculating the posterior distribution of the
model parameters by directly sampling the likelihood. We obtain the full
posterior distribution by multiplying the likelihood by the prior on the
parameters using Bayes theorem $p\left(\bm{\theta}|d\right)\propto
p(d|\bm{\theta})p(\bm{\theta})$, where we have dropped the evidence, as we
will not be concerned with issues related to model selection. We then estimate
errors on each parameter by marginalizing over the others, as customary. In
practice, we calculate the likelihood utilizing the lisabeta code introduced
in Ref. [55]. To produce samples from the posterior, we use ptemcee [56], an
ensemble sampler [57, 58] with adaptive parallel tempering.

While sampling the posterior provides a fully Bayesian procedure to calculate
any posterior distribution, regardless of its structure, it is in practice
computationally consuming, and it would become quickly too demanding to scale
the forecast to, e.g., a whole population of MBHBs. For this reason, in order
to quickly estimate the uncertainty in the parameter reconstruction, we also
pursue an alternative route and employ a simpler Fisher formalism. Defining
the Fisher Information Matrix (FIM) as

$F_{lk}\equiv\left(\frac{\partial\tilde{h}^{\rm
th}(f,\bm{\theta})}{\partial\theta_{l}}\Biggr{\rvert}\frac{\partial\tilde{h}^{\rm
th}(f,\bm{\theta})}{\partial\theta_{l}}\right)\Biggr{\rvert}_{\bm{\theta}=\bm{\theta}_{0}}\;.$
(5)

we can approximate the likelihood in the vicinity of its maximum
$\bm{\theta}_{0}$ as

$\ln\mathcal{L}\simeq-\frac{1}{2}F_{lk}\,\Delta\theta_{l}\Delta\theta_{k},$
(6)

where $\Delta\theta_{i}$ is the deviation of a parameter from its true value.
We can then invert the FIM to compute the covariance matrix
$\Sigma_{kl}=F_{ij}^{-1}$. As customary, the confidence intervals on the model
parameters $\bm{\theta}$ can be drawn from the covariance matrix
$\Sigma_{lk}$. In particular, we estimate the (marginalized) error on each
parameter using the Kramer-Rao bound
$\sigma(\theta_{i})\geq\sqrt{\Sigma_{ii}}$.

Let us emphasize that Fisher forecasts rely on the assumption of a Gaussian
likelihood in the parameters $\bm{\theta}$, which is not always the case,
especially when degeneracies are present. In general, this condition is more
likely to be met if the signal is sufficiently loud [59].

Inclusion of an SGWB. Let us proceed by discussing the inclusion of an SGWB in
the data stream and how it would affect the measurement of the individual
source parameters. The key point is to realize that, for individual event
characterization, an SGWB can effectively be interpreted as an additional
noise source that, due to the response function, affects the different TDI
channels in different ways. In practice, for what concerns the
characterization of signals from individual sources, in any of the TDI
channels, the instrumental noise is replaced by 444See also Refs. [60, 61] for
a similar approach.:

$\tilde{n}_{i}\mapsto\tilde{n}^{\rm
TOT}_{i}=\tilde{n}_{i}+\tilde{\sigma}_{i},$ (7)

where $\tilde{\sigma}_{i}$ is the SGWB contribution in the TDI channel $i$,
which is uncorrelated with $\tilde{n}_{i}$. Similarly to the noise, in the AET
basis, the SGWB satisfies

$\langle\tilde{\sigma}_{i}(f)\tilde{\sigma}_{j}^{*}(f^{\prime})\rangle\equiv\frac{1}{2}\delta\left(f-f^{\prime}\right)\delta_{ij}S_{{\rm
SGWB,i}}\left(f\right)\;,$ (8)

where $S_{{\rm SGWB,i}}$ is a combination of the SGWB PSD and the (sky average
of the square of) detector response function555For details on the projection
of the SGWB onto the data stream see, e.g., Refs. [62, 32, 33]. Using the
formula

$h^{2}\Omega_{\textrm{GW}}(f)\equiv\frac{4\pi^{2}f^{3}}{3H_{0}^{2}}S_{{\rm
SGWB}}(f)\;,$ (9)

where $H_{0}$ is the present value of the Hubble parameter, we can relate the
PSD of the SGWB to its energy density, which is often easier to calculate from
a given theoretical model. With these results, it is easy to extend our
likelihood to include the effect of an SGWB by simply replacing

$S_{n,\,i}\mapsto S_{n,\,i}+S_{{\rm SGWB}}$ (10)

in the likelihood (4). In the following, we will refer to “No SGWB” as the
case where $S_{{\rm SGWB}}$ only contains the contribution from the galactic
astrophysical SGWB, sometimes also referred to as foreground, modeled
according to Ref. [20].

Figure 1: Power spectral density $S_{n}$ with and without a SGWB. For
comparison, we also plot illustrative signals including MBHBs for different
total masses and redshifts, a GW150914-like event [63] and galactic and
verification binaries. The latter are taken from Ref. [64, 65], and we only
plot those with SNR$>5$. For MBHBs and GW150914, we only plot the +
polarization of GWs for definiteness. The parameters for the SGWBs are
reported in the main text and are chosen to saturate the CMB+BBN bound. All
the noise curves also include galactic foregrounds.

SGWB models. We close this section by presenting the modeling of the
cosmological SGWBs that we use in our forecast. The only current experimental
limits on the amplitude of the SGWB in the mHz range are set by BBN and CMB
data, which constrain the energy density of GWs at early times, through their
contribution to the number of extra relativistic degrees of freedom $\Delta
N_{\rm eff}$. The most updated limit is [34]:

$h^{2}\,\int_{f_{\rm min}}^{f_{\rm max}}\,\mathrm{d}\ln f\,\Omega_{\rm
GW}(f)\,<\,1.2\times 10^{-6}\,\,\,\,{\rm at}\,\,\,\,95\%{\rm CL},$ (11)

where we take $f_{\rm min}=10^{-5}\,{\rm Hz}$ and $f_{\rm max}=0.5\,{\rm Hz}$.
For definiteness, we choose our fiducials to saturate the BBN+CMB limit (11).
As such, the outcome of our analysis in the following Sections has to be
interpreted as the most pessimistic scenario for detecting deterministic GW
signals with LISA.

In this paper, we remain agnostic to the specific kind of production mechanism
generating the SGWB, and we consider the following phenomenological templates:

$\displaystyle\Omega_{\rm GW}^{\rm PL}$ $\displaystyle=$ $\displaystyle A_{\rm
PL}\left(\frac{f}{f_{*}}\right)^{\alpha}\;,$ (12) $\displaystyle\Omega_{\rm
GW}^{\rm Bump}$ $\displaystyle=$ $\displaystyle\frac{A_{\rm
Bump}}{\sqrt{2\pi}\Delta}\exp\left[-\frac{1}{2\Delta^{2}}\log^{2}\frac{f}{f_{*}}\right]\;.$
(13)

The Power Law template (12) illustrates the effect of a broad SGWB, degrading
the sensitivity of LISA over a large range of frequencies. On the other hand,
the Bump template (13) is a proxy for any SGWB with sufficiently limited
frequency support, and affects the sensitivity in a narrower range of
frequencies.

Specifically, we will consider the following benchmarks, whose impact on the
mission’s sensitivity is shown in Fig. 1, where we also plot the strain of
some typical events that are targeted by LISA:

  * •

PL. $A_{\rm PL}=1.27\times 10^{-7},\,\alpha=0$. This case illustrates a
background spanning several decades in frequency within the LISA band. Note
that while the energy density is flat $\Omega_{\rm GW}$ for $\alpha=0$, its
power spectral density scales as $f^{-3}$. We expect it to impact the
detection and parameter estimation of MBHB spanning a large range of masses,
as well as the detection of GBs and SOBHBs.

  * •

Bump Low $f$. $A_{\rm Bump}=1.27\times 10^{-7},\,\Delta=0.5,\,f_{*}=2.2\times
10^{-3}$ Hz. We choose a small width $\Delta$ in order for the SGWB to be
reasonably peaked. Our choice of $f_{*}$ roughly corresponds to the ISCO
frequency of binaries with $M\sim 10^{7}\,M_{\odot}$. As can be seen from Fig.
1, the SGWB peaks around the frequencies where GBs are expected to be detected
and has an amplitude comparable to or even larger than their signals. However,
since the SGWB signal dies off at higher frequencies, we do not expect it to
have a significant effect on SOBHB detections.

  * •

Bump High $f$. $A_{\rm Bump}=1.27\times 10^{-7},\,\Delta=0.5,\,f_{*}=2.2\times
10^{-2}$ Hz. Same as the previous benchmark, but now the choice of $f_{*}$
roughly corresponds to the ISCO frequency of binaries with $M\sim
10^{6}\,M_{\odot}$. In addition to that, and unlike Bump Low $f$, this signal
is expected to impact measurements of SOBHs, leaving almost unaltered
detections of GBs.

Figure 2: [First column.] We plot the ${\rm SNR}$ distribution averaged over
100 realizations of our catalogs of MBHBs populations. We emphasize the point
${\rm SNR}=10$ with a vertical black line and report the number of events
detected with ${\rm SNR}>10$ in the legend. [Second, third, and fourth
columns.] For a specific realization of our populations, we plot the
distribution of events as a function of redshift and the total mass of the
binary. The color bar indicates the ${\rm SNR}$ of each event, and we
emphasize events with ${\rm SNR}>10$ with a colored circle. We plot results
for PopIII , Q3delays , and Q3nodelays populations in the first, second, and
third rows respectively. Blue, orange, green, and red colors represent our
results without SGWB, for PL, Bump High $f$, and Bump Low $f$ respectively.

Figure 3:  Relative 68% CL errors on the total mass $M$ [top] and the
luminosity distance $d_{L}$ [center] and absolute 68% CL errors on the sky
localization $\Omega$ [bottom] from our FIM analysis. We color the region
where the percentage error is $>100\%$ with black. We only perform the Fisher
analysis for the events that are detected with ${\rm SNR}>10$ in the absence
of an SGWB. The errors are shown as function of the fiducial values in the
scatter plots, and we also provide histograms to visualize their distribution.
Our results are averaged over 15 random realizations of our populations. We
plot results for PopIII , Q3delays , and Q3nodelays populations in the first,
second, and third rows respectively. Blue, orange, green, and red colors
represent our results without SGWB, for PL, Bump High $f$, and Bump Low $f$
respectively.

##  III Massive Black Hole Binaries

As highlighted in the introduction, the primary focus of our analysis is on
MBHBs, aligning with one of the key science objectives of LISA. By measuring
hundreds of MBHBs, LISA will contribute significantly to answering questions
about the formation and evolution of these systems at high redshifts, a topic
that remains a subject of ongoing investigation. In this section, our goal is
to assess the impact of a significant SGWB on parameter estimation for
individual MBHBs and extend this analysis to their populations.

For clarity in our subsequent discussions, we introduce the signal-to-noise
ratio (SNR) of a GW event as ${\rm SNR}^{2}\equiv(s|s)$. Here, the scalar
product is computed using the noise power spectral density (PSD) defined in
Eq. (10), which accommodates the inclusion of the SGWB. Additionally, we
define the quantity $\bar{N}_{\rm events}$ as the number of events with an SNR
greater than or equal to 10, which we establish as a threshold for detecting a
GW event.

###  III.1 Population analysis

Let us begin by introducing our populations of MBHBs. We adopt the populations
presented in Ref. [66], based on the semi-analytical model described in Ref.
[67]. We note that slightly more updated predictions for MBHB populations
exist (see, e.g., [68, 69] and references therein). However, since the data
for the populations in [66] are publicly
available666https://people.sissa.it/barausse/catalogs, we employ those for our
analysis. Moreover, while results would quantitatively change for different
populations, the qualitative points raised by our work would still hold.
Following the methodology of Ref. [66], we generate three distinct populations
that primarily differ in the initial MBH seeds and the time delay between MBH
and galaxy mergers. For each of these populations, we illustrate in Fig. 2 the
distribution of the SNR of the generated events, along with their
distributions as functions of redshift and (redshifted) total binary mass. To
provide insight into the effects of introducing an SGWB, we plot these
quantities for each SGWB benchmark discussed in the previous section.

Next, we delve into the features of these generated populations and explore
the impact of an SGWB on the detectability of their sources:

  * •

PopIII : This population is based on the assumption of light MBH seeds
resulting from the explosion of PopIII stars and incorporates delays between
MBH and galaxy mergers, making the model realistic and conservative. Notably,
even if the delays were set to zero, the event rate would only change by
approximately a factor of two [66]. When averaging over 100 realizations of
our catalogs, we identify roughly 140 detectable events with ${\rm SNR}\geq
10$. As depicted in the upper panels of Fig. 2, this population predicts a
high merger rate extending up to $z\sim 20$, with many detectable events
occurring at high redshifts. Arising from light seeds, most of these events
involve binaries with masses ranging between $M_{z}\sim
10^{4}-10^{5}M_{\odot}$, while only a few dozen have larger masses of around
$M_{z}\sim 10^{6}-10^{8}M_{\odot}$. The SGWB exerts a substantial impact
across all three benchmarks, with the most pronounced effect observed for the
PL case, where only approximately one-seventh of the events remain detectable,
predominantly with larger masses. In contrast, in both the Bump Low and High
$f$ scenarios, nearly double the events remain detectable, albeit with a more
pronounced reduction in SNR around masses $M_{z}\sim 10^{7}M_{\odot}$ and
$M_{z}\sim 10^{5}M_{\odot}$ respectively compared to the PL case, consistent
with the discussion in the previous section.

  * •

Q3delays : This population is distinguished by its heavy MBH seeds originating
from the collapse of protogalactic disks. Similar to PopIII , it accounts for
the delays between MBH and galaxy mergers, lending it a realistic and
conservative nature. The name ”Q3” is derived from the so-called critical
Toomre parameter $Q_{c}$, which quantifies the stability of protogalactic
disks and the likelihood that a halo hosts a black hole seed. In this model,
$Q_{c}$ is set to 3 [66], although reasonable values typically range from 1.5
to 3 (see, e.g., Ref. [70]). The mass and redshift distribution of events in
this model differs from that of PopIII , with all events occurring at
redshifts $z\leq 10$ and involving large masses $M_{z}\gtrsim
10^{6}M_{\odot}$. We identified approximately 30 detectable events for this
model, a number significantly lower than that for PopIII . However, their SNRs
are consistently high, with nearly all events remaining detectable even in the
presence of a strong SGWB, although the SNR may degrade by a large extent—up
to an order of magnitude in the PL case.

  * •

Q3nodelays : This population mirrors Q3delays , but without accounting for
delays between galaxy and MBHB mergers, presenting a more optimistic scenario
for LISA’s MBHB merger rates. The event distribution extends to high
redshifts, reaching $z\sim 20$, similar to PopIII , and involves mainly large
masses $M_{z}\gtrsim 10^{6}M_{\odot}$, akin to Q3delays , but with some events
ranging from $M_{z}\sim 10^{5}M_{\odot}$ to $M_{z}\sim 10^{8}M_{\odot}$. We
identify nearly 500 detectable events, which is over three times the number
for PopIII and ten times that of Q3delays . Similar to Q3delays ,
$\bar{N}_{\rm events}$ does not change dramatically, although the average SNR
is reduced by up to an order of magnitude.

The distribution of events within each population heavily influences the
impact of the SGWB on the detectability of MBHBs. Populations with a
significant fraction of events at masses $M_{z}\lesssim
10^{5}-10^{6}M_{\odot}$ (like PopIII in our study) experience a notable
reduction in detectable events. This reduction could significantly impact the
reconstruction of the mass distribution of events. Conversely, populations
with a predominant fraction of events at higher masses, such as Q3delays and
Q3nodelays , fare better under the influence of the SGWB.

To quantify these effects more rigorously, we extend our analysis beyond
simple SNR considerations to assess the degradation of parameter estimation
due to the SGWB. We employ the Fisher formalism to derive constraints on three
key parameters: the total mass $M$, the luminosity distance $d_{L}$, and the
sky localization $\Delta\Omega$. Our results are depicted in Fig. 3. These
figures compare the (marginalized) error distributions for events with ${\rm
SNR}\geq 10$ in the absence and presence of our three benchmark SGWBs. Regions
where the forecast relative error exceeds $100\%$ indicate non-detection. For
sky localization, we use a non-detectability criterion of
$\log_{10}\Delta\Omega\,[{\rm deg}^{2}]\geq 2$. Our results are based on
averaging over 15 realizations, to average over Poisson fluctuations and
better represent the actual event distribution.

From Fig. 3, the error distribution for the total mass $M$ shifts by an order
of magnitude towards larger errors for all three populations. Notably, errors
for Q3delays and Q3nodelays populations remain relatively low compared to
PopIII , where errors can reach up to $10\%$.

In Fig. 3, the impact of the SGWB is even more pronounced, particularly for
the PopIII population. We observe a degradation of the relative error by
approximately one order of magnitude. While this poses no significant threat
to populations dominated by high-mass events, it severely impacts populations
like PopIII where low-mass events are prevalent, pushing many events into non-
detection territory.

Turning to the sky localization $\Delta\Omega$, whose forecast uncertainties
are shown in he bottom panels of Fig. 3, we find that constraints on this
parameter are also significantly degraded by the presence of an SGWB. Even the
Q3delays and Q3nodelays populations are affected. For Q3nodelays , nearly half
of the events with ${\rm SNR}\geq 10$ are impacted for the Bump Low and High
$f$ cases, and even more so for the PL case.

While our Fisher formalism provides a quantitative estimate of these impacts,
its assumptions may break down when approaching the detection threshold, as
indicated by the red regions in Fig. 3. Despite this limitation, we expect our
estimates to provide a reasonable order of magnitude for the errors and an
accurate description of their distribution across populations.

Figure 4: Results of simulated Bayesian parameter estimation for the MBHB 1 [left panels] and MBHB 2 [right panels] systems of Table 1. Shown is the posterior distribution on intrinsic parameters, masses, and spins [upper panels] and on the two angles for sky localization, as well as luminosity distance [lower panels]. Contours are shown at 90 % CL. Blue, orange, green, and red colors represent our results without SGWB, for PL, Bump High $f$, and Bump Low $f$ respectively. | MBHB 1 | MBHB 2
---|---|---
$M$ | $10^{5}M_{\odot}$ | $3\times 10^{6}M_{\odot}$
$q$ | $2.5$ | $1.2$
$\chi_{1}$ | $0.95$ | $0.9$
$\chi_{2}$ | $0.95$ | $0.95$
$z$ | $3.$ | $3.$
$t_{0}$ | $0.75$ | $0.65$
$\iota$ | $2.4$ | $2.8$
$\varphi$ | $-0.48$ | $0.$
$\lambda$ | $-1.48$ | $1.5$
$\beta$ | $-0.53$ | $0.2$
$\psi$ | $1.8$ | $1.66$
SNR no SGWB | $228.8$ | $5096.2$
SNR power law | $29.4$ | $239.6$
SNR bump low-f | $68.7$ | $1411.4$
SNR bump high-f | $84.9$ | $442.1$
Table 1:  Parameters of the MBHB systems chosen for a simulated Bayesian
parameter estimation.

###  III.2 In-depth discussion for illustrative injections

To complement the analysis of the previous section, we delve into the detailed
parameter degeneracies for a few illustrative events. We employ a direct
sampling approach of the likelihood distribution (4) to capture potential non-
Gaussian features and multimodalities of the posterior distributions. The
parameters for our representative events are listed in Table 1. We select two
injections to represent typical events predicted by our populations. While
MBHB 1 is more likely found in populations similar to PopIII , it can also
appear as an outlier in populations like Q3delays and Q3nodelays , and vice
versa for MBHB 2.

MBHB 1 and MBHB 2 are relatively close events, with $z=3$. Thus, they
exemplify strong signals that LISA would detect with a quite large SNR in the
absence of a large SGWB. The presence of an SGWB would still allow the
detection of these two events, albeit with a significantly reduced SNR
compared to the case without an SGWB. This reduction is most prominent for the
PL case, for which the SNR is reduced by a factor $\sim 8$ and $\sim 21$ for
MBHB 1 and MBHB 2 respectively. We illustrate the SGWB impact on the accuracy
of parameter constraints in Fig. 4, where we present the 1D and 2D posterior
distributions of some representative parameters.

Let us start by discussing the intrinsic parameters. We focus on the spin
projections along the angular momentum of the binary, $\chi_{-}$, $\chi_{+}$,
on the chirp mass $\mathcal{M}_{c}$, and the mass ratio $q$. Visual inspection
of the posterior distributions reveals Gaussian-like shapes for most
parameters, with slight non-Gaussianities observed in the contours involving
the spins. Interestingly, the position of the SGWB bump does not significantly
impact the estimation of the spins for both MBHB 1 and 2. On the other hand,
compared to Bump High $f$, the Bump Low $f$ scenario severely affects the
estimate of the chirp mass for MBHB 1. This can be attributed to the noise
level being comparable to the strain of signals like MBHB 1 during their
inspiral stage, thereby reducing the ability to constrain the chirp mass of
such events. Among the benchmark SGWBs considered, Bump Low $f$ provides the
worst-case scenario for constraining the chirp mass of events with relatively
small total masses. This situation may deteriorate further for events with
masses of the order of $10^{4}M_{\odot}$, where much of the information about
the inspiral phase could be obscured by a large SGWB peaked at small
frequencies. Except for this parameter, the PL case stands out as the most
pessimistic, providing the least favorable constraints, particularly evident
for MBHB 2, where constraints worsen by nearly an order of magnitude compared
to the no SGWB scenario, diminishing the expected accuracy for such bright
signals. Finally, we stress a non-trivial effect occurring in the $q-\chi_{+}$
plane in the MBHB 2 case. We notice that in the High $f$ Bump case, the
ellipse is oriented in a different direction with respect to the PL and Low
$f$ bump cases. This effect originates from the partial masking of the last
part of the GW signal, which effectively modifies the relative weights of the
contributions brought by the different frequencies to the total information.
Since most information on the spin parameters comes from the very last part of
the GW signal, the High $f$ Bump case degrades more substantially the
determination of the spins compared to all other parameters.

Moving on to the discussion of extrinsic parameters, we focus on the
luminosity distance to the source, $d_{L}$, the ecliptic latitude and
longitude, denoted with $\beta_{L}$, and $\lambda_{L}$, respectively. We start
by commenting on the degradation of $d_{L}$, which in both cases is more
significant in the PL case compared to the Low and High $f$ bumps. Notice the
nontrivial scaling of the degradation of the uncertainty in the determination
of $d_{L}$ for the two bumps in the MBHB 1 and MBHB 2 cases. This behavior
originates from the relative difference between the frequency of the mergers
and the bumps: for MBHB 1, the two bumps only mask part of the inspiral, while
for MBHB 2, the High $f$ bump masks part of the merger. The posterior
distributions of $d_{L}$ shown in Fig. 4 and the SNR values reported in Tab. 1
are thus consistent with simple scaling arguments, which suggest that the
uncertainty on $d_{L}$ scales inversely with the SNR.

In addition to the constraints on the distance of the events, we also notice a
degradation in the sky localization of the events. Information about the sky
localization comes from both the inspiral and merger stages of the signal. The
inspiral stage contributes more significantly for low-mass systems, while the
merger stage is more dominant for high-mass systems. In the presence of an
SGWB affecting either low or high frequencies, lighter or heavier systems will
experience a degradation in their localization information. This degradation
might also give rise to multimodality in the sky localization of the source,
making it challenging to precisely determine the true signal parameters even
for nearby sources. This effect is manifest in the panels of the two triangle
plots showing the constraints on the latitude $\beta_{L}$. Particularly for
MBHB 1 (and also in the PL case for MBHB 2), the injected event could not be
distinguished from a similar signal with a flipped sign of $\beta_{L}$.

The importance of this effect varies depending on the properties of the
population of MBHBs. We explore its relevance for the different populations
and SGWBs considered in this work in Fig. 5. To quantify this effect, we
compute the difference $\Delta\ln\mathcal{L}_{\rm reflected}$ between the log-
likelihood evaluated in the true parameters and the sky position most likely
to develop a secondary mode, corresponding to the flipped sign of $\beta_{L}$.
The $x$/$y$-axes of the panels in Fig. 5 correspond to this quantity evaluated
in the “No SGWB” case and in the presence of the SGWB, respectively. Although
the mapping from this simplified metric to the probability support of this
possible secondary mode in the sky is not simple, we can consider that for a
value of $\Delta\ln\mathcal{L}_{\rm reflected}$ smaller than $\sim 10$,
highlighted by a vertical/horizontal line, the secondary mode becomes
compatible with the injected values, meaning it cannot be excluded by the
measurement LISA will perform. We can see that, compatibly with the results
shown in Fig. 3, the number of sources featuring a multimodality in the sky
position depends on the population/SGWB benchmark, with the High $f$ bump and
Q3delays being consistently the less dramatic cases.

Figure 5: Scatter plot for our 3 populations showing the value of an
approximate measure for sky multimodality, $\Delta\ln\mathcal{L}_{\rm
reflected}$, computed in the absence of an SGWB (x-axes) versus its value
computed in the presence of our 3 benchmark SGWBs (y-axes). Events on the left
of the vertical dashed gray line are likely to feature multimodal sky
localizations even in the absence of an SGWB. Events in the bottom right
quadrants, i.e. to the right of the vertical and below the horizontal dashed
gray lines, are likely to feature a multimodal sky localization only when the
SGWB is included.

##  IV Other sources in the LISA band

Having carefully re-analyzed the science case for MBHBs, we now move on to the
other main astrophysical sources in the LISA band: SOBHBs and GBs. In the next
subsections, we will limit the analysis to a simple SNR computation. As we
will shortly see, this would be sufficient to draw important conclusions,
without employing more complex parameter estimation techniques like the one we
performed for MBHBs.

###  IV.1 Stellar Origin Black Hole Binaries

Figure 6: SNR histogram for a SOBHB system similar to GW150914, varying
inclination, phase, sky position and polarization, for the different SGWB
cases.

Figure 7: [Left] SNR histogram for Galactic Binaries from the LDC catalogue,
for the different SGWB cases. We start from a reduced catalogue with purely
instrumental $\mathrm{SNR}>5$, plotted in black. [Center] Percentage loss of
detached and interacting binaries with ${\rm SNR}>5$ (with respect to the No
SGWB case) as a function of frequency. [Right] Degradation of the SNR of the
verification binaries. We only plot those that have ${\rm SNR}>5$ in the No
SGWB case. In the left and right panels, dashed and solid lines indicate SNR
$=5$ and $8$ respectively.

As an illustrative case, we consider an event with characteristics similar to
GW150914 [63] and discuss the SGWB’s impact on its detectability prospects
with LISA. In practice, we choose the masses, $d_{L}$, and spins to match the
measured values, pick the initial frequency, say $f_{\rm i}$, such that the
signal exits the LISA band at the end after 4 years of data taking, and
randomize all the angular variables, i.e., inclination, sky position, and
phase. The results of this analysis are shown in Fig. 6. Notice that the peak
at low values in the SNR distribution, which appears in all cases, is mostly
due to the dependency on $\iota$. Indeed, the SNR is minimal for
$\iota=\pi/2$, which corresponds to the peak of the distribution for the
number of events.

Since with our choice for the initial frequency, the signal only appears in
the high-frequency part of the LISA band (for reference, see Fig. 1), the SNR
would largely be unaffected by the presence of a bump at low frequencies. For
the same reason, the high-frequency bump here constitutes the most pessimistic
case, which might reduce the SNR of a factor $\sim 3$. The PL case sits
between the two other cases with a degradation in the SNR of a factor $\sim
2$. In both cases, the degradation is rather significant and might completely
jeopardize the possibility of resolving individual SOBHBs with LISA.

In our analysis, we have assumed the initial frequency $f_{\rm i}$ for the
signal such that the system leaves the LISA sensitivity band after 4 years of
data taking. In principle, events could have different values of $f_{\rm i}$,
which would impact the results of our analysis. Increasing $f_{\rm i}$, the
system would spend less time in the LISA band, which reducing the SNR, but
without modifying the qualitative behaviors for the different SGWB benchmarks.
On the other hand, reducing $f_{\rm i}$ implies that the system does not leave
the LISA frequency band after 4 years of data taking, leading to more
interesting nontrivial considerations. In this case, the system only spans a
range of frequencies between $f_{\rm i}$ and a maximal frequency777Notice that
$f_{\rm e}$ decreases quite rapidly with $f_{\rm i}$. In particular, systems
with very low $f_{\rm i}$ are still deep in the inspiral regime. Similar to
GBs, these systems do not evolve much in frequency during the mission lifetime
and can be effectively treated as monochromatic signals., say $f_{\rm e}$. In
this case, depending on the specific choice of $f_{\rm i}$, the bump at low
frequency (or at some intermediate position) might impact the SNR more than
the High $f$ bump case. Similarly, assuming we stick to the choice of
determining $f_{\rm i}$ so that the signal leaves the LISA band after 4 years
of data taking, increasing the mass of the system would lower the merger
frequency and speed up the frequency evolution, shifting $f_{\rm i}$ to a
lower frequency. Even in this case, we expect some qualitative modifications
of the conclusions drawn in this section.

###  IV.2 Galactic Binaries

As a last example, we discuss the impact of large SGWBs on GB detections. For
this purpose, we use the LISA Data Challenge (LDC) catalog for GBs based on
the populations in [71, 72, 73]. As for SOBHBs, we restrict our analysis to
SNR considerations. Since the SNR of a specific source can only decrease in
the presence of an SGWB, we start by computing the SNR for all the events in
the catalog in the case where only the instrumental noise is present (i.e., we
neglect the confusion noise) and drop all the events with SNR $<5$, being
undetectable even in the most optimistic scenario. We then compute the SNR for
the reduced catalog for the “No SGWB” configuration and the three benchmarks
discussed throughout this paper. The results of the analysis are shown in the
left panel of Fig. 7, where we show variations of the SNR distribution and in
its center panel, where we show the percentage loss of GBs with SNR$>5$ as a
function of frequency, taking as a reference the No SGWB case. Since most GB
signals sit on the low end of the LISA frequency band (see Fig. 1), a bump at
low frequency has the greatest impact on the SNR distribution of the sources,
whereas a bump at high frequency has the least. Nevertheless, even in the
latter case nearly 90% of the GBs are lost at frequencies in the range
$\sim[3,\,15]$ mHz. The PL case, while not as severe as the low-frequency
bump, would still significantly reduce the number of sources LISA can resolve,
decreasing from tens of thousands to hundreds. As for SOBHBs, if the
degradation due to a large SGWB is present, it could dramatically threaten the
prospects of resolving individual GBs with LISA.

One major impact of this degradation would be on the detection of verification
binaries [74] plotted in Fig. 1. These are ultra-compact GBs with orbital
periods shorter than a few hours, usually composed of degenerate stellar
remnants such as white dwarfs, neutron stars, or SOBHs. Having already been
detected electromagnetically and expected to be individually detected due to
their strong GW signals, they are considered guaranteed sources for LISA that
are crucial for facilitating the functional tests of the instrument and
maximizing LISA’s scientific output. However, only a few of them would be
detectable in the presence of a significant SGWB with substantial support at
frequencies $\lesssim 0.01$ Hz. For the 3 benchmark SGWBs considered in this
work, we show explicitly the lost verification binaries in the right panel of
Fig. 7.

As a final comment, we stress that the strong reduction in the SNR of each
source would complicate the subtraction of signals from GBs from the data. As
a consequence, we expect a (possibly) significant impact on the overall
confusion noise due to these sources, which peaks at around $1$ mHz, see the
small bump in the blue sensitivity curve in Fig. 1. Qualitatively, due to a
less efficient removal of individual sources, we expect the signal to extend
to higher frequencies and to have a larger amplitude. An updated estimate of
the confusion noise would require techniques similar to the ones employed in
[20], which lay beyond the scopes of the current work.

##  V Discussion and conclusions

Gravitational wave (GW) detectors are sensitive to two types of signals:
deterministic sources, such as coalescing compact objects, and stochastic
gravitational wave backgrounds (SGWBs), produced by the incoherent
superposition of unresolvable sources or early Universe mechanisms. The latter
appear as additional noise in GW detectors, potentially dominating the
instrumental noise at certain frequencies and making deterministic sources
harder to resolve. Well-motivated models of the early Universe predict an SGWB
that could be significantly larger than LISA’s instrumental noise level. The
only existing constraints on the SGWB amplitude in the LISA frequency band are
experimental, coming from Cosmic Microwave Background (CMB) and Big Bang
Nucleosynthesis (BBN).

Our paper is based on the observation that current bounds on the SGWB still
allow for amplitudes large enough to threaten the detection of deterministic
sources. We conducted a comprehensive study on how a strong SGWB, saturating
the CMB+BBN bounds, could degrade LISA’s capability to detect deterministic
sources. Our results should be considered as the most pessimistic scenario for
the LISA mission. Specifically, we analyzed the impact on three major
astrophysical targets of LISA: Massive Black Hole Binaries (MBHBs), Stellar
Origin Black Hole Binaries (SOBHBs), and Galactic Binaries (GBs). We modeled
the SGWB energy density, $h^{2}\Omega_{\rm GW}(f)$, as either a narrow bump
peaking at low frequency ($2.2\times 10^{-3}$ Hz) or high frequency
($2.2\times 10^{-2}$ Hz), or a flat power law across LISA’s sensitivity band.
While this approach is phenomenological, it reflects realistic scenarios:
peaked SGWBs can be produced by various inflationary mechanisms [31] or First
Order Phase Transitions [32], and nearly flat SGWBs are typically generated by
cosmic strings [33], which are expected to fall in the LISA band if the recent
PTA detection of a SGWB is interpreted in this context [75, 76, 77, 78].

Our results are summarized as follows.

  * •

MBHBs. We conducted a comprehensive study on how the distribution of the SNR
and the errors on the binary parameters, estimated via a Fisher analysis, are
impacted by an SGWB. We examined three representative types of MBHB
populations, each exemplifying different features. We further complemented
this with a detailed Bayesian analysis of a pair of representative events,
which we identified among the three populations. Our main findings are that a
strong SGWB will reduce the number of detectable events, predominantly
affecting populations with relatively small masses ($M/M_{\odot}\sim
10^{4}-10^{5}$). The exact number of non-detectable events depends on the
spectral shape of the SGWB. We also find that the errors on many binary
parameters can increase by up to an order of magnitude. Interestingly, a large
SGWB can strongly impact the sky localization, as we have observed through the
appearance of bimodal distributions in the event longitude, even for
relatively close events that would be well-localized without the SGWB.

  * •

SOBHBs. For these sources, we only performed an SNR analysis, which was
sufficient for our purposes. These types of signals are only affected by SGWBs
with support at high frequencies. Given that these sources are already
relatively faint, we found that the effect of the background is either to
entirely prevent their detection or to have nearly no effect on them.

  * •

Galactic binaries. For this type of source, we found results similar to those
for SOBHBs. In this case, the effect of the SGWB is also of the ”on-off” kind.
However, as opposed to SOBHBs, these sources emit at lower frequencies, and
SGWBs only marginally affecting the detections are those peaking at high
frequencies.

In short, our results show that there would be a strong impact on LISA’s
scientific objectives in the pessimistic scenario (for detecting individual GW
signals with LISA) where a primordial mechanism produces a loud SGWB. We find
that the study of MBHBs would be significantly affected, although still
possible to a more limited extent. On the other hand, SOBHBs and GBs could be
rendered completely undetectable, potentially nullifying two of LISA’s major
objectives. While a detailed study of the dependence of our assessments on the
properties of the SGWB is beyond the scope of this work, we point out that
these considerations are often not taken into account when assessing the
mission’s capabilities. For example, it would be interesting to quantify, the
maximal SGWB amplitude leading to the minimal degradation in detecting the
sources mentioned above. Beyond the sources we have studied, these important
results warrant further investigation into the impact of the SGWB on the
detection and characterization of other types of signals, such as Extreme Mass
Ratio Inspirals or other burst-like signals.

To conclude, let us discuss the implications of our study on some of the major
astrophysical and cosmological tests that the community hopes to perform by
detecting the sources we studied. Beginning with MBHBs, the reduction in the
number of detections and the degradation of the forecast errors will impact
the astrophysical and cosmological tests which require a large statistical
sample of events. These studies include population analyses, which provide
insights into the formation and evolution of MBHs, and cosmological studies
aimed at testing the expansion rate of the Universe and potential
modifications to gravity. However, we expect that tests of the nature of BHs
in the strong gravity regime or those relying on detecting electromagnetic
counterparts could still be performed by detecting some very strong sources,
which would be less affected by the SGWB.

Moving on to SOBHBs, a loud SGWB may completely ruin the prospects of
detecting this class of sources with LISA. Thus, the most direct implication
would be the loss of multiband GW detection [79]. A secondary yet significant
consequence would be the partial or even total masking of the SGWB generated
by these objects. As demonstrated in [22], measuring such a signal could
provide substantial insights into the behavior of these objects at high
redshift, revealing properties of the population that complement those tested
by LVK. A large SGWB would make it challenging to achieve this intriguing
prospect.

Finally, our results on GBs suggest that a significant SGWB could potentially
disrupt the program of verification binaries, which are expected to play a
crucial role in the science verification of LISA data.

The findings of our paper highlight the crucial need to consider the potential
impact of an SGWB on the detection of deterministic sources. Given the vast
range of theoretical predictions, uncertainties related to the SGWB might
outweigh those in the characterization of LISA’s noise. We hope this study
spurs further research to devise strategies that minimize the impact of an
SGWB and fully exploit LISA’s capabilities. On an optimistic note, even if the
scientific goals related to the detection of deterministic sources would be
somewhat compromised, detecting an SGWB of the amplitude considered in this
paper would be groundbreaking. It would confirm its primordial origin and mark
a significant milestone in physics, showcasing the remarkable potential of the
LISA instrument.

###### Acknowledgements.

We thank Nicola Tamanini for interesting discussions in the early stages of
this project. We also thank Valeriya Korol and Luca Reali for very useful
comments on a draft of this paper. MP acknowledges the hospitality of Imperial
College London, which provided office space during some parts of this project.
MP also thanks the Center for Cosmology and Particle Physics (CCPP) at New
York University for its hospitality during the initial and final parts of the
project.

## References

  * [1] P. Amaro-Seoane et al. [LISA], [arXiv:1702.00786 [astro-ph.IM]].
  * [2] P. A. Seoane et al. [LISA], Living Rev. Rel. 26, no.1, 2 (2023) doi:10.1007/s41114-022-00041-y [arXiv:2203.06016 [gr-qc]].
  * [3] K. G. Arun et al. [LISA], Living Rev. Rel. 25, no.1, 4 (2022) doi:10.1007/s41114-022-00036-9 [arXiv:2205.01597 [gr-qc]].
  * [4] L. Barack and C. Cutler, Phys. Rev. D 75, 042003 (2007) doi:10.1103/PhysRevD.75.042003 [arXiv:gr-qc/0612029 [gr-qc]].
  * [5] E. Barausse, E. Berti, T. Hertog, S. A. Hughes, P. Jetzer, P. Pani, T. P. Sotiriou, N. Tamanini, H. Witek and K. Yagi, et al. Gen. Rel. Grav. 52, no.8, 81 (2020) doi:10.1007/s10714-020-02691-1 [arXiv:2001.09793 [gr-qc]].
  * [6] E. Berti, E. Barausse, V. Cardoso, L. Gualtieri, P. Pani, U. Sperhake, L. C. Stein, N. Wex, K. Yagi and T. Baker, et al. Class. Quant. Grav. 32, 243001 (2015) doi:10.1088/0264-9381/32/24/243001 [arXiv:1501.07274 [gr-qc]].
  * [7] T. Bogdanovic, M. C. Miller and L. Blecha, Living Rev. Rel. 25, no.1, 3 (2022) doi:10.1007/s41114-022-00037-8 [arXiv:2109.03262 [astro-ph.HE]].
  * [8] K. Breivik, K. Kremer, M. Bueno, S. L. Larson, S. Coughlin and V. Kalogera, Astrophys. J. Lett. 854, no.1, L1 (2018) doi:10.3847/2041-8213/aaaa23 [arXiv:1710.08370 [astro-ph.SR]].
  * [9] A. Lamberts, S. Garrison-Kimmel, P. Hopkins, E. Quataert, J. Bullock, C. A. Faucher-Giguère, A. Wetzel, D. Keres, K. Drango and R. Sanderson, Mon. Not. Roy. Astron. Soc. 480, no.2, 2704-2718 (2018) doi:10.1093/mnras/sty2035 [arXiv:1801.03099 [astro-ph.GA]].
  * [10] G. Nelemans, L. R. Yungelson and S. F. Portegies Zwart, Astron. Astrophys. 375, 890-898 (2001) doi:10.1051/0004-6361:20010683 [arXiv:astro-ph/0105221 [astro-ph]].
  * [11] A. J. Ruiter, K. Belczynski, M. Benacquista, S. L. Larson and G. Williams, Astrophys. J. 717, 1006-1021 (2010) doi:10.1088/0004-637X/717/2/1006 [arXiv:0705.3272 [astro-ph]].
  * [12] C. Caprini and N. Tamanini, JCAP 10, 006 (2016) doi:10.1088/1475-7516/2016/10/006 [arXiv:1607.08755 [astro-ph.CO]].
  * [13] R. G. Cai, N. Tamanini and T. Yang, JCAP 05, 031 (2017) doi:10.1088/1475-7516/2017/05/031 [arXiv:1703.07323 [astro-ph.CO]].
  * [14] M. Corman, A. Ghosh, C. Escamilla-Rivera, M. A. Hendry, S. Marsat and N. Tamanini, Phys. Rev. D 105, no.6, 064061 (2022) doi:10.1103/PhysRevD.105.064061 [arXiv:2109.08748 [gr-qc]].
  * [15] M. Corman, C. Escamilla-Rivera and M. A. Hendry, JCAP 02, 005 (2021) doi:10.1088/1475-7516/2021/02/005 [arXiv:2004.04009 [astro-ph.CO]].
  * [16] E. Belgacem et al. [LISA Cosmology Working Group], JCAP 07, 024 (2019) doi:10.1088/1475-7516/2019/07/024 [arXiv:1906.01593 [astro-ph.CO]].
  * [17] P. A. Seoane, M. A. Sedda, S. Babak, C. P. L. Berry, E. Berti, G. Bertone, D. Blas, T. Bogdanović, M. Bonetti and K. Breivik, et al. Gen. Rel. Grav. 54, no.1, 3 (2022) doi:10.1007/s10714-021-02889-x [arXiv:2107.09665 [astro-ph.IM]].
  * [18] M. Colpi, K. Danzmann, M. Hewitson, P. Jetzer, G. Nelemans, A. Petiteau, D. Shoemaker, C. Sopuerta, R. Stebbins and N. Tanvir, et al. [arXiv:2402.07571 [astro-ph.CO]].
  * [19] S. Nissanke, M. Vallisneri, G. Nelemans and T. A. Prince, Astrophys. J. 758, 131 (2012) doi:10.1088/0004-637X/758/2/131 [arXiv:1201.4613 [astro-ph.GA]].
  * [20] N. Karnesis, S. Babak, M. Pieroni, N. Cornish and T. Littenberg, Phys. Rev. D 104, no.4, 043019 (2021) doi:10.1103/PhysRevD.104.043019 [arXiv:2103.14598 [astro-ph.IM]].
  * [21] G. Boileau, A. Lamberts, N. J. Cornish and R. Meyer, Mon. Not. Roy. Astron. Soc. 508, no.1, 803-826 (2021) [erratum: Mon. Not. Roy. Astron. Soc. 508, no.4, 5554-5555 (2021)] doi:10.1093/mnras/stab2575 [arXiv:2105.04283 [gr-qc]].
  * [22] S. Babak, C. Caprini, D. G. Figueroa, N. Karnesis, P. Marcoccia, G. Nardini, M. Pieroni, A. Ricciardone, A. Sesana and J. Torrado, JCAP 08, 034 (2023) doi:10.1088/1475-7516/2023/08/034 [arXiv:2304.06368 [astro-ph.CO]].
  * [23] M. Bonetti and A. Sesana, Phys. Rev. D 102 (2020) no.10, 103023 doi:10.1103/PhysRevD.102.103023 [arXiv:2007.14403 [astro-ph.GA]].
  * [24] F. Pozzoli, S. Babak, A. Sesana, M. Bonetti and N. Karnesis, Phys. Rev. D 108 (2023) no.10, 103039 doi:10.1103/PhysRevD.108.103039 [arXiv:2302.07043 [astro-ph.GA]].
  * [25] S. Naoz and Z. Haiman, Astrophys. J. Lett. 955 (2023) no.2, L27 doi:10.3847/2041-8213/acf8c9 [arXiv:2307.11149 [astro-ph.HE]].
  * [26] N. Bartolo, C. Caprini, V. Domcke, D. G. Figueroa, J. Garcia-Bellido, M. C. Guzzetti, M. Liguori, S. Matarrese, M. Peloso and A. Petiteau, et al. JCAP 12, 026 (2016) doi:10.1088/1475-7516/2016/12/026 [arXiv:1610.06481 [astro-ph.CO]].
  * [27] C. Caprini, M. Chala, G. C. Dorsch, M. Hindmarsh, S. J. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No and K. Rummukainen, et al. JCAP 03, 024 (2020) doi:10.1088/1475-7516/2020/03/024 [arXiv:1910.13125 [astro-ph.CO]].
  * [28] P. Auclair, J. J. Blanco-Pillado, D. G. Figueroa, A. C. Jenkins, M. Lewicki, M. Sakellariadou, S. Sanidas, L. Sousa, D. A. Steer and J. M. Wachter, et al. JCAP 04, 034 (2020) doi:10.1088/1475-7516/2020/04/034 [arXiv:1909.00819 [astro-ph.CO]].
  * [29] C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, no.16, 163001 (2018) doi:10.1088/1361-6382/aac608 [arXiv:1801.04268 [astro-ph.CO]].
  * [30] P. Auclair et al. [LISA Cosmology Working Group], Living Rev. Rel. 26, no.1, 5 (2023) doi:10.1007/s41114-023-00045-2 [arXiv:2204.05434 [astro-ph.CO]].
  * [31] [LISA Cosmology Working Group], To appear.
  * [32] C. Caprini et al. [LISA Cosmology Working Group], [arXiv:2403.03723 [astro-ph.CO]].
  * [33] J. J. Blanco-Pillado et al. [LISA Cosmology Working Group], [arXiv:2405.03740 [astro-ph.CO]].
  * [34] L. Pagano, L. Salvati and A. Melchiorri, Phys. Lett. B 760, 823-825 (2016) doi:10.1016/j.physletb.2016.07.078 [arXiv:1508.02393 [astro-ph.CO]].
  * [35] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L8 (2023) doi:10.3847/2041-8213/acdac6 [arXiv:2306.16213 [astro-ph.HE]].
  * [36] G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L9 (2023) doi:10.3847/2041-8213/acda9a [arXiv:2306.16217 [astro-ph.HE]].
  * [37] J. Antoniadis et al. [EPTA and InPTA:], Astron. Astrophys. 678, A50 (2023) doi:10.1051/0004-6361/202346844 [arXiv:2306.16214 [astro-ph.HE]].
  * [38] J. Antoniadis et al. [EPTA], Astron. Astrophys. 678, A48 (2023) doi:10.1051/0004-6361/202346841 [arXiv:2306.16224 [astro-ph.HE]].
  * [39] J. Antoniadis et al. [EPTA], [arXiv:2306.16227 [astro-ph.CO]].
  * [40] D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane and J. Askew, et al. Astrophys. J. Lett. 951, no.1, L6 (2023) doi:10.3847/2041-8213/acdd02 [arXiv:2306.16215 [astro-ph.HE]].
  * [41] A. Zic, D. J. Reardon, A. Kapur, G. Hobbs, R. Mandow, M. Curyło, R. M. Shannon, J. Askew, M. Bailes and N. D. R. Bhat, et al. Publ. Astron. Soc. Austral. 40, e049 (2023) doi:10.1017/pasa.2023.36 [arXiv:2306.16230 [astro-ph.HE]].
  * [42] D. J. Reardon, A. Zic, R. M. Shannon, V. Di Marco, G. B. Hobbs, A. Kapur, M. E. Lower, R. Mandow, H. Middleton and M. T. Miles, et al. Astrophys. J. Lett. 951, no.1, L7 (2023) doi:10.3847/2041-8213/acdd03 [arXiv:2306.16229 [astro-ph.HE]].
  * [43] H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. N. Caballero, J. Yuan and Y. Xu, et al. Res. Astron. Astrophys. 23, no.7, 075024 (2023) doi:10.1088/1674-4527/acdfa5 [arXiv:2306.16216 [astro-ph.HE]].
  * [44] M. R. Adams and N. J. Cornish, Phys. Rev. D 82, 022002 (2010) doi:10.1103/PhysRevD.82.022002 [arXiv:1002.1291 [gr-qc]].
  * [45] T. A. Prince, M. Tinto, S. L. Larson and J. W. Armstrong, Phys. Rev. D 66, 122002 (2002) doi:10.1103/PhysRevD.66.122002 [arXiv:gr-qc/0209039 [gr-qc]].
  * [46] S. Babak, A. Petiteau and M. Hewitson, [arXiv:2108.01167 [astro-ph.IM]].
  * [47] C. Pitte, Q. Baghi, S. Marsat, M. Besançon and A. Petiteau, Phys. Rev. D 108 (2023) no.4, 044053 doi:10.1103/PhysRevD.108.044053 [arXiv:2304.03142 [gr-qc]].
  * [48] L. Speri, A. Antonelli, L. Sberna, S. Babak, E. Barausse, J. R. Gair and M. L. Katz, Phys. Rev. X 13 (2023) no.2, 021035 doi:10.1103/PhysRevX.13.021035 [arXiv:2207.10086 [gr-qc]].
  * [49] T. L. Smith and R. Caldwell, Phys. Rev. D 95, no.4, 044036 (2017) doi:10.1103/PhysRevD.95.044036 [arXiv:1609.05901 [gr-qc]].
  * [50] T. L. Smith, T. L. Smith, R. R. Caldwell and R. Caldwell, Phys. Rev. D 100, no.10, 104055 (2019) [erratum: Phys. Rev. D 105, no.2, 029902 (2022)] doi:10.1103/PhysRevD.100.104055 [arXiv:1908.00546 [astro-ph.CO]].
  * [51] R. Flauger, N. Karnesis, G. Nardini, M. Pieroni, A. Ricciardone and J. Torrado, JCAP 01, 059 (2021) doi:10.1088/1475-7516/2021/01/059 [arXiv:2009.11845 [astro-ph.CO]].
  * [52] C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena and R. Jaume, Phys. Rev. D 102 (2020) no.6, 064002 doi:10.1103/PhysRevD.102.064002 [arXiv:2001.10914 [gr-qc]].
  * [53] C. García-Quirós, S. Husa, M. Mateu-Lucena and A. Borchers, Class. Quant. Grav. 38 (2021) no.1, 015006 doi:10.1088/1361-6382/abc36e [arXiv:2001.10897 [gr-qc]].
  * [54] E. E. Flanagan and S. A. Hughes, New J. Phys. 7, 204 (2005) doi:10.1088/1367-2630/7/1/204 [arXiv:gr-qc/0501041 [gr-qc]].
  * [55] S. Marsat, J. G. Baker and T. Dal Canton, Phys. Rev. D 103, no.8, 083011 (2021) doi:10.1103/PhysRevD.103.083011 [arXiv:2003.00357 [gr-qc]].
  * [56] W. D. Vousden, W. M. Farr and I. Mandel, Mon. Not. Roy. Astron. Soc. 455, no.2, 1919-1937 (2016) doi:10.1093/mnras/stv2422
  * [57] J. Goodman, J. Weare, Comm. in App. Math. and Comp. Sci. 5 No. 1, 65-80 (2010) doi:10.2140/camcos.2010.5.65.
  * [58] D. Foreman-Mackey, D. W. Hogg, D. Lang and J. Goodman, Publ. Astron. Soc. Pac. 125, 306-312 (2013) doi:10.1086/670067 [arXiv:1202.3665 [astro-ph.IM]].
  * [59] M. Vallisneri, Phys. Rev. D 77, 042001 (2008) doi:10.1103/PhysRevD.77.042001 [arXiv:gr-qc/0703086 [gr-qc]].
  * [60] L. Reali, A. Antonelli, R. Cotesta, S. Borhanian, M. Çalışkan, E. Berti and B. S. Sathyaprakash, [arXiv:2209.13452 [gr-qc]].
  * [61] L. Reali, A. Maselli and E. Berti, [arXiv:2307.01264 [gr-qc]].
  * [62] N. Bartolo et al. [LISA Cosmology Working Group], JCAP 11, 009 (2022) doi:10.1088/1475-7516/2022/11/009 [arXiv:2201.08782 [astro-ph.CO]].
  * [63] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116, no.6, 061102 (2016) doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * [64] T. Kupfer, V. Korol, T. B. Littenberg, S. Shah, E. Savalle, P. J. Groot, T. R. Marsh, M. Le Jeune, G. Nelemans and A. F. Pala, et al. Astrophys. J. 963 (2024) no.2, 100 doi:10.3847/1538-4357/ad2068 [arXiv:2302.12719 [astro-ph.SR]].
  * [65] https://gitlab.in2p3.fr/LISA/lisa-verification-binaries
  * [66] A. Klein, E. Barausse, A. Sesana, A. Petiteau, E. Berti, S. Babak, J. Gair, S. Aoudia, I. Hinder and F. Ohme, et al. Phys. Rev. D 93, no.2, 024003 (2016) doi:10.1103/PhysRevD.93.024003 [arXiv:1511.05581 [gr-qc]].
  * [67] E. Barausse, Mon. Not. Roy. Astron. Soc. 423, 2533-2557 (2012) doi:10.1111/j.1365-2966.2012.21057.x [arXiv:1201.5888 [astro-ph.CO]].
  * [68] E. Barausse and A. Lapi, doi:10.1007/978-981-15-4702-7_18-1 [arXiv:2011.01994 [astro-ph.GA]].
  * [69] E. Barausse, I. Dvorkin, M. Tremmel, M. Volonteri and M. Bonetti, Astrophys. J. 904, no.1, 16 (2020) doi:10.3847/1538-4357/abba7f [arXiv:2006.03065 [astro-ph.GA]].
  * [70] M. Volonteri, G. Lodato and P. Natarajan, Mon. Not. Roy. Astron. Soc. 383, 1079 (2008) doi:10.1111/j.1365-2966.2007.12589.x [arXiv:0709.0529 [astro-ph]].
  * [71] V. Korol, E. M. Rossi, P. J. Groot, G. Nelemans, S. Toonen and A. G. A. Brown, Mon. Not. Roy. Astron. Soc. 470 (2017) no.2, 1894-1910 doi:10.1093/mnras/stx1285 [arXiv:1703.02555 [astro-ph.HE]].
  * [72] V. Korol, S. Toonen, A. Klein, V. Belokurov, F. Vincenzo, R. Buscicchio, D. Gerosa, C. J. Moore, E. Roebber and E. M. Rossi, et al. Astron. Astrophys. 638, A153 (2020) doi:10.1051/0004-6361/202037764 [arXiv:2002.10462 [astro-ph.GA]].
  * [73] G. Nelemans, ASP Conf. Ser. 467, 27-36 (2013) [arXiv:1302.0138 [astro-ph.HE]].
  * [74] A. Stroeer and A. Vecchio, Class. Quant. Grav. 23, S809-S818 (2006) doi:10.1088/0264-9381/23/19/S19 [arXiv:astro-ph/0605227 [astro-ph]].
  * [75] J. Ellis, M. Lewicki, C. Lin and V. Vaskonen, Phys. Rev. D 108, no.10, 103511 (2023) doi:10.1103/PhysRevD.108.103511 [arXiv:2306.17147 [astro-ph.CO]].
  * [76] G. Lazarides, R. Maji and Q. Shafi, Phys. Rev. D 108, no.9, 095041 (2023) doi:10.1103/PhysRevD.108.095041 [arXiv:2306.17788 [hep-ph]].
  * [77] D. G. Figueroa, M. Pieroni, A. Ricciardone and P. Simakachorn, Phys. Rev. Lett. 132, no.17, 171002 (2024) doi:10.1103/PhysRevLett.132.171002 [arXiv:2307.02399 [astro-ph.CO]].
  * [78] J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen and H. Veermäe, Phys. Rev. D 109 (2024) no.2, 023522 doi:10.1103/PhysRevD.109.023522 [arXiv:2308.08546 [astro-ph.CO]].
  * [79] A. Sesana, Phys. Rev. Lett. 116, no.23, 231102 (2016) doi:10.1103/PhysRevLett.116.231102 [arXiv:1602.06951 [gr-qc]].