# Can gravitational-wave memory help constrain binary black-hole parameters?
A LISA case study

Silvia Gasparotto<EMAIL_ADDRESS>Institut de Fisica d’Altes Energies
(IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193
Bellaterra (Barcelona), Spain  Rodrigo Vicente  Institut de Fisica d’Altes
Energies (IFAE), The Barcelona Institute of Science and Technology, Campus
UAB, 08193 Bellaterra (Barcelona), Spain  Diego Blas  Institut de Fisica
d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology,
Campus UAB, 08193 Bellaterra (Barcelona), Spain  Grup de Física Teòrica,
Departament de Física, Universitat Autònoma de Barcelona, 08193 Bellaterra,
Spain  Alexander C. Jenkins  Department of Physics and Astronomy, University
College London, London WC1E 6BT, United Kingdom  Enrico Barausse  SISSA, Via
Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste  IFPU -
Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy

###### Abstract

Besides the transient effect, the passage of a gravitational wave also causes
a persistent displacement in the relative position of an interferometer’s test
masses through the _nonlinear memory effect_. This effect is generated by the
gravitational backreaction of the waves themselves, and encodes additional
information about the source. In this work, we explore the implications of
using this information for the parameter estimation of massive binary black
holes with LISA. Based on a Fisher analysis for nonprecessing black hole
binaries, our results show that the memory can help to reduce the degeneracy
between the luminosity distance and the inclination for binaries observed only
for a short time ($\sim$ few hours) before merger. To assess how many such
short signals will be detected, we utilized state-of-the-art predictions for
the population of massive black hole binaries and models for the gaps expected
in the LISA data. We forecast from tens to few hundreds of binaries with
observable memory, but only $\sim\mathcal{O}(0.1)$ events in 4 years for which
the memory helps to reduce the degeneracy between distance and inclination.
Based on this, we conclude that the new information from the nonlinear memory,
while promising for testing general relativity in the strong field regime, has
probably a limited impact on further constraining the uncertainty on massive
black hole binary parameters with LISA.

##  I Introduction

The direct detection of gravitational waves (GWs), predicted by Einstein in
1916 [1], is one of the greatest accomplishments in modern physics, showing
(at present) a spectacular agreement with the theory of general relativity
(GR) [2, 3, 4]. By now, almost a hundred GW signals have been observed and
interpreted as resulting from the coalescence of compact binaries by
LIGO/Virgo [5, 6, 7]. As the sensitivity of current detectors improves and new
detectors become available, it will be possible to estimate the binary
parameters more accurately and to find GWs from new types of sources. This
will allow us to not only better test GR in its strong-field regime, but also
to probe astrophysics, cosmology and fundamental physics [8, 9, 10]. The
future space-borne detector Laser Interferometer Space Antenna (LISA) [11]
will play a key role in this quest, due to both its expected high signal-to-
noise ratio (SNR) measurements and the rich population of sources expected to
inhabit its frequency band (from 0.1 mHz to 0.1 Hz) [12, 13, 14].

Progress may still be hindered by the fact that some binary parameters — such
as the luminosity distance $d_{\mathrm{L}}$ and inclination of the orbital
plane with respect to the line of sight $\iota$ — may be highly correlated in
GW signals, limiting our ability to accurately estimate them [15, 16, 17].
This is especially important in the context of _standard sirens_ [18, 19],
where the precision with which one can estimate the present-day Hubble
parameter $H_{0}$ depends primarily on how accurately one can measure
$d_{\mathrm{L}}$ [20, 21, 17]. Indeed, this was the main contribution to the
large uncertainty on the first estimate of $H_{0}$ from GW170817 in Ref. [22].
The distance-inclination degeneracy can be simply understood by noting that,
at leading (Newtonian) order, an inspiralling binary sources the two GW
polarizations $h_{+}\propto(1+\cos^{2}\iota)/d_{\mathrm{L}}$ and
$h_{\times}\propto\cos\iota/d_{\mathrm{L}}$ [23]; if the detector network is
mostly sensitive to one particular combination of $h_{+}$ and $h_{\times}$,
the luminosity distance and inclination are therefore degenerate [15, 17]
(c.f. Appendix B). The degree of this degeneracy depends on the sky location
of the binary and the specific detector network [17], and can be greatly
reduced by the observation of the afterglow light curve of an electromagnetic
counterpart (which critically depends on $\iota$) [24, 25]. Interestingly, the
degeneracy may also be mitigated by using subleading effects in the waveform.
Examples include the effect of spins misaligned with the orbital angular
momentum [26, 27, 28] (which lead to the precession of the orbital plane), of
higher multipole modes (HMs) [29, 30, 31, 32, 33, 34, 35, 36, 37] (in
particular, for unequal component masses or eccentric orbits), or using binary
Love relations [38] (for neutron star binaries).111Nevertheless, in the $\sim
100$ GW signals observed to date there is only limited evidence for higher
multipole content (with no evidence at all beyond $\ell=3$) [39] and only one
measurement of strong-field precession has been claimed [40] (though some
doubt has been cast on this claim due to data-quality issues [41, 42]).

Another subleading effect in the GW signal with the potential to break the
distance-inclination degeneracy is the nonlinear (Christodoulou) GW _memory_
[43]. This is a well-grounded prediction of GR which originates from a change
in the radiative multipole moments of the gravitational field sourced by the
flux of gravitational radiation itself, resulting in a permanent displacement
of free-falling test masses upon the passage of GWs [44, 43, 45, 46, 47, 48,
49].222A similar effect sourced by the flux of matter or non-gravitational
radiation was actually discovered first and is called the _linear_ GW memory
[50, 51, 52]. Hereafter, we use “memory” to refer exclusively to
Christodoulou’s GW memory. While essentially any source of GWs will generate
nonlinear memory,333See, e.g., Refs. [53, 54], which studied the nonlinear
memory generated by cosmic string loops. our focus here (and in much of the
relevant literature) is on binary black holes (BBHs). The reason for this is
twofold: first, binaries are the one source of _detectable_ GWs we
definitively know to exist; and second, the amplitude of the memory scales
with the total GW energy radiated, which for binaries is $\sim 1\%-10\%$ of
the total mass [55], favoring BBHs over lighter binaries containing neutron
stars or white dwarfs.

The gravitational wave memory modifies the BBH waveform by introducing a
slowly-growing offset of the oscillations that builds over the whole
coalescence and whose time evolution follows that of the instantaneous orbital
frequency. This shift rises over the radiation-reaction timescale during the
inspiral [47, 48, 56] and accumulates rapidly during the merger before
saturating to its final value during the ringdown [57]. Although the memory
arises from a 2.5 PN nonlinear interaction in a post-Newtonian (PN) expansion
of Einstein equations, because it accumulates over the whole coalescence, it
affects the gravitational waveform at leading (Newtonian) order [57],
increasing the prospects of observing it in the near future.

Several searches for memory from BBHs have been performed using LIGO/Virgo
data, returning only null results thus far [58, 59, 60, 61]. This is in
agreement with forecasts for LIGO/Virgo, which show that the detection of
memory from a single event would require a much more massive and nearby binary
than any yet observed, and that to find collective evidence of memory in the
total population of observed binaries one would need $\sim 5$ yr of collected
data [62, 63, 64] (or $\sim 2.5$ yr for the LIGO/Virgo/KAGRA network, taking
into account the expected improvement of the network sensitivities [65]). The
difficulty in detecting the memory with current ground-based interferometers
resides mostly in the fact that, besides being responsible for only a small
amplitude offset, its power is larger at lower frequencies where the detectors
sensitivity is limited by several sources of noise [66]. However, the
prospects for memory detection in single events are considerably better for
third-generation ground-based detectors (e.g. Einstein Telescope [67] and
Cosmic Explorer [68]) and for the future space-based detectors LISA and
TianQin, due to their better sensitivity and low frequency coverage [57, 69,
70, 71, 72, 65].444Memory from the merger of supermassive BHs is also a target
of Pulsar Timing Arrays (PTAs) [73, 74], but searches in PTA data thus far
have returned only null results [75, 76].

In this work we investigate the impact of the nonlinear memory on parameter
estimation via a Fisher matrix analysis. In particular, we focus on how the
memory signal breaks the distance-inclination degeneracy, which is crucially
important if these binaries are to be used as standard sirens. We found that
the information of the memory can indeed reduce the uncertainty on the
luminosity distance by reducing its correlation with the inclination angle,
whereas it has almost no impact on the uncertainty of the other parameters.
For LISA sources its greatest effect involves cases where (i) the constituent
BHs are light enough [$M_{\mathrm{BH}}\lesssim 10^{5}\,M_{\odot}/(1+z)$, with
$z$ the source’s cosmological redshift] that the merger takes place near the
upper edge of the LISA band, and (ii) the information from the primary
waveform is limited to a few cycles before the merger.

The presence of gaps in the data stream and confusion noise from other sources
will reduce the effective duration of usable LISA data, and thus the observed
number of cycles for BBHs [77], making the memory potentially useful for the
distance estimation of some BBH events. Therefore, using the state-of-art
astrophysical BBH population models described in Refs. [78, 79] (and based
upon previous work presented in [80, 81, 82]), we assess quantitatively the
impact of the memory in the distance estimation of LISA sources, taking into
account the presence of gaps in the data stream. Considering the particular
gap model used in Ref. [83] we did _not_ find any significant enhancement in
the distance estimation by the inclusion of memory on the BBH waveforms. We do
however find a greater number of events with detectable memory at LISA as
compared to previous forecasts [72, 70], especially in our models with heavy
BH seeds.

This paper is organized as follows. In Sec. II we review the computation of
the memory and describe the phenomenology of the signal. In Sec. III we
describe our Fisher forecasting analysis. In Sec. IV we present our results
for the distance-inclination inference of individual BBHs, and discuss the
impact of the binary parameters and signal duration. In Sec. V we perform
population-level forecasts for LISA using synthetic BBH catalogues, and assess
the impact of including the memory on the luminosity distance estimation in
the presence of gaps in the data stream. We conclude in Sec. VI. Some
technical material is discussed in the appendices. We use geometric units
throughout ($c=G=1$).

##  II GW memory waveform

###  II.1 Computation scheme: Thorne’s formula

The most direct way to compute the memory contribution to waveforms would be
to extract it directly from numerical relativity simulations. However, most
simulations to date have struggled to accurately capture this information for
a number of reasons (see e.g. [47]). Some exceptions are Ref. [84], where the
dominant memory mode $(\ell,m)=(2,0)$ was first resolved, and the recent work
of Ref. [85], which used a Cauchy-characteristic extraction (CCE) technique to
extract the waveform. Alternatively, the Bondi, van der Burg, Metzner, and
Sachs (BMS) _balance laws_ [86, 87] have recently been used to add the memory
to waveforms [88, 89] (see also [60, 60, 72]).

Instead, in this work we use a perturbative approach to evaluate the memory
[49, 47], which we now briefly review. A GW strain $h_{0}$ (which we call the
“primary” signal) sources an additional memory strain $\delta h$, which can be
expressed in the transverse-traceless (TT) gauge using Thorne’s formula [90],

$\delta
h_{ij}^{\mathrm{TT}}(u)=\frac{4}{R}\int_{-\infty}^{u}\mathrm{d}u^{\prime}\int_{R}\mathrm{d}\Omega\,\frac{\mathrm{d}^{2}E_{\mathrm{GW}}}{\mathrm{d}u^{\prime}\mathrm{d}\Omega}\left[\frac{n_{i}n_{j}}{1-n_{k}N^{k}}\right]^{\mathrm{TT}},$
(1)

where the angular integral is over the solid angle $\mathrm{d}\Omega$ of a
(large) sphere of radius $R$ surrounding the source, and $n^{i}$ and $N^{i}$
are the unit radial vectors pointing, respectively, to $\mathrm{d}\Omega$ and
the detector. The TT superscript represents a TT projection with respect to
the direction of the detector. The time integral is over the entire history of
the source up to retarded time $u$, which shows that the memory is a
hereditary effect. The GW energy flux carried by the primary GW is [47]

$\frac{\mathrm{d}^{2}E_{\mathrm{GW}}}{\mathrm{d}t\mathrm{d}\Omega}=\frac{R^{2}}{16\pi}(\dot{h}_{0,+}^{2}+\dot{h}_{0,\times}^{2}),$
(2)

where $\dot{h}\equiv\mathrm{d}h/\mathrm{d}t$ and $h_{+,\times}\equiv
h_{ij}^{\mathrm{TT}}e^{ij}_{+,\times}$. We use the same choice of TT-
polarization tensors $e^{ij}_{+,\times}$ as Ref. [91]. From the spin-weighted
spherical harmonic decomposition

$h_{+}-ih_{\times}\equiv\sum_{\ell\geq 2}\,\sum_{|m|\leq\ell}h^{\ell
m}(u,r)\,{}_{-2}Y_{\ell m}(\iota,\phi),$ (3)

it is possible to show that the sourced memory can be expressed as [59]

$\displaystyle\delta h^{\ell
m}(u)=-R\sum_{\ell^{\prime},\ell^{\prime\prime}\geq
2}\,\sum_{m^{\prime},m^{\prime\prime}}\sqrt{\frac{(\ell-2)!}{(\ell+2)!}}$
$\displaystyle\times\int\mathrm{d}\Omega\,Y_{\ell
m}^{*}\,{}_{-2}Y^{*}_{\ell^{\prime}m^{\prime}}\,{}_{-2}Y_{\ell^{\prime\prime}m^{\prime\prime}}\int_{-\infty}^{u}\mathrm{d}u^{\prime}\dot{h}_{0}^{*\ell^{\prime}m^{\prime}}\dot{h}_{0}^{\ell^{\prime\prime}m^{\prime\prime}},$
(4)

which allows us to straightforwardly compute the memory modes $\delta h^{\ell
m}$ from the primary waveform modes $h_{0}^{\ell m}$. We can then reconstruct
$\delta h_{+}$ and $\delta h_{\times}$ from Eq. (3) to give the total strain
$h\approx h_{0}+\delta h$. In principle, this process should be iterated to
give higher-order contributions (the “memory of the memory” [62]). In
practice, these extra terms are subleading, and it is sufficient for our
purposes to consider just the leading-order memory effect, $\delta h$.

To generate our primary waveforms, we use the surrogate NRHybSur3dq8 model
[92], which has only aligned-spin waveforms and includes all higher spherical
harmonic modes up to $(\ell,|m|)=(4,4)$ and is considerably more accurate than
other often-used phenomenological models in modeling the merger stage of BBH
coalescence [64]. We use the publicly available GWmemory package [93] to
implement the calculation scheme described above for the corresponding memory
signal.

Equation (II.1) is valid on a background Minkowski spacetime. It can be
extended to a spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW)
spacetime using the fact that, for sources at the same luminosity distance
$d_{\mathrm{L}}$, the memory amplitude in FLRW is enhanced over the Minkowski
case by the redshift factor $(1+z)$ [94, 95, 96]. Additionally, we shall use
the time at the detector $t\equiv
t_{\mathrm{peak}}-(1+z)(u_{\mathrm{peak}}-u)$, where $t_{\mathrm{peak}}$ is
the instant when the primary strain reaches its peak amplitude. Summarizing,
in this work we use

$\delta h^{\ell m}_{\mathrm{FLRW}}(t)=(1+z)\delta h^{\ell
m}_{\mathrm{Mink}}\big{(}u(t)\big{)}_{R\to d_{\mathrm{L}}}.$ (5)

This can be shown to be equivalent to using redshifted component masses
$M_{i,z}\equiv(1+z)M_{i}$, with $i\in\\{1,2\\}$, and luminosity distance
$d_{\mathrm{L}}$ to generate the primary signal (e.g., with NRHybSur3dq8),
plugging this primary directly into Eq. (II.1), and identifying
$(R,u_{\mathrm{peak}}-u)\to(d_{\mathrm{L}},t_{\mathrm{peak}}-t)$. We omit the
subscript “FLRW” throughout, but a spatially flat FLRW is implicitly assumed
in all our expressions.

###  II.2 Phenomenology of the signal

Equation (II.1) shows how the memory modes are sourced by pairs of primary
modes. For a BBH coalescence occurring in the $x$-$y$ plane, the primary modes
have the form $h_{0}^{\ell m}\propto e^{-im\varphi(t)}$ with $\varphi(t)$ the
orbital phase. So, from the time integral in Eq. (II.1) it is clear that the
leading contribution to the memory modes at low frequency (i.e.,
$u\to+\infty$) comes from the nonoscillatory (DC) terms with
$m^{\prime}-m^{\prime\prime}=0$ which accumulate in time [47, 48]; these
source $m=0$ memory modes (from the angular integral in Eq. (II.1)). Note,
however, that oscillatory $m\neq 0$ memory modes do become dominant at high
frequencies (cf. Fig. 1).

Figure 1: Memory characteristic strain $h_{\mathrm{c}}^{\ell m}(f)\equiv
2f|\widetilde{\delta h}{}^{\ell m}|$ of the most important modes computed from
Eq. (II.1). We consider the non-spinning “heavy” BBH studied in Figs. 2 and 3,
with total mass $M=2\times 10^{5}\,M_{\odot}$, mass ratio $q=1.2$ and redshift
$z=2$. The non-oscillatory $(\ell,m)=(2,0)$ mode dominates at low frequencies,
but is suppressed at $f\gtrsim 1/60M_{z}$, where oscillatory $m\neq 0$ modes
start becoming important (in particular, at their maxima $f\sim
mf_{\mathrm{peak}}$). We can also see the presence of the ringdown in the
memory spectrum in the form of high-frequency peaks [$\mathrm{QNM}_{1}$ is at
$f=(\tau_{221}^{-1}+\tau_{222}^{-1})/2\pi M_{z}$, and $\mathrm{QNM}_{2}$ at
$f=\tau_{222}^{-1}/\pi M_{z}$]. (QNMs were computed using the Python package
qnm [97], and the final mass and spin of the remnant BH via surfinBH [98],
whose fitting procedure is described in [99].)

The memory sourced in the quasi-circular inspiral of non-spinning BBHs is
known analytically up to 3 PN [47, 48]. Due to the accumulation over the
inspiral, the non-oscillatory contribution to the memory enters at Newtonian
(0 PN) order in the waveform,

$\displaystyle\delta h^{(0\mathrm{PN})}_{+}$ $\displaystyle=\frac{\eta
M_{z}}{48\,d_{\mathrm{L}}}[M\omega(t)]^{\frac{2}{3}}\sin^{2}\iota\,(17+\cos^{2}\iota),$
(6)

and $\delta h^{(0\mathrm{PN})}_{\times}=0$, with the conventional choice of
polarization triad [47]. The orbital frequency is
$\omega(t)\equiv\dot{\varphi}(t)$ and the symmetric mass ratio $\eta\equiv
M_{1}M_{2}/M^{2}$, where $M\equiv M_{1}+M_{2}$ is the total mass and
$M_{z}\equiv(1+z)M$ is the total redshifted mass. The memory has the same
scaling as the Newtonian primary waveform (cf. Eq. (13)), but rather than the
main time dependence coming from the oscillatory term, its time evolution is
captured by the instantaneous orbital frequency $\omega(t)$. This explains the
typical step shape of the memory, which has a steep increase in the merger-
plunge phase and a saturation during the ringdown. Moreover, the memory is
characterized by a _different_ dependence on the inclination angle $\iota$ and
an overall amplitude $\sim 20$ times weaker than the primary. In particular,
the two signals have opposite monotonic dependence on the inclination angle,
and while the primary signal is maximized for face-on binaries ($\iota=0$),
the memory is instead maximized for edge-on binaries ($\iota=\pi/2$). This
behavior is maintained when using primary waveforms generated by NRHybSur3dq8.
The different dependence on $\iota$ is what makes the memory helpful in
reducing the ($\iota,d_{\mathrm{L}}$) correlation (c.f. Fig. 4).

Figure 1 shows the spectral shape of memory mode characteristic strains
$h^{\ell m}_{\mathrm{c}}(f)\equiv 2f|\widetilde{\delta h}{}^{\ell m}(f)|$,
with $\widetilde{\delta h}(f)\equiv\int\mathrm{d}t\,e^{-i2\pi ft}\delta h(t)$
the Fourier transform (FT) of $\delta h$. All modes exhibit a plateau at low
frequencies, but the spectral content of the non-oscillatory ($m=0$) modes is
clearly distinct from the oscillatory ($m\neq 0$) modes. The plateau of the
$m=0$ modes is easily understood from the approximation $\delta h^{\ell
0}\approx\Delta h^{\ell 0}\mathcal{H}(t-t_{\mathrm{peak}})$, with
$\mathcal{H}$ the Heaviside step function and where $\Delta h^{\ell 0}$ scales
with the fraction of radiated $E_{\mathrm{GW}}$ that sources $\delta h^{\ell
0}$; this results in a constant $h_{\mathrm{c}}\approx\Delta h^{\ell 0}/\pi$
[56].

Taking into account that the memory growth is not instantaneous, but happens
in $\tau\sim 60M_{z}$ (the timescale over which most of $E_{\mathrm{GW}}$ is
radiated [100]), one can understand the suppression at $f\gtrsim 1/60M_{z}$ in
Fig. 1. For the $m\neq 0$ modes the low-frequency plateau has a similar
origin, but its value is always subdominant because, due to the oscillations
in the integral in Eq. (II.1), the memory does not accumulate, averaging out
to a net small value that depends strongly on the value of the orbital phase
at which the BHs merge; this is also expected from PN results [47, 48]. We
note that for $m\neq 0$ the maximum of the strain scales as
$mf_{\mathrm{peak}}$ where the peak frequency $f_{\mathrm{peak}}\sim 0.1/\pi
M_{z}$ roughly corresponds to the moment in which most of the energy is
radiated [100]. The oscillations at the left of these maxima are numerically
stable (in particular, they are not artifacts of our FT) and come from
interference between different $\ell$-modes in Eq. (II.1). On the other hand,
the high-frequency peaks seen in all modes are associated with the ringdown
stage and are located at
$f\approx(\tau_{\ell^{\prime}m^{\prime}n^{\prime}}^{-1}+\tau_{\ell^{\prime\prime}m^{\prime\prime}n^{\prime\prime}}^{-1})/2\pi
M_{z}$, where the complex quasi-normal modes are $\sigma_{\ell
mn}\equiv(\omega_{\ell mn}+i\tau_{\ell mn}^{-1})/M_{z}$ [101]. It is
interesting to note that these peaks occur at a frequency set by the QNM
_decay rate_ $\tau_{\ell mn}^{-1}$ (i.e. imaginary frequency), not the
oscillatory part $\omega_{\ell mn}$ (i.e. real frequency). This can be
confirmed analytically using Favata’s minimal waveform model (MWM; Eq. (14) of
[57]).

Figure 2: Characteristic strain $h_{\mathrm{c}}(f,\iota,\varphi_{c})\equiv
2f|\widetilde{h}_{+}-i\widetilde{h}_{\times}|$ of the primary (solid curve)
and memory (dashed curve) computed from Eqs. (3) and (II.1), seen from a fixed
direction $\iota=40\,\mathrm{deg}$ and $\varphi_{c}=0$. We consider two
fiducial non-spinning BBHs, which will also be used in the following sections:
a “light” binary (in violet) with total mass $M=2\times 10^{4}M_{\odot}$ at
redshift $z=0.5$, and an “heavy” binary (in blue) with $M=2\times
10^{5}M_{\odot}$ at $z=2$. Both BBH have mass ratio $q=1.2$. We consider the
last $25$ cycles before merger for both BBHs (which corresponds to $\sim 6$
minutes for the “light” source, and $\sim 2$ hours for the “heavy”).

Figure 2 shows the characteristic strain
$h_{\mathrm{c}}(f,\iota,\varphi_{c})\equiv
2f|\widetilde{h}_{+}-i\widetilde{h}_{\times}|$ of the primary and memory –
containing all modes up to $(\ell,|m|)=(4,4)$ – for two binaries with
different total masses and redshifts, seen from a fixed direction
$\iota=40\,\mathrm{deg}$ and $\varphi_{c}=0$. At this particular direction,
the primary characteristic strain is $\mathcal{O}(10^{2})$ greater than the
memory. Thus, the memory adds information to parameter estimation only if the
number of cycles that can be observed during inspiral is limited (c.f. Fig.
5); this may happen due to gaps in the data stream and/or confusion noise from
other sources. Indeed, whereas truncating the primary waveform at some minimum
$f_{\mathrm{in}}$ (related to the time/cycles prior to merger) significantly
reduces the SNR of the primary, the SNR of the memory is almost unchanged (as
also noted in Refs. [72, 65]). As long as the memory is observed in LISA for a
period of at least $10^{3}\,\mathrm{s}\approx 15\,\mathrm{min}$ after merger,
its SNR is approximately independent of the observation time (c.f. Fig. 6).

In this work we focus on quasi-circular and non-spinning BBHs. For the
dependence of the memory on the spins, mass ratio and eccentricity, we refer
the reader to Refs. [102, 84, 93, 60, 103, 71].

##  III Fisher Matrix and Parameter Estimation

We use a Fisher matrix analysis (commonly used in GW astronomy [104, 15, 105,
106, 107, 108]) to quantify the impact of GW memory on the parameter
estimation of the BBH source. We only consider events with merger occurring
during the mission lifetime of LISA, since the memory is mainly generated
during this phase of the binary evolution. The total signal $h=h_{0}+\delta h$
is generated in time-domain and is given by the sum of the (primary) surrogate
NRHybSur3dq8 waveform ${h}_{0}$, and the memory $\delta{h}$ computed through
the GWMemory package [93] from the primary waveform. In this work we focus on
non-spinning binaries, so the total number of parameters is reduced compared
to the general spinning case. The intrinsic parameters, defined as those that
do not depend on the relative position and orientation of the orbital plane
and the detector, are simply the redshifted mass $M_{z}$ and the mass ratio
$q\equiv M_{1}/M_{2}$ (in our convention $M_{1}>M_{2}$). On the other hand,
the extrinsic parameters are the luminosity distance $d_{\mathrm{L}}$, the
inclination angle $\iota$, and the coalescence phase $\varphi_{\mathrm{c}}$.
Thus, the input parameters for our waveforms are
$\Theta=\\{\ln{M_{z}},q,\ln{d_{\mathrm{L}}},\iota,\varphi_{\mathrm{c}}\\}$. We
only investigate the dependence on these parameters, meaning that we ignore
the particular sky position of the source and the LISA response function.
Accounting for them should not affect our results, since the LISA orbital
motion is a year time scale, whereas the longest inspirals where we see the
effect of the memory are of the order of hours.

In the strong-signal limit the probability distribution of the parameters is a
multivariate Gaussian distribution centred on the true values
$\Theta=\bar{\Theta}$, with the covariance matrix $\Sigma_{ij}$ described by
the inverse of the Fisher information matrix $\Gamma_{ij}$, up to corrections
of order of the inverse signal-to-noise ratio,

$\displaystyle\Sigma_{ij}$
$\displaystyle=(\Gamma^{-1})_{ij}[1+\mathcal{O}(\mathrm{SNR}^{-1})].$ (7)

The SNR $\rho$ and the Fisher matrix are computed as

$\displaystyle\rho^{2}$
$\displaystyle=(\tilde{h}|\tilde{h}),\qquad\Gamma_{ij}\equiv\Big{(}\frac{\partial\tilde{h}}{\partial\Theta_{i}}\Big{|}\frac{\partial\tilde{h}}{\partial\Theta_{j}}\Big{)}\Big{|}_{\Theta=\bar{\Theta}},$
(8)

where we have used the standard inner product

$(\tilde{a}|\tilde{b})\equiv
4\mathrm{Re}\int_{f_{\mathrm{min}}}^{f_{\mathrm{max}}}\mathrm{d}f\frac{\tilde{a}^{*}(f)\tilde{b}(f)}{S_{n}(f)}.$
(9)

The $S_{n}(f)$ is the sky-position- and polarization-averaged, but not
inclination-averaged, LISA power spectral density, as found in [109], with the
confusion noise corresponding to $T_{\text{obs}}=4$ years. Computing the
Fisher matrix thus allows us to find the variance of each parameter and
correlation of each pair of parameters due to measurement errors,

$\sigma^{2}_{i}=(\Gamma^{-1})_{ii},\qquad
c_{ij}=\frac{(\Gamma^{-1})_{ij}}{\sigma_{i}\sigma_{j}},$ (10)

(with no summation implied over repeated indices here).

Figure 3: The $1\sigma$ confidence ellipses computed for the primary waveform
without memory (in magenta), with memory (in blue), and with only the memory
(in green). For each of those we compare the effect of including higher modes
(solid lines) and excluding them (dashed lines). The left panel shows the
result for the “light” binary, whereas the right panel shows it for the
“heavy” binary (c.f. Fig. 2); in both cases we consider the last $25$ cycles
before merger and a line of sight $\iota=40\,\mathrm{deg}$. For the “heavy”
binary, the memory has negligible impact, since we cannot distinguish the
purple and the blue lines (the green contours of the memory fall beyond the
panel).

The total frequency-domain signal $\tilde{h}=\tilde{h}_{0}+\delta\tilde{h}$ is
given by the sum of the FT of the primary and the memory signals, which we
compute numerically via the fast Fourier transform (FFT) implemented in NumPy
[110]. In Appendix A we explain in detail our choices for manipulating the
signal, such as the padding and the window function applied, while in Appendix
C we elaborate on the numerical computation of the Fisher matrix.

Figure 4: Correlation matrices $c_{ij}$ for the “light” binary of Fig. 2. The
upper panel shows the results for the primary signal without memory, while the
lower panel is for the total waveform including memory. The major effect of
including the memory is decreasing the $(d_{\mathrm{L}},\iota)$ correlation.

##  IV Results for the distance-inclination estimation

The goal of this section is to study the impact of the memory on the parameter
estimation of intermediate-mass to supermassive BBHs with LISA, with a
particular focus on breaking the distance-inclination degeneracy. Our
motivation is twofold: (i) although nonlinear memory is expected to be
observed in single events at LISA, no attention (to our knowledge) has been
paid on its impact on the parameter estimation, and (ii) the memory signal has
an opposite correlation between distance and inclination
$c_{\ln{d_{L}},\iota}$ with respect to the primary signal, which can make it
useful to break the distance-inclination degeneracy, thus improving the
accuracy of distance estimations (c.f. Fig. 3).

Firstly, using our Fisher analysis we noted that the primary signal alone
already constrains the binary intrinsic parameters quite well, and that the
additional information provided by the memory does not constrain them further.
On the other hand, as already anticipated, we found that the dependence of the
memory on the extrinsic parameters is complementary to that of the primary
signal, and can mitigate the uncertainty on the luminosity distance
$d_{\mathrm{L}}$ and the inclination $\iota$ estimations. This is demonstrated
in Fig. 3 where we explicitly show the constraints coming from the primary
signal and the memory, as well as the effect of including all higher modes
(HMs) up to $(\ell,|m|)=(4,4)$. These results are also summarized in Tab. 1.

As it can be seen in Fig. 3, the ellipses computed from the memory signal and
the primary waveform are orthogonal, this can be intuitively expected from the
opposite monotonic dependence on the inclination of the two components, as
explained in Sec. II.2. We analytically derive the opposite correlation for
the two signals in Appendix B, where we compute the $2\times 2$ Fisher matrix
of this pair of parameters $\\{\iota,\ln d_{\text{L}}\\}$ for the dominant
mode (DM) of the primary waveform and of the memory. This property leads to a
lower correlation $c_{\iota,\ln d_{\text{L}}}$ of the overall Fisher matrix,
as shown in Fig. 4, which can mitigate the error on the luminosity distance
estimation. The estimation of the other parameters is less affected, thus
justifying our approach to focus on this particular pair of parameters.

Light |  $h_{0,\mathrm{DM}}$ |  $h_{0,\mathrm{HM}}$ |  $\delta h_{\mathrm{DM}}$ |  $\delta h_{\mathrm{HM}}$
---|---|---|---|---
SNR | 20.5 | 20.6 | 1.9 | 2.3
| $h_{0,\mathrm{DM}}$ | $h_{0,\mathrm{HM}}$ | $h_{\mathrm{DM}}$ | $h_{\mathrm{HM}}$
$\sigma_{d_{\mathrm{L}}}/d_{\mathrm{L}}$ | 0.56 | 0.55 | 0.20 | 0.18
$\sigma_{\iota}\,\mathrm{[rad]}$ | 0.75 | 0.73 | 0.27 | 0.23

Heavy |  $h_{0,\mathrm{DM}}$ |  $h_{0,\mathrm{HM}}$ |  $\delta h_{\mathrm{DM}}$ |  $\delta h_{\mathrm{HM}}$
---|---|---|---|---
SNR | 1001.4 | 1006.3 | 3.5 | 4.3
| $h_{0,\mathrm{DM}}$ | $h_{0,\mathrm{HM}}$ | $h_{\mathrm{DM}}$ | $h_{\mathrm{HM}}$
$\sigma_{d_{\mathrm{L}}}/d_{\mathrm{L}}$ | 0.0116 | 0.0108 | 0.0115 | 0.0107
$\sigma_{\iota}\,\mathrm{[rad]}$ | 0.0158 | 0.0140 | 0.0157 | 0.0139

Table 1: Summary of the results shown in Fig. 3. The first two columns
correspond to the results coming from the primary signal alone, whereas the
last two correspond to the total (including the memory) waveform. While the
memory improves significantly the estimation of the distance-inclination of
the “light” binary, it does not help much with the “heavy” binary parameters.

For signals with sufficiently large SNR or whose sources are at sufficiently
high redshifts, the uncertainty on the luminosity distance estimation becomes
dominated by weak lensing effects (see, e.g., Ref. [111, 21]), and our Fisher
analysis (which neglects lensing) ceases to be valid. However, as we show
below, the memory is helpful to parameter estimation only for binaries whose
primary signal is not very loud, and whose total redshifted mass is in the
range $10^{4}M_{\odot}\lesssim M_{z}\lesssim 10^{5}M_{\odot}$. The memory of
these sources will only be observed if they are sufficiently close $z\lesssim
1.5$, where lensing effects are not too strong. Therefore, our results
indicate that the memory has an impact on the distance estimation only in
cases where its intrinsic uncertainty (even after adding the memory) is much
larger than the contribution due to lensing.

It is natural to compare the effect of the memory with that of HMs, as the
latter can also improve the BBH parameter estimation and partially break the
aforementioned degeneracy [29, 30, 31, 32, 112, 113, 114, 115, 116, 117, 118].
Their main contribution to the SNR comes from extending the signal to higher
frequencies as compared to the dominant $(2,2)$-mode, thus being more relevant
when the merger falls well inside the sensitivity curve of the detector.
Importantly, HMs have a different dependence on the inclination angle than the
DM, which is again why they may be useful in breaking the distance-inclination
degeneracy. The memory signal, on the other hand, is always subdominant, but
can play a significant role, as we will see, when the information from the
primary at lower frequencies (i.e., from the inspiral) is absent/degraded, and
the memory becomes the only contribution at those frequencies. In the
following we explain in detail the dependence of memory-assisted parameter
estimation on the binary mass, duration of the signal, and line of sight, and
we discuss the impact of HMs.

Figure 5: Effect of the memory on the luminosity distance estimation when the
primary signal is truncated at increasing number of cycles prior to merger.
The memory becomes less important for more massive binaries, for which the
merger occurs before or close to the peak of LISA sensitivity and the memory
adds negligible contribution to the total SNR. The sources are kept at
redshift $z=0.5$ and line of sight $\iota=40\,\mathrm{deg}$.

###  IV.1 Dependence on the binary mass

By evaluating the SNR defined in Eq. (8), we find that LISA can detect the
memory (i.e., can distinguish the memory over the detector noise) when its SNR
is $>1$,555Indeed, it has been shown that this is a reasonable criteria to
ensure the distinguishability of two waveforms that differ by $\delta h$, over
the detector noise [119, 120]. It is equivalent to require that, if we
parametrize the model waveform as $h=h_{0}+\lambda\delta h$ with the fudge
parameter $\lambda\in[0,1]$, the statistical error of the fudge parameter is
$\sigma_{\lambda}<1$ and, thus, distinguishable from zero. which happens for
binary mergers with total redshifted mass in the range
$[10^{4},10^{8}]\,M_{\odot}$, thus confirming previous results [57, 70, 72].
For lower masses, the memory strain is too weak to be detected, whereas for
larger masses the turnover frequency at which the memory drops off falls below
the LISA band. Because of the dependence of the memory characteristic strain
on the binary mass, LISA will be most sensitive to the low-frequency plateau
part of the signal for light binaries, and to the subsequent high-frequency
features associated with the merger/ringdown stage for more massive binaries
(cf. Fig. 1).666See also Fig. 4 of Ref. [69]. We found that nearly all the SNR
of the memory accumulates during the merger and it is maximized for $M_{z}\sim
10^{6}M_{\odot}$, in which case the memory can be detectable up to redshift
$z\sim 14$.

For short/degraded primary signals (with fixed number of cycles), we found
that the memory is most helpful in parameter estimation for binaries with
total redshifted mass $M_{z}\lesssim 10^{5}M_{\odot}$, in which case the
memory falls in the most sensitive part of the LISA sensitivity band, whereas
the merger covers only the high-frequency edge of the spectrum. For this
reason, the hierarchy of the SNR between the memory and the primary signals is
less severe than for more massive binaries $M_{z}\gtrsim 10^{6}M_{\odot}$,
whose mergers occur in the middle/low-frequency part of the LISA band.

The difference between “light” and “heavy” binaries is clear from Fig. 3,
which shows the effect of considering short signals (in this case, the last 25
cycles) for two different binary masses. In the left panel, for the “light”
binary, there is a manifest reduction in the parameter uncertainties, which is
due to the intersection of the primary and memory confidence ellipses that
results in the shrinking of the total signal’s confidence ellipse. Such an
effect is not present in the right panel, for the “heavy” binary, because, due
to the large SNR of the primary signal, its confidence ellipse is already
entirely enclosed within the memory’s confidence ellipse. The main factor
controlling the relative sizes of the confidence ellipses, and thus the impact
of the memory in parameter estimation, is the ratio between the memory and the
primary signal SNRs (cf. Fig. 5).

Figure 6: Upper panel: the effect of the memory on the estimation of the
luminosity distance as a function of the inclination. The parameters of the
source are those of the “light” binary of Fig. 3 with mass-ratio $q=1.2$, and
the number of cycles before merger is $N_{\text{cycles}}=25$. In green we show
the relative uncertainties $\sigma_{d_{\mathrm{L}}}$ computed from the
waveform without the memory $h_{0}$, with only the dominant mode DM and with
higher modes HM. In purple the same but for the total waveform $h=h_{0}+\delta
h$; in the DM case the memory is computed using only the dominant modes. Lower
panel: the $\rho$-ratio required to achieve a 10% improvement in the
luminosity distance estimation (i.e., such that
$\sigma_{d_{\mathrm{L}},\mathrm{wm}}/\sigma_{d_{\mathrm{L}}}=0.9$) as a
function of inclination for two different values of the mass-ratio
$q=\\{1.2,3\\}$.

###  IV.2 Dependence on the signal duration

As just discussed, our results show that the most reliable indicator for how
much the memory contributes to constrain the binary parameters (for a fixed
line of sight) is the ratio between the SNR of the memory and of the primary
signal, $\rho_{\mathrm{m}}/\rho_{0}$. This ratio depends on the total mass of
the binary (cf. Tab. 1), but it is also a function of the inclination and the
time duration of the data taken before the merger. Decreasing the mass and the
duration of the signal prior to merger tends to increase this ratio. We
quantify the improvement in parameter estimation in terms of the ratio between
the standard deviation of the luminosity distance with memory
$\sigma_{d_{\mathrm{L}},\mathrm{wm}}$ and without $\sigma_{d_{\mathrm{L}}}$ –
the memory contributes most when the ratio
$\sigma_{d_{\mathrm{L}},\mathrm{wm}}/\sigma_{d_{\mathrm{L}}}$ is minimized.

Our results are presented in Fig. 5 for different masses in the range
$[10^{4},10^{5}]M_{\odot}$. In the upper panel, we show the $\sigma$-ratio as
a function of the number of cycles $N_{\text{cycles}}$ observed prior to
merger. We note that for any given number of observed cycles, the memory of
more massive binaries leads to a smaller improvement in distance estimation as
compared to lighter binaries (as explained in Sec. IV.1). In the lower panel,
we show that there exists a monotonic, one-to-one relationship between the
$\sigma$-ratio and the $\rho$-ratio, for a fixed inclination, which is
approximately insensitive to the mass of the binary. This observation will
allows us to determine, in the next section, the “critical” $\rho$-ratio
needed to achieve a given improvement in the distance estimation as a function
of the inclination.

###  IV.3 Dependence on the inclination

Fixing the other parameters, the impact of the memory depends strongly on the
inclination of the binary, as we show in the upper panel of Fig. 6. In this
plot we present the relative uncertainty in the luminosity distance at
different inclination angles, considering the “light” binary (discussed
previously) with mass-ratio $q=1.2$ and fixing the observed signal to 25
cycles prior to merger. We restrict the plot to values of $\iota\in[0,\pi/2]$,
since its behavior is symmetric with respect to $\iota=\pi/2$. Note that, in
the absence of the memory and the HMs, the uncertainty on the luminosity
distance diverges for face-on $\iota\to 0$ (and face-off $\iota\to\pi$)
configurations. As we discuss in App. B, this is due to the fact that the
Fisher matrix is singular at this inclination and, thus, the Fisher analysis
for the primary DM ceases to be valid in a neighborhood of $\iota=0$ (and of
$\iota=\pi$). However, for unequal-mass binaries the HMs regularize the Fisher
matrix, and even a slightly asymmetric binary as $q=1.2$ has a relative
uncertainty in the distance smaller than $\sim 0.8$ for all inclination
angles.

In the upper panel of Fig. 6 we compare the impact of HMs on the distance
estimation with that of the memory for different inclination angles, keeping
the mass-ratio $q=1.2$ fixed. However, we also verified that HMs become more
relevant the more asymmetric the binary is, whereas the impact of the memory
is larger for smaller $q\sim 1$.777The influence of the mass-ratio on the
impact of the memory can be understood from Eq. (6), $\delta
h^{(0\textrm{PN})}\propto q/(1+q)^{2}$, which results in the memory
characteristic strain  $f\widetilde{\delta h}^{(0\textrm{PN})}\propto
q/(1+q)^{2},$ whereas the primary characteristic strain is [15]
$f\widetilde{h}_{0}^{(0\textrm{PN})}\propto\sqrt{q}/(1+q).$ The large effect
of the memory observed in the plot – improving the distance estimation by more
than a factor of $4$ for some inclinations – is partially related to the short
duration of the signal considered (truncated at 6 minutes prior to merger); as
discussed in Sec. IV.2, for longer signals the information in the primary
(which accumulates over the inspiral) eventually becomes dominant, with the
memory contributing negligibly to parameter estimation (cf. Fig. 5).

In the lower panel of Fig. 6 we show the critical $\rho$-ratio needed to
achieve a $10\%$ reduction in $\sigma_{d_{\mathrm{L}}}$ (i.e., such that
$\sigma_{d_{\mathrm{L}},\mathrm{wm}}/\sigma_{d_{\mathrm{L}}}=0.9$) as function
of the inclination, for two different mass-ratios $q=\\{1.2,3\\}$. From the
discussion in Sec. IV.2, these curves are approximately independent of the
binary’s total mass. Moreover, we see from this plot that they are also only
mildly dependent on the mass-ratio $q$. Another important observation is that
closer to face-on a smaller $\rho$-ratio is needed to achieve a given
improvement in the distance estimation compared to edge-on. So, close to face-
on the memory can be relevant even if the primary is observed for relatively
long periods ($\sim$ few hours). For example, if our “light” binary (with
$q=1.2$) is seen close to face-on, the inclusion of the memory information
leads to a $10\%$ reduction in $\sigma_{d_{\mathrm{L}}}$ if the inspiral is
observed over less than 6 hours prior to merger.

###  IV.4 Impact of higher modes

Here we discuss the contribution of HMs to the SNR of the primary and memory
signals, and its consequent influence on the estimation of the luminosity
distance. We consider the “light” and “heavy” binaries of Fig. 2, both of them
with mass ratio $q=1.2$ and seen at an intermediate inclination angle
$\iota=40\,\mathrm{deg}$. The outcomes are shown in Fig. 3 and Tab. 1. These
results show that, for fixed mass-ratio, the relative increase of the SNR due
to the inclusion of HMs is always larger for the memory than for the primary
signal. However, the overall SNR of the memory decreases faster with
increasing mass-ratio as compared to the SNR of the primary. Regarding
parameter estimation, Tab. 1 shows that HMs have a greater impact for the
“heavy” binary than for the “light” one. This can be understood from the fact
that the HMs add information mostly close to merger (i.e., at high
frequencies), which is therefore masked for the “light” binary. This is
opposed to the memory which, as we have seen, contributes the most to
parameter estimation for the “light” binary. Moreover, we confirmed that,
although the overall SNR is suppressed with increasing mass-ratio, the
increase in power on HMs for more asymmetric binaries can considerably improve
the distance-inclination estimation with the primary waveform, reducing the
impact of including the memory component; thus, the equal-mass scenario
presented in Tab. 1 corresponds to the maximum impact of the memory in
parameter estimation over mass-ratio. For instance, for the parameters of the
“light” binary of Tab. 1, but with $q=4$, the uncertainty on the luminosity
distance is $\sigma_{d_{\text{L}}}/d_{\text{L}}=0.93$ using only the DM of the
primary waveform, which is improved to
$\sigma_{d_{\text{L}}}/d_{\text{L}}=0.48$ considering the HMs of the primary,
and to $\sigma_{d_{\text{L}}}/d_{\text{L}}=0.38$ by including the memory.
Interestingly, even though it has a smaller $\rm{SNR}=12.53$, this binary is
better localized in space than the equal-mass one. For the parameters of the
“heavy” binary, but with $q=4$, the uncertainty
$\sigma_{d_{\text{L}}}/d_{\text{L}}=0.017$ is improved to
$\sigma_{d_{\text{L}}}/d_{\text{L}}=0.008$ with the inclusion of HMs in the
primary waveform, and, as for the close to equal-mass case, we do not see
further improvements by adding the memory.

Therefore we conclude that, in some cases, the information contained in the
memory and in the HMs is complementary for parameter estimation, due to their
different dependence on the mass-ratio and total mass.

##  V Population forecasts

Figure 7: The redshift and total mass distributions of mergers with observable
memory $\rho_{\mathrm{m}}>1$ for the various HS population models. The models
with “delays” have mergers typically at lower redshift, explaining the
corresponding higher fraction of events with observable memory (cf. Table 2).
The different redshift distribution translates into slightly different peak
locations in the total mass distributions.

Figure 8: Left panel: the mass-ratio distribution for the events of Fig. 7.
Right panel: the inclination angle distribution of the events with detectable
memory of the “SN short delays” model (with $N_{\text{th}}=418$) as compared
to the prior isotropic (cosine) distribution.

In this section we study the detectability of the nonlinear memory for
realistic population models of massive black holes, and assess its potential
impact on parameter estimation considering the presence of gaps in the data
stream. We update the previous forecasts of Refs. [72, 70] for the
measurability of the memory in single events with space-based detectors, by
using the recent population models described in Refs. [78, 79] (and recently
used in Ref. [83]).

###  V.1 Framework

In our analysis we consider all 8 models described in Refs. [78, 79], each of
those corresponding to different assumptions about some of the main
uncertainties in the cosmological evolution of massive BHs. These involve [78,
79]: (i) the high-redshift mass function of the “seeds” of the massive black
hole population [“light seeds” (LS) originating from population III stars, or
“heavy seeds” (HS) originating from direct collapse of protogalactic gaseous
disks], (ii) the time delay between the galaxy merger and the corresponding
BBH mergers (realistic “delays”, or “short delays” neglecting the contribution
from scales of the order of hundreds of pc), and (iii) the presence or not of
supernova feedback (“SN” or “noSN”) on the accretion disk of massive black
holes. For each of the 8 models we compute the memory random realizations of
mergers corresponding to 4 years of LISA mission, and present average results
over many such realizations.

We consider only the final $20$ cycles prior to merger, which is enough to
give us a reliable estimate of the SNR of the memory (cf. Sec. II.2). Due to
the limited parameter space covered by the waveforms NRHybSur3dq8, we restrict
the mass-ratio to $q\leq 8$ (i.e., we artificially fix $q=8$ for every merger
with higher values of $q$). This assumption is stronger for the LS case than
for the HS one, since the two cases have quite different mass-ratio
distributions, with the latter having a sharper peak close to equal mass (see
Fig. 11 of Ref. [78]). Since the waveforms NRHybSur3dq8 do not cover
precessing BBHs, we consider only the spin components orthogonal to the
orbital plane and restrict them to $|\chi_{i,\hat{z}}|\leq 0.8$,
$i\in\\{1,2\\}$.888Using the waveforms NrSur7dq2 which allow for precessing
BBHs, we found that introducing the spin components along the orbital plane
generally leads to an increase of the memory SNR. We removed from the catalogs
the sources with $M_{z}\geq 10^{8}\,M_{\odot}$, because evaluating their SNR
is computationally very expensive and, as discussed in Sec. IV.1, they fall
outside the parameter space of interest for memory observation with LISA. We
take the inclination angle and the coalescence phase to be uniformly
distributed in $\cos{\iota}\in[-1,1]$ and $\varphi_{\mathrm{c}}\in[0,2\pi]$.

Astrophysical Catalogues
---
|  Light seeds |  Heavy seeds
SN-delays |  $N_{\mathrm{tot}}=47$ |  $N_{\mathrm{tot}}=27.3$
|  $N_{\mathrm{th}}=0.4\,(0.1)$ |  $N_{\mathrm{th}}=21.2\,(10)$
|  $\langle\rho\rangle=0.04$ |  $\langle\rho\rangle=6$
|  $\rho_{\mathrm{max}}=7$ |  $\rho_{\mathrm{max}}=97$
noSN-delay |  $N_{\mathrm{tot}}=191$ |  $N_{\mathrm{tot}}=10$
|  $N_{\mathrm{th}}=6\,(1)$ |  $N_{\mathrm{th}}=7.5\,(4)$
|  $\langle\rho\rangle=0.17$ |  $\langle\rho\rangle=6.9$
|  $\rho_{\mathrm{max}}=11.64$ |  $\rho_{\mathrm{max}}=68.7$
SN-short |  $N_{\mathrm{tot}}=149$ |  $N_{\mathrm{tot}}=1245$
Delays |  $N_{\mathrm{th}}=1\,(1)$ |  $N_{\mathrm{th}}=418\,(33)$
|  $\langle\rho\rangle=0.04$ |  $\langle\rho\rangle=1$
|  $\rho_{\mathrm{max}}=5.01$ |  $\rho_{\mathrm{max}}=43$
noSN-short |  $N_{\mathrm{tot}}=1203$ |  $N_{\mathrm{tot}}=1251$
Delays |  $N_{\mathrm{th}}=12\,(2)$ |  $N_{\mathrm{th}}=392\,(29)$
|  $\langle\rho\rangle=0.06$ |  $\langle\rho\rangle=1.1$
|  $\rho_{\mathrm{max}}=17$ |  $\rho_{\mathrm{max}}=51$
Table 2: SNR of the memory for the astrophysical models of Refs. [78, 79]. For
each model we consider a random realization of 4 years of events. We denote
the total number of events by $N_{\mathrm{tot}}$ and the number of those with
SNR above the threshold value $\rho_{\mathrm{m}}>1\,(5)$ by $N_{\mathrm{th}}$.
The $\langle\rho\rangle$ and $\rho_{\mathrm{max}}$ are, respectively, the
average and the maximum $\rho_{\mathrm{m}}$ in the particular realization of
events. For the SN-delays models and the NoSN-delay HS model, the reported
values are the averages over 10 realizations of 4-year events, since for these
models the total number of mergers are much smaller than for the others.

###  V.2 Detectability of memory in single events

In the different 4-year realizations, we looked for events with SNR of the
memory above the threshold value $\rho_{\mathrm{th}}\equiv 1$. Table 2
summarizes our results. For each population model we denote by
$N_{\mathrm{tot}}$ the total number of events considered in the 4-year
realizations, and by $N_{\mathrm{th}}$ the number of those with memory SNR
above the threshold $\rho_{\mathrm{th}}=1$ (inside parentheses
$\rho_{\mathrm{th}}=5$). We also indicate the average SNR of the memory
$\langle\rho\rangle$ and its maximum value $\rho_{\mathrm{max}}$ in the
particular 4-year realization of events.

The number of events with significant memory SNR depends strongly on the
astrophysical population model, with clear differences between the LS and the
HS models. Our results are especially promising for the HS scenario, as it
suggests that about $75\%-78\%$ of events for the “delays” model and
$31\%-33\%$ for the “short-delays” model will have observable memory with
$\rho_{\text{m}}>1$. The numbers inside the parenthesis can be directly
compared with the results of Ref. [72], where the “Q3d” and “Q3nod” models
considered there (and presented in Ref. [121], based on Refs. [80, 81, 82])
can be compared, respectively, with “delays” and “short-delays” HS models. For
those, we found that about $36\%-40\%$ and $25\%$, respectively, of the total
events have detectable memory with $\rho_{\mathrm{m}}>5$, as opposed to
$3.7\%$ and $1\%$ for the “Q3d” and “Q3nod” models. The main reason for this
large mismatch is that the new population models of Refs. [78, 79] that we
used in this work have a different (more realistic) delay model, which shifts
the mergers to lower redshift.

Figure 7 shows the distribution of the redshift and the total mass of the
events with $\rho_{\mathrm{m}}>1$ for the various HS models. The distribution
of models with “delays” peaks at lower redshift, while that of “short-delays”
models extends up to $z\sim 13$, which explains why the former have a bigger
fraction of events with memory SNR above the threshold than the latter.
Interestingly, we found some events with particularly high SNR
($\rho_{\mathrm{m}}\gtrsim 50$), as can be seen in Tab. 2; these belong to the
low redshift tail of the distributions ($z\lesssim 1$). The location of the
peak of the total mass distribution changes slightly for the various models,
but it is such that the total redshifted mass is about $\sim 10^{6}M_{\odot}$.
In Fig. 8 we show the mass-ratio distribution for the same events of Fig. 7.
Most of the events with detectable memory have a mass-ratio close to unity
with a sharp suppression at higher values, so that the restriction to $q<8$
turns out not to affect our results for the HS models. The “noSN delays” model
is the only one presenting a mild accumulation of events with $q=8$ due to
this restriction.

The situation is quite different for the LS models, which have a much broader
distribution in the mass-ratio, resulting in a substantial (fictitious)
accumulation of events at $q=8$. So, we repeated our analysis removing
directly the binaries with $q>8$ from the catalog. In this more conservative
approach, for “noSN-short delays” we noted a reduction from 12 to 9 events
with $\rho_{\mathrm{m}}>1$, but no change in the number of events with
$\rho_{\mathrm{m}}>5$. For “SN-short delays” and “SN delays” models we found
no events with detectable memory, and for “noSN-short delays” a reduction from
6 to 4 events with $\rho_{\mathrm{m}}>1$, and no events at all with
$\rho_{\mathrm{m}}>5$. Therefore, for LS population models our results
indicate that the prospects of observing the memory with LISA do not seem
promising.

We repeated our analysis of LS population models also for the future
generation ground-based detectors Cosmic Explorer (CE) [68] and Einstein
Telescope (ET) [122], which have better sensitivity at higher frequencies, and
thus to lower masses.999We used the sensitivity curve of the configurations
ET-D of Ref. [123] and CE_40km_lf of
https://cosmicexplorer.org/sensitivity.html, which have the best sensitivity
at low frequencies. We have found almost no events with observable memory
$\rho_{\text{m}}>1$ in 4 years of observation (there was just one event with
$\rho_{\mathrm{m}}\simeq 2$ for the “SN short delays” model with ET), and an
average memory SNR within $10^{-1}-10^{-3}$.

###  V.3 Impact of the memory on distance estimation

Given that HS models predict such a large number of events with observable
memory at LISA (which can be almost up to $80\%$ of the total number of
events, in the most optimistic scenario), we consider here the impact of the
memory on the luminosity distance estimation for these sources. As we have
shown in Sec. IV.2 the impact of including the memory is highly dependent on
the ratio of the SNR of the memory and the primary signals, with the memory
helping substantially to constrain the distance when the information (or the
duration) of the primary signal is somehow limited. Among other possibilities,
this could happen due to the presence of gaps in the data stream, which causes
a partial loss of signal.

Two kinds of gaps are expected at LISA: the scheduled ones, related to the
regular maintenance of the detector, and the unscheduled ones, due to
unexpected problems/events. In Ref. [83] it was shown that the scheduled gaps
have little or no impact at all on the parameter estimation of massive BBHs,
but the unscheduled ones could degrade significantly the parameter estimation.
Thus, there is the intriguing possibility that, in the presence of unscheduled
gaps, the memory may add useful information to constrain the binary
parameters.101010Note that the LISA data analysis will be further complicated
by the presence of many overlapping signals [124].

To quantify this effect we consider the particular gap model used in Ref.
[83], which is consistent with a 75% duty cycle, as expected for LISA [77]. We
simulate the presence of gaps by windowing the signals as in Ref. [83],
considering scheduled gaps with a typical duration of 3.5 hours every week and
unscheduled gaps with a duration of 3 days. The time interval between two
unscheduled gaps is treated as a random variable following an exponential
probability distribution $p(\Delta T)=\lambda\exp{(-\lambda\Delta T)}$ with
$1/\lambda=9$ days. With these choices, we simulate the effective data taking
of the mission and we distribute the merger times uniformly over the 4-year
mission duration. From our study in Sec. IV, we expect the memory to be
helpful in constraining the binary (extrinsic) parameters for a particular
chunk of data if the merger happens within the first few hours from the last
gap.

For concreteness, let us compute the (average) total number of mergers
occurring within 6 hours from the last gap. The number of unscheduled gaps can
be estimated by $N_{\text{gap}}\sim T_{\text{mission}}/T_{\text{gap}}$, where
$T_{\text{gap}}$ is the sum of the average time interval between gaps and the
gap duration, $\langle\Delta T\rangle+3\simeq 12$ days, thus
$N_{\text{gap}}\sim 120$. We focus on the “SN-short delays” HS population
model, the most optimistic scenario with the highest number of events with
observable memory, $N_{\text{th}}=418$; note that due to the presence of gaps
this number is reduced by $75\%$. Thus, we can estimate the number of mergers
by multiplying the probability of having at least one merger in 6 hours by
$0.75N_{\text{th}}N_{\text{gap}}$, which gives $\sim 6.4$ events. We checked
this result numerically by simulating 50 times the distribution of mergers
over the gap realization and we found consistent results.

To find the number of events for which the inclusion of the memory decreases
by more than 5% the uncertainty on the luminosity distance, we computed the
$\rho-$ratio (i.e., $\rho_{\textrm{m}}/\rho_{0}$) for each event occurring
within 6 hours from the last gap in 50 numerical realizations, neglecting the
information accumulated in the inspiral before the gap. Subsequently, we
compared those $\rho$-ratios with the critical values needed to achieve a 5%
improvement on the luminosity distance estimation, which depend on the
particular binary inclination (as in the lower panel of Fig. 6, but for
$\sigma_{d_{L}\text{wm}}/\sigma_{{d_{L}}}=0.95$). We found that, on average,
only 0.14 events of the 6.4 occurring close after a gap have an improvement of
more than 5% on the distance estimation from including the memory; this
corresponds to 0.04% of the total number of events with observable memory in
this population model.

We believe that this low value is due to the fact that most of the BBH mergers
with observable memory correspond to configurations relatively close to edge-
on (c.f. right panel of Fig. 8), where the critical $\rho-$ratio is much
higher. Another reason is that the majority of the BBHs with observable memory
in the population considered have a redshifted total mass $\gtrsim
10^{6}M_{\odot}$, whereas as discussed in sec IV.1 the memory is more helpful
for lighter binaries $M_{z}\lesssim 10^{5}M_{\odot}$. The “noSN-short delays”
HS population model, which is the second most optimistic in terms of number of
events with detectable memory, suffers from these same issues and is, thus,
expected to give a similarly small result. The other population models have
far fewer events with observable memory, thus it is very unlikely that any of
these mergers will happen sufficiently close to a gap to have a sufficiently
large $\rho-$ratio.

In summary, applying our Fisher analysis to state-of-the-art synthetic
catalogs of massive BBHs indicates that the memory will not help constraining
further the binary parameters at LISA, even in the presence of gaps in the
data stream. However, there is substantial uncertainty on the assumptions
adopted in this analysis, in particular, regarding the population and gap
models, and we cannot exclude the possibility that there may exist additional
effects leading to a larger degradation of the primary signal than those
considered here.

##  VI Conclusion

In this work we have investigated the prospects of using the nonlinear GW
memory to help infer the parameters of merging BBHs. In particular, we have
focused on massive BBHs detections with the future space-based interferometer
LISA, as these are the most promising individual sources of memory. Our
motivation is to use the additional source of information provided by the
memory signal to break the degeneracy between inclination $\iota$ and
luminosity distance $d_{\mathrm{L}}$, which is present in the leading-order GW
signal. This is especially important for attempts to use these BBHs as
standard sirens (either via statistical identification of the host galaxy
[125], or possibly using an electromagnetic counterpart due to the merger
taking place in a gas-rich environment [21, 24, 25]), as the uncertainty on
the Hubble constant $H_{0}$ crucially depends on the uncertainty on
$d_{\mathrm{L}}$.

We find that the memory can indeed play a significant role in breaking this
inclination–luminosity distance degeneracy. This occurs in cases where the
redshifted total mass is relatively small ($\lesssim 10^{5}\,M_{\odot}$), the
binary is seen not very close to edge-on, and the observation time is limited
to a few hours prior to merger. The limitation on the observation time could
occur due to, e.g., gaps in the data stream caused by interferometer downtime,
or confusion noise from the presence of many other simultaneous signals in the
LISA frequency band.

In order to understand the relevance of these results for the LISA mission, we
started by performing a population study using new synthetic catalogs of
massive BBHs to forecast the number of BBH events with observable memory
($\rho_{\text{m}}>1$). While there are currently large theoretical
uncertainties on the astrophysical processes leading to these mergers, we find
a substantially larger number of events with significant memory as compared to
previous forecasts. The prospects are particularly bright for the heavy seed
model with “short delays” [78, 79], which presents about 400 memory events for
a 4-year mission time. On the other hand, most of the mergers coming from
light seed models [78, 79] are undetectable by LISA (and so is their memory).

Finally, we considered a commonly used gap model, which includes both the
scheduled and unscheduled types, to quantify the benefit of the memory in the
estimation of the luminosity distance. For the most optimistic “short delays”
heavy seed models [78, 79], we found that, out of the $\sim 0.75\times 400$
observable memory events in a 4-year mission time, just 0.14 events will
produce a larger than 5% decrease in $\sigma_{d_{L}}$. Thus, our analysis
indicates that the information in the memory signal will not help constraining
further the binary parameters at LISA, even in the presence of gaps in the
data stream. This is due to the fact that most of the events with observable
memory are seen close to edge-on, in which case the luminosity distance and
inclination are only slightly correlated in the primary signal and, thus, the
information added by the memory is negligible for parameter estimation.

Our study, based on a Fisher matrix analysis, could be further improved by
performing a full Bayesian analysis and by investigating the effect of priors
on the luminosity distance estimation, which is especially important when the
parameters are not well constrained. Another interesting extension of our work
would be to consider the impact of the memory on parameter estimation for
binaries with precession. However, we expect that our key finding — that the
memory signal can only play an important role in BBH parameter estimation when
there is limited information from the inspiral — holds generically, due to the
different orders of magnitude of the primary and memory signal characteristic
strains. We also leave open the possibility that some currently unforeseen
effects may lead to a much larger degradation of the primary signal than the
one due to the presence of gaps, which could make the memory information more
relevant to parameter estimation.

As a final remark, we note that even if the information in the memory turns
out not to be very useful in constraining the binary parameters, the amount of
events with detectable memory we found for LISA (c.f. Table 2) suggests that
it may still play a significant role as a test of GR in the strong-gravity
nonlinear regime, since most of the memory is generated close to merger. We
leave these questions for future work.

###### Acknowledgements.

The authors would like to thank Juan Calderon Bustillo, Xisco Jimenez Forteza,
Giada Caneva and Marc Andrés for their technical help in the first stage of
this project. We are also grateful to Neil Cornish for his valuable comments
on a draft of the paper. In this study we used the software packages
matplotlib [126], numpy [110], scipy [127], LISA Sensitivity [109], gwmemory
[62], GWsurrogate [128], surfinBH [98], and qnm [97]. RV is supported by grant
no. FJC2021-046551-I funded by MCIN/AEI/10.13039/501100011033 and by the
European Union NextGenerationEU/PRTR. RV also acknowledges support by grant
no. CERN/FIS-PAR/0023/2019. DB is supported by a ‘Ayuda Beatriz Galindo
Senior’ from the Spanish ‘Ministerio de Universidades’, grant BG20/00228. The
research leading to these results has received funding from the Spanish
Ministry of Science and Innovation
(PID2020-115845GB-I00/AEI/10.13039/501100011033). IFAE is partially funded by
the CERCA program of the Generalitat de Catalunya. This work was partly
enabled by the UCL Cosmoparticle Initiative. EB acknowledges support from the
European Union’s H2020 ERC Consolidator Grant “GRavity from Astrophysical to
Microscopic Scales” (Grant No. GRAMS-815673) and the EU Horizon 2020 Research
and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement No.
101007855.

##  Appendix A Signal processing

In this section we provide more details about our choices in manipulating the
BBH waveforms. We first generate the primary signal with a sampling time
$\Delta t=1/4\,\mathrm{s}$, and we subsequently generate its memory via the
GWmemory package [93]. As explained in Sec. III, we compute the total signal
in frequency domain, summing the individual FFTs of the primary waveform and
of the memory. However, we find that a spurious contribution of the primary
waveform at frequencies $f<f_{\text{in}}$ generates cross-terms between the
primary and the memory signal of order $\mathcal{O}(h_{\mathrm{c}})$ which
affect the computation of the SNR and the Fisher matrix. To prevent these
artifacts from affecting our results, we removed the contribution from
$f<f_{\text{in}}$ of the primary signal before summing the individual FFTs.

We follow a standard procedure to manipulate the primary waveform, namely,
applying a window function, padding the signal, and taking the FFT. We apply
the following window function to the primary signal:

$z(t)=\frac{1}{4}\Big{[}1+\tanh\big{(}\tfrac{t-t_{0}}{\sigma_{0}/4}\big{)}\Big{]}\Big{[}1-\tanh\big{(}\tfrac{t-t_{\mathrm{h}}}{\sigma_{\mathrm{h}}/4}\big{)}\Big{]},$
(11)

with $Mt_{0}=150$ and $Mt_{\mathrm{h}}=110$, respectively, at the beginning
and at the end of the time series. The duration of the windowing is set by
$M\sigma_{0}=50$ and $M\sigma_{\mathrm{h}}=20$.

For the memory signal we follow a different procedure. We first extend the
generated $\delta h(t)$ (evaluated numerically for $t_{0}\leq t\leq
t_{\mathrm{f}}$) at the beginning and the end with the constant values $\delta
h(t_{0})$ and $\delta h(t_{\mathrm{f}})$, respectively, using the same padding
length as for the primary signal. Subsequently, we apply the following window
(as in Ref. [129]):

$w(t)=\begin{cases}1,\qquad&t-t_{\mathrm{d}}<0\\\
\frac{1}{2}\big{(}1+\cos[2\pi
f_{\mathrm{d}}(t-t_{\mathrm{d}})]\big{)},\qquad&t-t_{d}\geq 0\\\
0,\qquad&t-t_{\mathrm{d}}\geq\frac{1}{2f_{\mathrm{d}}}\end{cases}$ (12)

and choose $Mt_{\mathrm{d}}=140$. The choice of the decay frequency of the
window function $f_{\mathrm{d}}$ greatly impacts the spectral shape of the
memory at low frequencies, as can be seen in Fig. 9, where we show the
characteristic strain $h_{c}$ of the memory for different values of
$f_{\mathrm{d}}\in\\{10^{-2},10^{-3},10^{-4},10^{-5}\\}\,\mathrm{Hz}$. Note
that higher values of $f_{\mathrm{d}}$ inject spurious power at the
frequencies for which LISA is most sensitive, thus leading to an artificial
increase of the respective memory SNR
$\rho_{\mathrm{m}}=\\{7.37,5.35,5.13,5.12\\}$. Thus, while taking a lower
$f_{d}$ is more reliable, in the sense that it does not overestimate the
memory SNR, it implies a corresponding longer observation time of the memory,
which can become inconsistent with the maximum observation time considered for
the primary signal. However, this does not pose a real problem since most of
our analysis applies to cases where the primary signal is observed for a few
hours, whereas the SNR of the memory does not greatly change as long as the
memory is observed for more than 15 minutes ($f_{d}\lesssim
10^{-3}\,\mathrm{Hz}$).

In computing the SNR and the Fisher matrix we take $f_{\text{min}}=1/T$, where
$T$ is the total length of the signal, and
$f_{\mathrm{max}}=\mathrm{min}\\{1\,\mathrm{Hz},f_{440}\\}$, since we find
that the QNM $f_{440}\equiv\omega_{440}/2\pi M_{z}$ is a good measure of the
maximum frequency present in the signal (note that our waveforms include modes
up to $\ell_{\mathrm{max}}=4$). For total masses between
$[10^{4},10^{5}]M_{\odot}$ we find that fixing the minimum frequency at
$f_{\text{min}}=10^{-4}\,\mathrm{Hz}$ does not change our SNR and Fisher
forecasts.

Figure 9: Characteristic strain $h_{\mathrm{c}}(f,\iota,\Phi)$ of the memory
for different choices of decay frequency of the window function (12). The
parameters of the binary are the same as those of the “light” binary in Fig.
2, but with an inclination $\iota=90\,\mathrm{deg}$.

##  Appendix B Analytic considerations on the distance-inclination Fisher
matrix

Here we review the Fisher matrix derivation of the $(d_{\mathrm{L}},\iota)$
degeneracy by computing the relative $2\times 2$ matrix analytically.
Including other parameters have little effect close to the degenerate points,
since the main source of error comes from this submatrix. Subsequently, we
show that the results represented in Fig. 3 can be understood by taken into
account simply the main angular dependence of the primary and the memory
signals. In order to estimate the effect of the memory on this degeneracy it
is enough to focus on the part of the Fisher matrix regarding the extrinsic
parameters $\\{d_{\mathrm{L}},\cos(\iota),\phi,\varphi_{\mathrm{c}}\\}$ since,
at linear order in $\mathrm{SNR}^{-1}$, it is decoupled from the one
associated to the intrinsic parameters
$\\{M_{z},q,t_{\mathrm{c}},\psi,\mathrm{Spins}\\}$ [15], where $\psi$ is the
polarization angle. Indeed, while the intrinsic parameters are mainly
extracted from the phase evolution of the waveform, the extrinsic parameters
depend on the amplitudes $h_{+}$ and $h_{\times}$ [16]. Here we ignore the
dependence on the coalescence phase $\varphi_{\mathrm{c}}$, since (at leading
order) it does not affect the memory. At Newtonian (0 PN) order the primary
waveform is

$\displaystyle h_{+,0}$ $\displaystyle=2\frac{\eta
M_{z}}{d_{\mathrm{L}}}[M\omega(t)]^{\frac{2}{3}}(1+\cos^{2}\iota)\cos[2\varphi(t)],$
(13) $\displaystyle h_{\times,0}$ $\displaystyle=4\frac{\eta
M_{z}}{d_{\mathrm{L}}}[M\omega(t)]^{\frac{2}{3}}\cos\iota\sin[2\varphi(t)],$
(14)

using the polarization conventions of Ref. [91]. The GW amplitude in the
detector can be written in the frequency domain as [109]

$\tilde{h}(f)=F^{+}(f)\tilde{h}_{+}(f)+F^{\times}(f)\tilde{h}_{\times}(f),$
(15)

where $F^{+,\times}(\iota,\phi,\psi,f)$ are the frequency-dependent detector
response functions, which also depend on the source sky-location and
polarization angle. Substituting the primary waveform FTs in the last
expression we find

$\tilde{h}_{0}=\frac{\kappa_{0}}{d_{\mathrm{L}}}\Big{\\{}F^{+}(1+\cos^{2}\iota)-2iF^{\times}\cos\iota\Big{\\}},$
(16)

where $\kappa_{0}$ is independent of both the luminosity distance (for fixed
$M_{z}$) and the inclination angle. The sky- and polarization-averaged Fisher
matrix (8) has then the form111111Where we used $\langle
F^{+}(f)F^{\times*}(f)\rangle=0$ for the sky- and polarization-averaging of
the cross terms [109].

$\Gamma^{0}_{ij}=\left(\frac{\rho_{\kappa_{0}}}{d_{\mathrm{L}}}\right)^{2}\hat{\Gamma}^{0}_{ij},$
(17)

with $i,j\in\\{\log{d_{\mathrm{L}}},\iota\\}$, where
$\rho_{\kappa_{0}}^{2}=(\kappa_{0}|\kappa_{0})$ and the matrix

$\hat{\Gamma}^{0}=\begin{pmatrix}(1+\cos^{2}\iota)^{2}+4\cos^{2}\iota&(3+\cos^{2}\iota)\sin(2\iota)\\\
(3+\cos^{2}\iota)\sin(2\iota)&4(1+\cos^{2}\iota)\sin^{2}\iota\end{pmatrix}\,,$

depends only on the inclination. This matrix is clearly singular for face-
on/off binaries, and it is diagonal for edge-on ones (implying that the two
parameters are uncorrelated),

$\displaystyle\hat{\Gamma}^{0}(\iota\in\\{0,\pi\\})=\begin{pmatrix}1&0\\\
0&0\end{pmatrix},\qquad\hat{\Gamma}^{0}(\iota=\tfrac{\pi}{2})=\begin{pmatrix}1&0\\\
0&1\end{pmatrix}.$

It is easy to see that for inclination angles $0<\iota<\frac{\pi}{2}$
($\frac{\pi}{2}<\iota<\pi$) the distance and the inclination are _negatively_
(_positively_) correlated.

This result shows that the degeneracy we focused in in this work is driven by
the dependence of the (leading) quadrupole waveform on the inclination, and
the particular combination of plus and cross populations measured by the
detector, which, in particular, lead to a singular Fisher matrix for face-
on/off configurations. This is a well-known issue in the literature and
special care must be taken close to the singular points, where one should use
a beyond-Gaussian analysis [107, 106]. However, for these face-on/off
configurations the memory is almost vanishing, so that in this work our focus
is on intermediate inclination angles that are not too close to the singular
points. Despite that, in the main text our analysis includes higher modes,
which break the complete degeneracy (“regularizing” the Fisher matrix) (c.f.
Fig. 6).

Now we repeat the above computation, but for the memory signal. Using the 0 PN
waveform in Eq. (6) we find that (in the frequency domain) the GW memory at
the detector is

$\widetilde{\delta
h}=\frac{\kappa_{\mathrm{m}}F^{+}}{d_{\mathrm{L}}}\sin^{2}\iota(17+\cos^{2}\iota),$
(18)

with a factor $\kappa_{\mathrm{m}}$ independent of both the distance and the
inclination, and such that
$\kappa_{0}/\kappa_{\mathrm{m}}\sim\mathcal{O}(100)$. The sky- and
polarization-averaged Fisher matrix of the memory is

$\Gamma^{\mathrm{m}}_{ij}=\left(\frac{\rho_{\kappa_{\mathrm{m}}}}{d_{\mathrm{L}}}\right)^{2}\hat{\Gamma}^{\mathrm{m}}_{ij},$
(19)

where
$\rho_{\kappa_{\mathrm{m}}}^{2}=(\kappa_{\mathrm{m}}|\kappa_{\mathrm{m}})$ and
the matrix elements

$\displaystyle\hat{\Gamma}^{\mathrm{m}}_{\log{d_{\mathrm{L}}},\log{d_{\mathrm{L}}}}=\frac{\sin^{4}\iota}{2}(17+\cos^{2}\iota)^{2},$
$\displaystyle\hat{\Gamma}^{\mathrm{m}}_{\log{d_{\mathrm{L}}},\iota}=-\sin(2\iota)\sin^{2}\iota\,(8+\cos^{2}\iota)(17+\cos^{2}\iota),$
$\displaystyle\hat{\Gamma}^{\mathrm{m}}_{\iota,\iota}=2\sin^{2}(2\iota)(8+\cos^{2}\iota)^{2},$

depend only on the inclination. Because of the simple structure of the memory
signal (at leading order), its Fisher matrix is singular for all inclination
angles. This is not an issue, since this singularity is cured through the
inclusion of (subleading) higher modes of the memory. Contrarily to what
happens with the primary signal, here the distance and the inclination are
_positively_ (_negatively_) correlated for inclination angles
$0<\iota<\frac{\pi}{2}$ ($\frac{\pi}{2}<\iota<\pi$). This opposite behavior is
nicely illustrated by the orthogonality of the two confidence ellipses in Fig.
3. Note that although these results were derived for the 0 PN waveforms, this
picture still holds generically, since it relies mostly on the leading
dependence on the inclination.

The Fisher matrix for the total (primary + memory) waveform
$\Gamma^{\mathrm{tot}}$ includes additional cross-terms,

$\Gamma^{\mathrm{tot}}_{ij}=\Gamma^{0}_{ij}+\Gamma^{\mathrm{m}}_{ij}+(\partial_{i}\widetilde{\delta
h}|\partial_{j}\tilde{h}_{0})+(\partial_{i}\tilde{h}_{0}|\partial_{j}\widetilde{\delta
h}).$ (20)

Above we focused on $\Gamma^{0}_{ij}$ and $\Gamma^{\mathrm{m}}_{ij}$. This is
because we verified that, due to the rapid oscillations of the integrands in
the cross-terms, the individual Fisher matrices dominate with respect to
those.

##  Appendix C Numerical Fisher matrix

To calculate the Fisher matrix elements we need to numerically compute
derivatives of the waveform. We do so using a second-order finite differences,

$\frac{\partial\tilde{h}}{\partial\Theta_{i}}\approx\frac{\partial\tilde{h}(\Theta_{i}+\Delta\Theta_{i})-\partial\tilde{h}(\Theta_{i}-\Delta\Theta_{i})}{2\Delta\Theta_{i}},$
(21)

except for the luminosity distance, for which we have the exact result

$\frac{\partial\tilde{h}}{\partial\log d_{\mathrm{L}}}=-\tilde{h},$ (22)

since $\tilde{h}\propto 1/d_{\mathrm{L}}$ (keeping $M_{z}$ fixed). We checked
that our Fisher matrices are numerically stable in our region of interest in
the parameter space for the increments:
$\Delta\Theta_{M_{z}}=10^{-3}M_{\odot}$, $\Delta\Theta_{q}=10^{-4}$,
$\Delta\Theta_{\iota,\varphi_{\mathrm{c}}}=10^{-6}\,\mathrm{rad}$. Varying by
an order of magnitude the finite increments gives just a few per cent change
in the final matrix elements. We also checked that the Fisher matrices are
stable by computing the derivatives with a higher-order finite differences
method. We verified the reliability of our Fisher matrix inversion by
confirming that in all cases,

$\mathrm{max}(|\Gamma_{ij}\Gamma^{-1}_{ij}-I_{ij}|)<10^{-6},$ (23)

where $I_{ij}$ is the identity matrix.

## References

  * Einstein [1916] A. Einstein, Näherungsweise Integration der Feldgleichungen der Gravitation, Sitzber. Preuss. Akad. Wiss. , 688 (1916).
  * Abbott  _et al._ [2019a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1, Phys. Rev. D 100, 104036 (2019a), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D 103, 122002 (2021a), arXiv:2010.14529 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Tests of General Relativity with GWTC-3, arXiv e-prints (2021b), arXiv:2112.06861 [gr-qc] .
  * Abbott  _et al._ [2019b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019b), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021c), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2021d] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, arXiv e-prints (2021d), arXiv:2111.03606 [gr-qc] .
  * Berti  _et al._ [2015] E. Berti _et al._ , Testing General Relativity with Present and Future Astrophysical Observations, Class. Quant. Grav. 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
  * Barack  _et al._ [2019] L. Barack _et al._ , Black holes, gravitational waves and fundamental physics: a roadmap, Class. Quant. Grav. 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Barausse  _et al._ [2020a] E. Barausse _et al._ , Prospects for Fundamental Physics with LISA, Gen. Rel. Grav. 52, 81 (2020a), arXiv:2001.09793 [gr-qc] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ , Laser Interferometer Space Antenna, arXiv e-prints (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Seoane  _et al._ [2013] P. A. Seoane _et al._ (eLISA), The Gravitational Universe, arXiv e-prints (2013), arXiv:1305.5720 [astro-ph.CO] .
  * Baibhav  _et al._ [2021] V. Baibhav _et al._ , Probing the nature of black holes: Deep in the mHz gravitational-wave sky, Exper. Astron. 51, 1385 (2021), arXiv:1908.11390 [astro-ph.HE] .
  * Arun  _et al._ [2022] K. G. Arun _et al._ (LISA), New horizons for fundamental physics with LISA, Living Rev. Rel. 25, 4 (2022), arXiv:2205.01597 [gr-qc] .
  * Cutler and Flanagan [1994] C. Cutler and E. E. Flanagan, Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral wave form?, Phys. Rev. D 49, 2658 (1994), arXiv:gr-qc/9402014 .
  * Usman  _et al._ [2019] S. A. Usman, J. C. Mills, and S. Fairhurst, Constraining the Inclinations of Binary Mergers from Gravitational-wave Observations, Astrophys. J. 877, 82 (2019), arXiv:1809.10727 [gr-qc] .
  * Chassande-Mottin  _et al._ [2019] E. Chassande-Mottin, K. Leyde, S. Mastrogiovanni, and D. A. Steer, Gravitational wave observations, distance measurement uncertainties, and cosmology, Phys. Rev. D 100, 083514 (2019), arXiv:1906.02670 [astro-ph.CO] .
  * Schutz [1986] B. F. Schutz, Determining the Hubble Constant from Gravitational Wave Observations, Nature 323, 310 (1986).
  * Holz and Hughes [2005] D. E. Holz and S. A. Hughes, Using gravitational-wave standard sirens, Astrophys. J. 629, 15 (2005), arXiv:astro-ph/0504616 .
  * Nissanke  _et al._ [2010] S. Nissanke, D. E. Holz, S. A. Hughes, N. Dalal, and J. L. Sievers, Exploring short gamma-ray bursts as gravitational-wave standard sirens, Astrophys. J. 725, 496 (2010), arXiv:0904.1017 [astro-ph.CO] .
  * Tamanini  _et al._ [2016] N. Tamanini, C. Caprini, E. Barausse, A. Sesana, A. Klein, and A. Petiteau, Science with the space-based interferometer eLISA. III: Probing the expansion of the Universe using gravitational wave standard sirens, JCAP 04, 002, arXiv:1601.07112 [astro-ph.CO] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), A gravitational-wave standard siren measurement of the Hubble constant, Nature 551, 85 (2017), arXiv:1710.05835 [astro-ph.CO] .
  * Poisson and Will [2014] E. Poisson and C. M. Will, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge University Press, 2014).
  * Hotokezaka  _et al._ [2019] K. Hotokezaka, E. Nakar, O. Gottlieb, S. Nissanke, K. Masuda, G. Hallinan, K. P. Mooley, and A. T. Deller, A Hubble constant measurement from superluminal motion of the jet in GW170817, Nature Astron. 3, 940 (2019), arXiv:1806.10596 [astro-ph.CO] .
  * Chen  _et al._ [2019] H.-Y. Chen, S. Vitale, and R. Narayan, Viewing angle of binary neutron star mergers, Phys. Rev. X 9, 031028 (2019), arXiv:1807.05226 [astro-ph.HE] .
  * Lang and Hughes [2006] R. N. Lang and S. A. Hughes, Measuring coalescing massive binary black holes with gravitational waves: The Impact of spin-induced precession, Phys. Rev. D 74, 122001 (2006), [Erratum: Phys.Rev.D 75, 089902 (2007), Erratum: Phys.Rev.D 77, 109901 (2008)], arXiv:gr-qc/0608062 .
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * Lang  _et al._ [2011] R. N. Lang, S. A. Hughes, and N. J. Cornish, Measuring parameters of massive black hole binaries with partially aligned spins, Phys. Rev. D 84, 022002 (2011), arXiv:1101.3591 [gr-qc] .
  * Arun  _et al._ [2007] K. G. Arun, B. R. Iyer, B. S. Sathyaprakash, and S. Sinha, Higher harmonics increase LISA’s mass reach for supermassive black holes, Phys. Rev. D 75, 124002 (2007), arXiv:0704.1086 [gr-qc] .
  * Porter and Cornish [2008] E. K. Porter and N. J. Cornish, The Effect of Higher Harmonic Corrections on the Detection of massive black hole binaries with LISA, Phys. Rev. D 78, 064005 (2008), arXiv:0804.0332 [gr-qc] .
  * Trias and Sintes [2008] M. Trias and A. M. Sintes, LISA observations of supermassive black holes: Parameter estimation using full post-Newtonian inspiral waveforms, Phys. Rev. D 77, 024030 (2008), arXiv:0707.4434 [gr-qc] .
  * Klein  _et al._ [2009] A. Klein, P. Jetzer, and M. Sereno, Parameter estimation for coalescing massive binary black holes with LISA using the full 2PN gravitational waveform and spin-orbit precession, Phys. Rev. D 80, 064027 (2009), arXiv:0907.3318 [astro-ph.CO] .
  * London  _et al._ [2018] L. London, S. Khan, E. Fauchon-Jones, C. García, M. Hannam, S. Husa, X. Jiménez-Forteza, C. Kalaghatgi, F. Ohme, and F. Pannarale, First higher-multipole model of gravitational waves from spinning and coalescing black-hole binaries, Phys. Rev. Lett. 120, 161102 (2018), arXiv:1708.00404 [gr-qc] .
  * Borhanian  _et al._ [2020] S. Borhanian, A. Dhani, A. Gupta, K. G. Arun, and B. S. Sathyaprakash, Dark Sirens to Resolve the Hubble–Lemaître Tension, Astrophys. J. Lett. 905, L28 (2020), arXiv:2007.02883 [astro-ph.CO] .
  * Yang  _et al._ [2022a] T. Yang, R.-G. Cai, Z. Cao, and H. M. Lee, Eccentricity of Long Inspiraling Compact Binaries Sheds Light on Dark Sirens, Phys. Rev. Lett. 129, 191102 (2022a), arXiv:2202.08608 [gr-qc] .
  * Yang  _et al._ [2022b] T. Yang, R.-G. Cai, and H. M. Lee, Space-borne atom interferometric gravitational wave detections. Part III. Eccentricity on dark sirens, JCAP 10, 061, arXiv:2208.10998 [gr-qc] .
  * Yang  _et al._ [2022c] T. Yang, R.-G. Cai, Z. Cao, and H. M. Lee, Parameter Estimation of Eccentric Gravitational Waves with Decihertz Observatory and Its Cosmological Implications, (2022c), arXiv:2212.11131 [gr-qc] .
  * Xie  _et al._ [2022] Y. Xie, D. Chatterjee, G. Holder, D. E. Holz, S. Perkins, K. Yagi, and N. Yunes, Breaking Bad Degeneracies with Love: Improving gravitational-wave measurements through universal relations, (2022), arXiv:2210.09386 [gr-qc] .
  * Hoy  _et al._ [2022] C. Hoy, C. Mills, and S. Fairhurst, Evidence for subdominant multipole moments and precession in merging black-hole-binaries from GWTC-2.1, Phys. Rev. D 106, 023019 (2022), arXiv:2111.10455 [gr-qc] .
  * Hannam  _et al._ [2021] M. Hannam, C. Hoy, J. E. Thompson, S. Fairhurst, and V. Raymond (LIGO/Virgo), Measurement of general-relativistic precession in a black-hole binary, (2021), arXiv:2112.11300 [gr-qc] .
  * Payne  _et al._ [2022] E. Payne, S. Hourihane, J. Golomb, R. Udall, D. Davis, and K. Chatziioannou, The curious case of GW200129: interplay between spin-precession inference and data-quality issues, (2022), arXiv:2206.11932 [gr-qc] .
  * Ng  _et al._ [2022] K. K. Y. Ng _et al._ , Measuring properties of primordial black hole mergers at cosmological distances: effect of higher order modes in gravitational waves, (2022), arXiv:2210.03132 [astro-ph.CO] .
  * Christodoulou [1991] D. Christodoulou, Nonlinear nature of gravitation and gravitational wave experiments, Phys. Rev. Lett. 67, 1486 (1991).
  * Payne [1983] P. N. Payne, SMARR’S ZERO FREQUENCY LIMIT CALCULATION, Phys. Rev. D 28, 1894 (1983).
  * Wiseman and Will [1991] A. G. Wiseman and C. M. Will, Christodoulou’s nonlinear gravitational wave memory: Evaluation in the quadrupole approximation, Phys. Rev. D 44, R2945 (1991).
  * Blanchet and Damour [1992] L. Blanchet and T. Damour, Hereditary effects in gravitational radiation, Phys. Rev. D 46, 4304 (1992).
  * Favata [2009a] M. Favata, Post-Newtonian corrections to the gravitational-wave memory for quasi-circular, inspiralling compact binaries, Phys. Rev. D 80, 024002 (2009a), arXiv:0812.0069 [gr-qc] .
  * Blanchet  _et al._ [2008] L. Blanchet, G. Faye, B. R. Iyer, and S. Sinha, The third post-Newtonian gravitational wave polarisations and associated spherical harmonic modes for inspiralling compact binaries in quasi-circular orbits, Class. Quant. Grav. 25, 165003 (2008), [Erratum: Class.Quant.Grav. 29, 239501 (2012)], arXiv:0802.1249 [gr-qc] .
  * Favata [2010] M. Favata, The gravitational-wave memory effect, Class. Quantum Gravity 27, 084036 (2010), arXiv:1003.3486 [gr-qc] .
  * Zel’dovich and Polnarev [1974] Y. B. Zel’dovich and A. G. Polnarev, Radiation of gravitational waves by a cluster of superdense stars, Sov. Astron. 18, 17 (1974).
  * Braginsky and Grishchuk [1985] V. B. Braginsky and L. P. Grishchuk, Kinematic Resonance and Memory Effect in Free Mass Gravitational Antennas, Sov. Phys. JETP 62, 427 (1985).
  * Braginsky and Thorne [1987] V. B. Braginsky and K. S. Thorne, Gravitational-wave bursts with memory and experimental prospects, Nature 327, 123 (1987).
  * Aurrekoetxea  _et al._ [2020] J. C. Aurrekoetxea, T. Helfer, and E. A. Lim, Coherent Gravitational Waveforms and Memory from Cosmic String Loops, Class. Quant. Grav. 37, 204001 (2020), arXiv:2002.05177 [gr-qc] .
  * Jenkins and Sakellariadou [2021] A. C. Jenkins and M. Sakellariadou, Nonlinear gravitational-wave memory from cusps and kinks on cosmic strings, Class. Quant. Grav. 38, 165004 (2021), arXiv:2102.12487 [gr-qc] .
  * Barausse  _et al._ [2012] E. Barausse, V. Morozova, and L. Rezzolla, On the mass radiated by coalescing black-hole binaries, Astrophys. J. 758, 63 (2012), [Erratum: Astrophys.J. 786, 76 (2014)], arXiv:1206.3803 [gr-qc] .
  * Favata [2011] M. Favata, The Gravitational-wave memory from eccentric binaries, Phys. Rev. D 84, 124013 (2011), arXiv:1108.3121 [gr-qc] .
  * Favata [2009b] M. Favata, Nonlinear gravitational-wave memory from binary black hole mergers, Astrophys. J. Lett. 696, L159 (2009b), arXiv:0902.3660 [astro-ph.SR] .
  * Hübner  _et al._ [2020] M. Hübner, C. Talbot, P. D. Lasky, and E. Thrane, Measuring gravitational-wave memory in the first LIGO/Virgo gravitational-wave transient catalog, Phys. Rev. D 101, 023011 (2020), arXiv:1911.12496 [astro-ph.HE] .
  * Ebersold and Tiwari [2020] M. Ebersold and S. Tiwari, Search for nonlinear memory from subsolar mass compact binary mergers, Phys. Rev. D 101, 104041 (2020), arXiv:2005.03306 [gr-qc] .
  * Zhao  _et al._ [2021] Z.-C. Zhao, X. Liu, Z. Cao, and X. He, Gravitational wave memory of the binary black hole events in GWTC-2, Phys. Rev. D 104, 064056 (2021), arXiv:2111.13882 [gr-qc] .
  * Hübner  _et al._ [2021] M. Hübner, P. Lasky, and E. Thrane, Memory remains undetected: Updates from the second LIGO/Virgo gravitational-wave transient catalog, Phys. Rev. D 104, 023004 (2021), arXiv:2105.02879 [gr-qc] .
  * Lasky  _et al._ [2016] P. D. Lasky, E. Thrane, Y. Levin, J. Blackman, and Y. Chen, Detecting gravitational-wave memory with LIGO: implications of GW150914, Phys. Rev. Lett. 117, 061102 (2016), arXiv:1605.01415 [astro-ph.HE] .
  * Yang and Martynov [2018] H. Yang and D. Martynov, Testing Gravitational Memory Generation with Compact Binary Mergers, Phys. Rev. Lett. 121, 071102 (2018), arXiv:1803.02429 [gr-qc] .
  * Boersma  _et al._ [2020] O. M. Boersma, D. A. Nichols, and P. Schmidt, Forecasts for detecting the gravitational-wave memory effect with Advanced LIGO and Virgo, Phys. Rev. D 101, 083026 (2020), arXiv:2002.01821 [astro-ph.HE] .
  * Grant and Nichols [2022] A. M. Grant and D. A. Nichols, Outlook for detecting the gravitational wave displacement and spin memory effects with current and future gravitational wave detectors, (2022), arXiv:2210.16266 [gr-qc] .
  * Buikema  _et al._ [2020] A. Buikema _et al._ , Sensitivity and performance of the Advanced LIGO detectors in the third observing run, Phys. Rev. D 102, 062003 (2020), arXiv:2008.01301 [astro-ph.IM] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The Einstein Telescope: A third-generation gravitational wave observatory, Class. Quant. Grav. 27, 194002 (2010).
  * Reitze  _et al._ [2019] D. Reitze _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Johnson  _et al._ [2019] A. D. Johnson, S. J. Kapadia, A. Osborne, A. Hixon, and D. Kennefick, Prospects of detecting the nonlinear gravitational wave memory, Phys. Rev. D 99, 044045 (2019), arXiv:1810.09563 [gr-qc] .
  * Islo  _et al._ [2019] K. Islo, J. Simon, S. Burke-Spolaor, and X. Siemens, Prospects for Memory Detection with Low-Frequency Gravitational Wave Detectors, (2019), arXiv:1906.11936 [astro-ph.HE] .
  * Islam  _et al._ [2021] T. Islam, S. E. Field, G. Khanna, and N. Warburton, Survey of gravitational wave memory in intermediate mass ratio binaries, (2021), arXiv:2109.00754 [gr-qc] .
  * Sun  _et al._ [2022] S. Sun, C. Shi, J.-d. Zhang, and J. Mei, Detecting the gravitational wave memory effect with TianQin, (2022), arXiv:2207.13009 [gr-qc] .
  * Seto [2009] N. Seto, Search for Memory and Inspiral Gravitational Waves from Super-Massive Binary Black Holes with Pulsar Timing Arrays, Mon. Not. Roy. Astron. Soc. 400, L38 (2009), arXiv:0909.1379 [astro-ph.CO] .
  * Cordes and Jenet [2012] J. M. Cordes and F. A. Jenet, Detecting gravitational wave memory with pulsar timing, Astrophys. J. 752, 54 (2012).
  * Wang  _et al._ [2015] J. B. Wang _et al._ , Searching for gravitational wave memory bursts with the Parkes Pulsar Timing Array, Mon. Not. Roy. Astron. Soc. 446, 1657 (2015), arXiv:1410.3323 [astro-ph.GA] .
  * Aggarwal  _et al._ [2020] K. Aggarwal _et al._ (NANOGrav), The NANOGrav 11 yr Data Set: Limits on Gravitational Wave Memory, Astrophys. J. 889, 38 (2020), arXiv:1911.08488 [astro-ph.HE] .
  * Seoane  _et al._ [2022] P. A. Seoane _et al._ , The effect of mission duration on LISA science objectives, Gen. Rel. Grav. 54, 3 (2022), arXiv:2107.09665 [astro-ph.IM] .
  * Barausse and Lapi [2020] E. Barausse and A. Lapi, Massive Black Hole Mergers, (2020), arXiv:2011.01994 [astro-ph.GA] .
  * Barausse  _et al._ [2020b] E. Barausse, I. Dvorkin, M. Tremmel, M. Volonteri, and M. Bonetti, Massive Black Hole Merger Rates: The Effect of Kiloparsec Separation Wandering and Supernova Feedback, Astrophys. J. 904, 16 (2020b), arXiv:2006.03065 [astro-ph.GA] .
  * Barausse [2012] E. Barausse, The evolution of massive black holes and their spins in their galactic hosts, Mon. Not. Roy. Astron. Soc. 423, 2533 (2012), arXiv:1201.5888 [astro-ph.CO] .
  * Sesana  _et al._ [2014] A. Sesana, E. Barausse, M. Dotti, and E. M. Rossi, Linking the spin evolution of massive black holes to galaxy kinematics, Astrophys. J. 794, 104 (2014), arXiv:1402.7088 [astro-ph.CO] .
  * Antonini  _et al._ [2015] F. Antonini, E. Barausse, and J. Silk, The Coevolution of Nuclear Star Clusters, Massive Black Holes, and their Host Galaxies, Astrophys. J. 812, 72 (2015), arXiv:1506.02050 [astro-ph.GA] .
  * Dey  _et al._ [2021] K. Dey, N. Karnesis, A. Toubiana, E. Barausse, N. Korsakova, Q. Baghi, and S. Basak, Effect of data gaps on the detectability and parameter estimation of massive black hole binaries with LISA, Phys. Rev. D 104, 044035 (2021), arXiv:2104.12646 [gr-qc] .
  * Pollney and Reisswig [2011] D. Pollney and C. Reisswig, Gravitational memory in binary black hole mergers, Astrophys. J. Lett. 732, L13 (2011), arXiv:1004.4209 [gr-qc] .
  * Mitman  _et al._ [2020] K. Mitman, J. Moxon, M. A. Scheel, S. A. Teukolsky, M. Boyle, N. Deppe, L. E. Kidder, and W. Throwe, Computation of displacement and spin gravitational memory in numerical relativity, Phys. Rev. D 102, 104007 (2020), arXiv:2007.11562 [gr-qc] .
  * Ashtekar  _et al._ [2020] A. Ashtekar, T. De Lorenzo, and N. Khera, Compact binary coalescences: Constraints on waveforms, Gen. Rel. Grav. 52, 107 (2020), arXiv:1906.00913 [gr-qc] .
  * Compère  _et al._ [2020] G. Compère, R. Oliveri, and A. Seraj, The Poincaré and BMS flux-balance laws with application to binary systems, JHEP 10, 116, arXiv:1912.03164 [gr-qc] .
  * Khera  _et al._ [2021] N. Khera, B. Krishnan, A. Ashtekar, and T. De Lorenzo, Inferring the gravitational wave memory for binary coalescence events, Phys. Rev. D 103, 044012 (2021), arXiv:2009.06351 [gr-qc] .
  * Mitman  _et al._ [2021] K. Mitman _et al._ , Adding gravitational memory to waveform catalogs using BMS balance laws, Phys. Rev. D 103, 024031 (2021), arXiv:2011.01309 [gr-qc] .
  * Thorne [1992] K. S. Thorne, Gravitational-wave bursts with memory: The Christodoulou effect, Phys. Rev. D 45, 520 (1992).
  * Kidder [2008] L. E. Kidder, Using full information when computing modes of post-Newtonian waveforms from inspiralling compact binaries in circular orbit, Phys. Rev. D 77, 044016 (2008), arXiv:0710.0614 [gr-qc] .
  * Varma  _et al._ [2019a] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E. Kidder, and H. P. Pfeiffer, Surrogate model of hybridized numerical relativity binary black hole waveforms, Phys. Rev. D 99, 064045 (2019a), arXiv:1812.07865 [gr-qc] .
  * Talbot  _et al._ [2018] C. Talbot, E. Thrane, P. D. Lasky, and F. Lin, Gravitational-wave memory: waveforms and phenomenology, Phys. Rev. D 98, 064031 (2018), arXiv:1807.00990 [astro-ph.HE] .
  * Tolish and Wald [2016] A. Tolish and R. M. Wald, Cosmological memory effect, Phys. Rev. D 94, 044009 (2016), arXiv:1606.04894 [gr-qc] .
  * Bieri  _et al._ [2017] L. Bieri, D. Garfinkle, and N. Yunes, Gravitational wave memory in $\Lambda$CDM cosmology, Class. Quant. Grav. 34, 215002 (2017), arXiv:1706.02009 [gr-qc] .
  * Jokela  _et al._ [2023] N. Jokela, K. Kajantie, and M. Sarkkinen, Gravitational wave memory in conformally flat spacetimes, (2023), arXiv:2301.07680 [gr-qc] .
  * Stein [2019] L. C. Stein, qnm: A Python package for calculating Kerr quasinormal modes, separation constants, and spherical-spheroidal mixing coefficients, J. Open Source Softw. 4, 1683 (2019), arXiv:1908.10377 [gr-qc] .
  * Varma  _et al._ [2018] V. Varma, L. C. Stein, and D. Gerosa, vijayvarma392/surfinbh: Surrogate final bh properties 10.5281/zenodo.1435832 (2018).
  * Varma  _et al._ [2019b] V. Varma, D. Gerosa, L. C. Stein, F. Hébert, and H. Zhang, High-accuracy mass, spin, and recoil predictions of generic black-hole merger remnants, Phys. Rev. Lett. 122, 011101 (2019b), arXiv:1809.09125 [gr-qc] .
  * Buonanno  _et al._ [2007] A. Buonanno, G. B. Cook, and F. Pretorius, Inspiral, merger and ring-down of equal-mass black-hole binaries, Phys. Rev. D 75, 124018 (2007), arXiv:gr-qc/0610122 .
  * Berti  _et al._ [2006] E. Berti, V. Cardoso, and C. M. Will, On gravitational-wave spectroscopy of massive black holes with the space interferometer LISA, Phys. Rev. D 73, 064030 (2006), arXiv:gr-qc/0512160 .
  * Reisswig  _et al._ [2009] C. Reisswig, S. Husa, L. Rezzolla, E. N. Dorband, D. Pollney, and J. Seiler, Gravitational-wave detectability of equal-mass black-hole binaries with aligned spins, Phys. Rev. D 80, 124026 (2009), arXiv:0907.0462 [gr-qc] .
  * Liu  _et al._ [2021] X. Liu, X. He, and Z. Cao, Accurate calculation of gravitational wave memory, Phys. Rev. D 103, 043005 (2021).
  * Finn [1992] L. S. Finn, Detection, measurement and gravitational radiation, Phys. Rev. D 46, 5236 (1992), arXiv:gr-qc/9209010 .
  * Poisson and Will [1995] E. Poisson and C. M. Will, Gravitational waves from inspiraling compact binaries: Parameter estimation using second postNewtonian wave forms, Phys. Rev. D 52, 848 (1995), arXiv:gr-qc/9502040 .
  * Cutler [1998] C. Cutler, Angular resolution of the LISA gravitational wave detector, Phys. Rev. D 57, 7089 (1998), arXiv:gr-qc/9703068 .
  * Vallisneri [2008] M. Vallisneri, Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects, Phys. Rev. D 77, 042001 (2008), arXiv:gr-qc/0703086 .
  * Berti  _et al._ [2005] E. Berti, A. Buonanno, and C. M. Will, Estimating spinning binary parameters and testing alternative theories of gravity with LISA, Phys. Rev. D 71, 084025 (2005), arXiv:gr-qc/0411129 .
  * Robson  _et al._ [2019] T. Robson, N. J. Cornish, and C. Liu, The construction and use of LISA sensitivity curves, Class. Quant. Grav. 36, 105011 (2019), arXiv:1803.01944 [astro-ph.HE] .
  * Harris  _et al._ [2020] C. R. Harris _et al._ , Array programming with NumPy, Nature 585, 357 (2020), arXiv:2006.10256 [cs.MS] .
  * Hirata  _et al._ [2010] C. M. Hirata, D. E. Holz, and C. Cutler, Reducing the weak lensing noise for the gravitational wave Hubble diagram using the non-Gaussianity of the magnification distribution, Phys. Rev. D 81, 124046 (2010), arXiv:1004.3988 [astro-ph.CO] .
  * Calderón Bustillo  _et al._ [2016] J. Calderón Bustillo, S. Husa, A. M. Sintes, and M. Pürrer, Impact of gravitational radiation higher order modes on single aligned-spin gravitational wave searches for binary black holes, Phys. Rev. D 93, 084019 (2016), arXiv:1511.02060 [gr-qc] .
  * Marsat  _et al._ [2021] S. Marsat, J. G. Baker, and T. Dal Canton, Exploring the Bayesian parameter estimation of binary black holes with LISA, Phys. Rev. D 103, 083011 (2021), arXiv:2003.00357 [gr-qc] .
  * Graff  _et al._ [2015] P. B. Graff, A. Buonanno, and B. S. Sathyaprakash, Missing Link: Bayesian detection and measurement of intermediate-mass black-hole binaries, Phys. Rev. D 92, 022002 (2015), arXiv:1504.04766 [gr-qc] .
  * Payne  _et al._ [2019] E. Payne, C. Talbot, and E. Thrane, Higher order gravitational-wave modes with likelihood reweighting, Phys. Rev. D 100, 123017 (2019), arXiv:1905.05477 [astro-ph.IM] .
  * Varma  _et al._ [2014] V. Varma, P. Ajith, S. Husa, J. C. Bustillo, M. Hannam, and M. Pürrer, Gravitational-wave observations of binary black holes: Effect of nonquadrupole modes, Phys. Rev. D 90, 124004 (2014), arXiv:1409.2349 [gr-qc] .
  * Abbott  _et al._ [2020] R. Abbott _et al._ (LIGO Scientific, Virgo), GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses, Phys. Rev. D 102, 043015 (2020), arXiv:2004.08342 [astro-ph.HE] .
  * London  _et al._ [2018] L. London, S. Khan, E. Fauchon-Jones, C. García, M. Hannam, S. Husa, X. Jiménez-Forteza, C. Kalaghatgi, F. Ohme, and F. Pannarale, First Higher-Multipole Model of Gravitational Waves from Spinning and Coalescing Black-Hole Binaries, Phys. Rev. Lett. 120, 161102 (2018), arXiv:1708.00404 [gr-qc] .
  * Hu and Veitch [2022] Q. Hu and J. Veitch, Assessing the model waveform accuracy of gravitational waves, Phys. Rev. D 106, 044042 (2022), arXiv:2205.08448 [gr-qc] .
  * Lindblom  _et al._ [2008] L. Lindblom, B. J. Owen, and D. A. Brown, Model Waveform Accuracy Standards for Gravitational Wave Data Analysis, Phys. Rev. D 78, 124020 (2008), arXiv:0809.3844 [gr-qc] .
  * Klein  _et al._ [2016] A. Klein _et al._ , Science with the space-based interferometer eLISA: Supermassive black hole binaries, Phys. Rev. D 93, 024003 (2016), arXiv:1511.05581 [gr-qc] .
  * Maggiore  _et al._ [2020] M. Maggiore _et al._ , Science Case for the Einstein Telescope, JCAP 03, 050, arXiv:1912.02622 [astro-ph.CO] .
  * Hild  _et al._ [2011] S. Hild _et al._ , Sensitivity Studies for Third-Generation Gravitational Wave Observatories, Class. Quant. Grav. 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Littenberg and Cornish [2023] T. B. Littenberg and N. J. Cornish, Prototype Global Analysis of LISA Data with Multiple Source Types, (2023), arXiv:2301.03673 [gr-qc] .
  * Palmese  _et al._ [2020] A. Palmese _et al._ , A statistical standard siren measurement of the hubble constant from the LIGO/virgo gravitational wave compact object merger GW190814 and dark energy survey galaxies, The Astrophysical Journal 900, L33 (2020).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing in Science & Engineering 9, 90 (2007).
  * Virtanen  _et al._ [2020] P. Virtanen _et al._ , SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Field  _et al._ [2014] S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, and M. Tiglio, Fast prediction and evaluation of gravitational waveforms using surrogate models, Phys. Rev. X 4, 031006 (2014), arXiv:1308.3565 [gr-qc] .
  * Richardson  _et al._ [2022] C. J. Richardson, M. Zanolin, H. Andresen, M. J. Szczepańczyk, K. Gill, and A. Wongwathanarat, Modeling core-collapse supernovae gravitational-wave memory in laser interferometric data, Phys. Rev. D 105, 103008 (2022), arXiv:2109.01582 [astro-ph.HE] .