# Measuring Supermassive Black Hole Properties via Gravitational Radiation
from Eccentrically Orbiting Stellar Mass Black Hole Binaries

Andrew Laeuger  0000-0002-8212-6496<EMAIL_ADDRESS>California Institute
of Technology, Pasadena, CA 91126  Brian Seymour  0000-0002-7865-1052
California Institute of Technology, Pasadena, CA 91126  Yanbei Chen
0000-0002-9730-9463 California Institute of Technology, Pasadena, CA 91126
Hang Yu  0000-0002-6011-6190 Montana State University, Bozeman, MT 59717

(September 2023)

###### Abstract

There may exist stellar-mass binary black holes (BBH) which merge while
orbiting nearby a supermassive black hole (SMBH). In such a triple system, the
SMBH will modulate the gravitational waveform of the BBH through orbital
Doppler shift and de Sitter precession of the angular momentum. Future space-
based GW observatories focused on the milli– and decihertz band will be
uniquely poised to observe these waveform modulations, as the GW frequency
from stellar-mass BBHs varies slowly in this band while modulation effects
accumulate. In this work, we apply the Fisher information matrix formalism to
estimate how well space-borne GW detectors can measure properties of BBH+SMBH
hierarchical triples using the GW from orbiting BBH. We extend previous work
by considering the more realistic case of an eccentric orbit around the SMBH,
and notably include the effects of orbital pericenter precession. We find that
for detector concepts such as LISA, B-DECIGO, and TianGO, we can extract the
SMBH mass and semimajor axis of the orbit with a fractional uncertainty below
the 0.1% level over a wide range of triple system parameters. Furthermore, we
find that the effects of pericenter precession and orbital eccentricity
significantly improve our ability to measure this system. We also find that
while LISA could measure these systems, the decihertz detector concepts
B-DECIGO and TianGO would enable better sensitivity to the triple’s
parameters.

##  I Introduction

Since the first detection of gravitational waves (GWs), GW astronomy by
ground-based detectors has cemented itself as an advantageous method for
studying binary systems of compact objects, the majority of which are binary
black holes (BBHs) [1, 2, 3]. Within the population of observed BBHs, there
are systems with progenitors whose masses exceed the predictions of stellar
evolution [4, 5, 6, 7]. One possible explanation of this detection could be
that the progenitors themselves were themselves products of previous mergers
[8, 9, 10, 11]. The deep potential wells created by supermassive black holes
(SMBHs) and their host galactic nuclei could trap the products of stellar mass
BBH mergers, making galactic nuclei ideal locations for generating many
repeated compact object mergers [12, 10, 11]. Numerical simulations of BBH
formation in galactic nuclei due to gas friction [13, 14] and dynamic capture
through gravitational interactions [15] suggest that the cosmological merger
rate of BBH near galactic nuclei could be of order $\sim$ a few
$\text{Gpc}^{-3}\text{yr}^{-1}$. Studying the properties of these repeated
merger systems and of the SMBHs which encourage their formation could open a
new window on understanding the dynamics of galactic nuclei and the processes
which drive galaxy evolution. The most recent analysis of the BBH population
in GWTC-3 is consistent with contributions from both isolated and AGN
formations [16], though more observations are needed.

In a hierarchical triple system consisting of a stellar-mass BBH orbiting an
SMBH, as depicted in Fig. 1, the presence of the SMBH would modulate the BBH
GW signal through many effects. For example, the velocity of the BBH in its
orbit will produce a Doppler shift in the waveform [17, 18, 19]. Allowing the
BBH to take an eccentric orbit around the SMBH introduces relativistic effects
such as pericenter precession as the outer orbital path approaches near the
SMBH [20]. Furthermore, the presence of the SMBH will cause the orbital
angular momentum of the inner binary $\hat{L}_{i}$ to experience de Sitter
precession about the orbital angular momentum of the outer binary
$\hat{L}_{o}$ [21]. This effect modulates the inclination angle of the BBH
angular momentum relative to an observatory in the Solar System. The Lidov-
Kozai and Lense-Thirring effects also play a role in the evolution of
hierarchical triples [22, 23].

By measuring the effects of Doppler shifts, pericenter precession, and de
Sitter precession on the stellar-mass BBH gravitational waveform, one can
measure the properties of this triple system, including the SMBH mass,
semimajor axis of the outer orbit, and various angles describing the system
geometry [24, 23, 19, 17, 18, 20, 22]. These effects accumulate substantially
over time scales roughly on the order of an orbital period, which for typical
BBH+SMBH triple systems can range from months to years. But because the
current ground detectors in LIGO/Virgo/KAGRA are most sensitive between 10 Hz
to a few kHz, which correspond to only the final seconds before merger for a
stellar-mass BBH, current GW observatories are not optimal for extracting
hierarchical triple system parameters through the influence of the SMBH on the
waveform [25, 26].

However, the coming decades could see the construction of a number of proposed
space-based detectors which would be sensitive to frequencies below $\sim$1-10
Hz. Building low frequency detectors in space is necessary due to technical
challenges from seismic noise [27, 28] and the need to create arms which are
large compared to the curvature of the Earth. The LISA [29], TianQin [30], and
Taiji [31, 32] detectors will target the millihertz GW band, while detector
concepts such as B-DECIGO [33, 34] and TianGO [35, 36] will focus on the
decihertz band. Since the instantaneous orbital decay timescale due to GW
emission during inspiral scales roughly with $\omega_{orb}^{-8/3}$ [37],
space-based low-frequency detectors could observe stellar mass BBH for much
longer times than ground detectors, making them more favorable for
measurements of SMBH-driven effects in the BBH waveform.

Measuring a SMBH with an orbiting binary’s GW would be be useful for studying
the environment at the center of galaxies. In a recent work by Yu and Chen, it
is shown that these proposed low-frequency GW observatories could feasibly
measure properties of interest to the few percent level over a wide range of
possible BBH+SMBH systems [24]. Current observational methods for measuring
properties of SMBHs and their local environments include tracking the orbital
dynamics of nearby test masses, like stars, and reverberation mapping of the
emission line fluxes from the accretion disk, if the SMBH is active [38].
Recent advances in observational technology and modeling active galactic
nuclei have enabled constraints of the masses of their central SMBHs to
roughly 10% precision [39, 40, 41, 42], though the results obtained by each
method do not always agree [43]. Adding a GW-based technique to this toolkit
could expand the set of observable SMBHs with well-constrained properties to
those which may have few electromagnetic radiation sources nearby [24] or
foster improvements in established electromagnetic techniques through
comparisons of joint measurements. Indeed, there has been significant progress
in understanding how space-based GW observatories may be able to measure
properties of SMBHs and the objects orbiting them through a variety of triple
system phenomena [44, 23, 45, 46].

The initial work of Yu and Chen assumes a circular Newtonian outer orbit in
the BBH+SMBH triple system [24]; however, it is expected that formation
channels for these systems, especially those which are dynamical in nature,
should produce a sizeable population of triples with eccentric outer orbits
[47]. In this work, we examine how adding a nonzero eccentricity to the outer
orbit affects parameter measurement uncertainties. We demonstrate that a
nonzero outer eccentricity can significantly improve these uncertainties
compared to the circular case, primarily through the inclusion of outer orbit
pericenter precession. In order to estimate parameter uncertainties, we rely
on the Fisher information matrix, a method which has been frequently used in
the past to gauge the measurability of compact binary parameters by ground-
based GW observatories [48]. In short, we find that uncertainties in triple
system parameters can consistently fall below the $0.1\%$ level, and that
these parameters are measured more precisely with larger $e_{o}$ and by
detectors targeting the decihertz band. We also find that the general trends
in parameter measurement are influenced almost entirely by pericenter and de
Sitter precession.

In Sec. II, we outline the mathematical description of the gravitational
waveform emitted from a BBH in a hierarchical triple and detected by a space-
borne observatory. In Sec. III, we outline the Fisher matrix calculation as
applied to parameter estimation and explain some simplifications we make to
the computation. In Sec. IV, we present the results of our Fisher matrix
computations, and in Sec. V, we offer conclusions and possible directions for
this work to proceed in the future. In this work, we use geometrized units
$G=c=1$.

##  II Mathematical Description of the SMBH+BBH Triple System

###  II.1 Geometry

We first describe the full geometry of the SMBH+BBH triple system with an
eccentric outer orbit. Table 1 below outlines the set of relevant parameters
used in calculating the waveform measured by a space-borne GW observatory. In
Fig. 1, the barred coordinates demarcate a Solar System centered coordinate
system, while the unbarred coordinates demarcate a coordinate system based on
the orientation of the observatory.

Figure 1: Top: Geometry of the SMBH+BBH triple system. Bottom, inset: View of the triple system normal to the plane of the outer orbit. The outer orbit angular momentum $\hat{L}_{o}$ points out of the page. See the discussion below and Table 1 for definition of all parameters. Figure dimensions are not an indication of true scale. $\boldsymbol{\theta^{a}}$ | Definition
---|---
$\log\mathcal{M}_{z}$ | Detector Frame Chirp Mass: $\mu^{3/5}(m_{1}+m_{2})^{2/5}$
$q$ | Mass Ratio $M_{2}/M_{1}$
$\log D_{L}$ | Luminosity Distance
$t_{c}$ | Coalescence Time
$\phi_{c}$ | Coalescence Phase
$\overline{\theta}_{S}$,$\overline{\phi}_{S}$ | Line of Sight of BBH+SMBH Triple
$\overline{\theta}_{J}$,$\overline{\phi}_{J}$ | Orientation of Total Angular Momentum $\boldsymbol{J}$
$\lambda_{L}$ | Angle Between $\boldsymbol{L}_{i}$ and $\boldsymbol{L}_{o}$
$\alpha_{0}$ | Initial Phase of $\boldsymbol{L}_{i}$ Around $\boldsymbol{L}_{o}$
$\log M_{3}$ | SMBH Mass
$\log a_{o}$ | Outer Orbit Semimajor Axis
$\gamma_{o}$ | Initial Outer Orbit Argument of Pericenter (See Note 1)
$e_{o}$ | Outer Orbit Eccentricity
$\varphi_{0}$ | Initial BBH Azimuthal Coordinate
Table 1: Relevant parameters in BBH+SMBH triple system for GW observed by
detectors. Bars over angles indicate the Solar System coordinate frame.

In order to compute the antenna response, we need to be able to convert from
the unbarred coordinates to the barred coordinates, which for a constellation-
preserving observatory such as LISA, is as follows [49]:

$\displaystyle\hat{x}=-\frac{1}{4}\sin(2\phi_{d})\hat{\overline{x}}+\frac{3+\cos(2\phi_{d})}{4}\hat{\overline{y}}+\frac{\sqrt{3}}{2}\sin(\phi_{d})\hat{\overline{z}}$
(1)
$\displaystyle\hat{y}=\frac{-3+\cos(2\phi_{d})}{4}\hat{\overline{x}}+\frac{1}{4}\sin(2\phi_{d})\hat{\overline{y}}-\frac{\sqrt{3}}{2}\cos(\phi_{d})\hat{\overline{z}}$
(2)
$\displaystyle\hat{z}=-\frac{\sqrt{3}}{2}\cos(\phi_{d})\hat{\overline{x}}-\frac{\sqrt{3}}{2}\sin(\phi_{d})\hat{\overline{y}}+\frac{1}{2}\hat{\overline{z}}.$
(3)

We note that even though B-DECIGO will posses a different detector geometry
than LISA during its orbit, we use the same configuration to simplify the
analysis. The sky location of the hierarchical triple is
$(\overline{\theta}_{S},\overline{\phi}_{S})$, which points along the vector
$\hat{N}$, and has a luminosity distance of $D_{L}$. The triple itself
consists of a BBH with black holes of masses $M_{1}$ and $M_{2}$, or
equivalently, a chirp mass of
$\mathcal{M}=\frac{(M_{1}M_{2})^{3/5}}{(M_{1}+M_{2})^{1/5}}$ and mass ratio of
$q=M_{2}/M_{1}$, and an SMBH of mass $M_{3}$. The shape of the BBH’s orbit
around the SMBH can be determined by the semimajor axis $a_{o}$, the
eccentricity $e_{o}$, the angle $\gamma_{o}$, analogous to the initial
Keplerian argument of pericenter111Of course, the outer orbit is not strictly
Keplerian. A rigorous definition of the instantaneous argument of pericenter
is subtle, though the picture of an elliptical orbital path with a pericenter
that rotates in space at the 1PN-accurate angular velocity of
$\frac{3M_{3}}{a_{o}(1-e_{o}^{2})}$ is appropriate as a rough approximation.
Within the mathematical framework of [50], $\gamma_{o}$ is implemented as a
simple arbitrary rotation of the orbital plane, as in Eq. (12). , and the
initial BBH azimuthal coordinate $\varphi_{o}$.

The unit vector of the angular momentum of the two lighter black holes in the
binary system is $\hat{L}_{i}$, and the unit vector of the angular momentum of
the binary’s orbit about the SMBH is $\hat{L}_{o}$. The opening angle
$\lambda_{L}$ is defined by

$\cos\lambda_{L}=\hat{L}_{o}\cdot\hat{L}_{i}.$ (4)

For $|\vec{L}_{o}|>>|\vec{L}_{i}|$ and neglecting long time scale orbital
effects as well as the spin of the SMBH (see Sec. II.4), the opening angle
stays constant in time, but the orientation of $\hat{L}_{i}$ traces a cone
around $\hat{L}_{o}$ due to de Sitter precession, with

$\frac{d\hat{L}_{i}}{dt}=\Omega_{dS}\hat{L}_{o}\times\hat{L}_{i}.$ (5)

Based on Eq. (9.200) of [51], we use the instantaneous de Sitter precession
frequency222Eq. (1) of the previous work [24] gave the orbit-averaged de
Sitter precession rate, which agrees with Eq. (6).

$\Omega_{dS}(t)=\frac{3}{2}\frac{M_{3}}{r(t)}\dot{\varphi}(t),$ (6)

where $r$ is the distance from the SMBH to the center of the BBH and and
$\varphi(t)$ is the azimuthal coordinate of the BBH in its orbit (as shown in
the inset of Fig. 1). The orbit-averaged precession rate is

$\langle\Omega_{dS}\rangle=\frac{3}{2}\frac{M_{3}}{a_{o}(1-e_{o}^{2})}\Omega_{o},$
(7)

where $\Omega_{o}\equiv\sqrt{M_{3}/a_{o}^{3}}$ is the Newtonian orbital
frequency. The phase of $\hat{L}_{i}$ in this cone, as shown in the inset of
Fig. 1, can be found by integrating the time-dependent de Sitter precession
rate:

$\alpha(t)=\alpha_{0}+\int_{t}^{t_{c}}\Omega_{dS}(t^{\prime})dt^{\prime},$ (8)

where $\alpha_{0}$ is the phase at the time of the binary coalescence $t_{c}$.

It is also useful to define the inclination angle $\iota_{J}$ of the outer
orbit angular momentum, given by

$\cos\iota_{J}=\hat{N}\cdot\hat{L}_{o}.$ (9)

###  II.2 BBH Orbit in Schwarzschild Spacetime

Despite the fact that there does not exist an analytic description of an
elliptical orbit in Schwarzschild spacetime, there are well-established
methods for computing Schwarzschild geodesics which can be applied to
numerically calculate the BBH orbital trajectory [50, 52, 53, 54, 55]. In
particular, we follow the procedure of [50]. Defining
$p=\frac{a_{o}}{M_{3}}(1-e_{o}^{2})$ for semimajor axis $a_{o}$ and $e_{o}$,
we find a minimum and maximum orbital radius

$r_{\text{min}}=\frac{pM_{3}}{1+e_{o}},\quad
r_{\text{max}}=\frac{pM_{3}}{1-e_{o}}$ (10)

Stable orbits only exist for $p>6+2e_{o}$ [50], and we will exclude unstable
systems from this analysis.

A relativistic anomaly $\chi$, which ranges from 0 to $2\pi$, is defined so
that

$r(\chi)=\frac{pM_{3}}{1+e_{o}\cos\chi},$ (11)

Furthermore, the azimuthal coordinate is given by

$\varphi(\chi)=2\Bigr{(}\frac{p}{p-6+2e_{o}}\Bigr{)}^{1/2}\Bigr{[}F\Bigr{(}\frac{\chi}{2}+\frac{\pi}{2},k^{2}\Bigr{)}-F\Bigr{(}\frac{\pi}{2},k^{2}\Bigr{)}\Bigr{]}+\gamma_{o},$
(12)

where $k^{2}=\frac{4e_{o}}{p-6+2e_{o}}$, $F$ is the incomplete elliptic
integral of the first kind, and $\gamma_{o}$ denotes the initial argument of
pericenter for the outer orbit (see Note 1).

The relationship between time and the relativistic anomaly is given by

$t(\chi)=p^{2}M_{3}(p-2-2e_{o})^{1/2}(p-2+2e_{o})^{1/2}\\\
\times\int_{0}^{\chi}d\chi^{\prime}\Bigr{\\{}(p-2-2e_{o}\cos\chi^{\prime})^{-1}(1+e_{o}\cos\chi^{\prime})^{-2}\\\
\times(p-6-2e_{o}\cos\chi^{\prime})^{-1/2}\Bigr{\\}}.$ (13)

In the end, the geodesic has a doubly periodic structure, and the radius has a
period of $r(\chi)$ has a period of $P_{r}=t(2\pi)$. During a time of $P_{r}$,
however, the azimuthal variable travels further than $2\pi$, which is the
relativistic pericenter precession. It is useful to define the shift in angle
over a radial period. This is equal to

$\Delta\varphi=4\Bigr{(}\frac{p}{p-6+2e_{o}}\Bigr{)}^{1/2}F(\pi/2,k^{2}).$
(14)

We note that this matches the 1PN GR result [56] for the amount of precession
during a radial period in the limit $p\gg 1$

$\Delta\varphi\approx 2\pi(1+3/p)=6\pi/p+2\pi\,.$ (15)

Defining the azimuthal frequency $\Omega_{\varphi}\equiv\Delta\varphi/P_{r}$,
it is shown that $\varphi(t)-\Omega_{\varphi}t$ is $P_{r}$–periodic [50]. We
note that $\varphi(t)$ itself is _not_ periodic – since the orbit precesses,
it takes $<P_{r}$ time for $\varphi$ to move through $2\pi$ radians. Even
though the precession angle over a full orbit remains constant, the time it
takes to move through the precession angle will depend on the BBH distance
from the SMBH (conserving angular momentum), so for an eccentric orbit, the
time to complete a full $2\pi$ in $\varphi$ will depend on the starting value
of $\varphi$ itself.

To find $r(t)$ and $\varphi(t)$ numerically over many full orbits, we
calculate the orbit over $\chi\in[0,2\pi]$ and utilize the periodicity of
$r(\chi)$ and $\varphi(t)-\Omega_{\varphi}t$. We furthermore choose some
$\chi_{0}\equiv\chi(t=0)$ so that $\varphi(\chi_{0})=\varphi_{0}$, where
$\varphi_{0}$ is the initial azimuthal coordinate of the BBH in the plane of
the outer orbit (see the bottom of Fig. 1). Furthermore, $\dot{r}(t)$ and
$\dot{\varphi}(t)$ can be calculated by application of the chain rule to the
expressions relating $r$, $\varphi$, and $t$ to $\chi$ above.

###  II.3 Waveform

We can now proceed to calculate the strain detected by the space-based
observatory, using the formalism of [21]. The overall measured signal is

$\tilde{h}(f)=\tilde{h}_{C}\sqrt{(A_{+}F_{+})^{2}+(A_{\times}F_{\times})^{2}}\\\
\times\exp\\{-i[\Phi_{P}+2\Phi_{T}+\Phi_{D}]\\},$ (16)

where $\tilde{h}_{C}$ is the carrier waveform of the BBH, $A_{+,\times}$ and
$F_{+,\times}$ are the polarization amplitude and antenna response,
respectively, and $\Phi_{P}$, $\Phi_{D}$, and $\Phi_{T}$ are the polarization,
Thomas, and Doppler phases. The carrier waveform in the frequency domain to
leading post-Newtonian (PN) order is [57]

$\tilde{h}_{C}(f)=\Bigr{(}\frac{5}{96}\Bigr{)}^{1/2}\frac{\mathcal{M}^{5/6}}{\pi^{2/3}D_{L}}f^{-7/6}\\\
\times\exp\\{i[2\pi
ft_{c}-\phi_{c}-\frac{\pi}{4}+\frac{3}{4}(8\pi\mathcal{M}f)^{-5/3}]\\},$ (17)

where $t_{c}$ and $\phi_{c}$ are the time and phase at coalescence. To the
leading PN order, the relationship between GW frequency and time is given by

$t(f)\approx t_{c}-\frac{5}{256\pi^{8/3}}\frac{1}{\mathcal{M}^{5/3}f^{8/3}}.$
(18)

The two polarizations of the strain, $h_{+}$ and $h_{\times}$, are modified by
the amplitude factors

$\displaystyle A_{+}=1+(\hat{L}_{i}\cdot\hat{N})^{2}$ (19) $\displaystyle
A_{\times}=-2\hat{L}_{i}\cdot\hat{N},$ (20)

and furthermore, the antenna responses for a 90-degree detector are

$F_{+}(\theta_{S},\phi_{S},\psi_{S})=\frac{1}{2}(1+\cos^{2}\theta_{S})\cos
2\phi_{S}\cos 2\psi_{S}\\\ -\cos\theta_{S}\sin 2\phi_{S}\sin 2\psi_{S}\,,$
(21)

$F_{\times}(\theta_{S},\phi_{S},\psi_{S})=\frac{1}{2}(1+\cos^{2}\theta_{S})\cos
2\phi_{S}\sin 2\psi_{S}\\\ +\cos\theta_{S}\sin 2\phi_{S}\cos 2\psi_{S},$ (22)

where

$\tan\psi_{S}(t)=\frac{\hat{L}_{i}\cdot\hat{z}-(\hat{L}_{i}\cdot\hat{N})(\hat{z}\cdot\hat{N})}{\hat{N}\cdot(\hat{L}_{i}\times\hat{z})}.$
(23)

Note the use of the detector-frame coordinates in Eqs. (21) and (22). For a
triangular detector such as LISA or B-DECIGO, the antenna pattern acquires a
factor of $\sqrt{3}/2$ and there are two effective detectors [58].

Let us now specify the phases in Eq. (16). Since the phases are slowly varying
functions of time, the stationary phase approximation is used to convert them
into frequency-dependent components via Eq. (18) – i.e., for some function
$g(t)$ appearing in the time-domain waveform $h(t)$, $g(f)\approx g(t(f))$
[59]. The polarization phase is given by

$\tan\Phi_{P}(t)=-\frac{A_{\times}(t)F_{\times}(t)}{A_{+}(t)F_{+}(t)}.$ (24)

The Thomas phase arises from the evolution of the principle $+$–polarization
axis [21], and thus the inner orbital phase of the two stellar mass BH in the
BBH, as the angular momentum $\hat{L}_{i}$ precesses. It is given by

$\Phi_{T}(t)=-\int_{t}^{t_{c}}dt\Bigr{[}\frac{\hat{L}_{i}\cdot\hat{N}}{1-(\hat{L}_{i}\cdot\hat{N})^{2}}\Bigr{]}(\hat{L}_{i}\times\hat{N})\cdot\frac{d\hat{L}_{i}}{dt}\,.$
(25)

The final phase term is the Doppler phase shift, the phase shift induced by
the changing distance between the detector and the GW source. There are two
contributions to this phase. The first is the contribution from the detector,
given at a particular time $t$ by

$\Phi_{D,\text{det}}=2\pi f\times(1\text{
AU})\sin\theta_{S}\cos(\phi_{det}-\phi_{S}).$ (26)

The other is from the source, which is modulated by the changing orbital
radius as well as the inclination of the outer orbit and the position of the
BBH in that orbit:

$\Phi_{D,\text{src}}=2\pi f\times r\sin\iota_{J}\sin\varphi.$ (27)

Gravitational lensing from the SMBH and its host galactic nucleus is neglected
in this waveform, though its effects on parameter estimation have been studied
in [44, 60].

###  II.4 Neglected Orbital Dynamics

A three body system is complicated, and exhibits some interesting
phenomenology. We will now discuss several additional well-known behaviors,
and why we neglect them. A useful benchmark for comparison is that the
characteristic frequency for de Sitter precession scales as

$\Omega_{\text{dS}}=\frac{1}{1.1\text{
yr}}\bigg{(}\frac{100}{a_{o}/M_{3}}\bigg{)}^{5/2}\bigg{(}\frac{10^{8}M_{\odot}}{M_{3}}\bigg{)}\bigg{(}\frac{1-0.3^{2}}{1-e^{2}}\bigg{)}\,.$
(28)

We consider the implications of non-zero BH spins on the orbital dynamics. The
precession of $\hat{L}_{o}$ around the spin of the SMBH $\hat{S}_{3}$ with
$S_{3}=\chi_{3}M_{3}^{2}$ has characteristic frequency [23]

$\Omega_{L_{o},S_{3}}=\frac{S_{3}(4+3(M_{1}+M_{2})/M_{3})}{2a_{o}^{3}(1-e_{o}^{2})^{3/2}}.$
(29)

If we consider the case $M_{3}\gg M_{1}+M_{2}$

$\frac{1}{t_{L_{o},S_{3}}}=\frac{1}{9.7\text{
yr}}\bigg{(}\frac{\chi_{3}}{0.7}\bigg{)}\bigg{(}\frac{100}{a_{o}/M_{3}}\bigg{)}^{3}\bigg{(}\frac{1-0.3^{2}}{1-e^{2}}\bigg{)}^{3/2}.$
(30)

Even for rapidly spinning SMBHs, this effect is about one order of magnitude
slower than de Sitter precession, so for now, we neglect it. It is worth
noting that each successive effect included in the waveform modulation
generally increases the amount of Fisher information. As such, we expect that
future inclusion of this effect will lead to further improved parameter
estimation uncertainties.

Lense-Thirring precession of $\hat{L}_{i}$ around $\hat{S}_{3}$ also
contributes to the orbital dynamics, with

$\Omega_{L_{i},S_{3}}=\frac{S_{3}}{2a_{o}^{3}(1-e_{o}^{2})^{3/2}}.$ (31)

This precession frequency is one-quarter of $\Omega_{L_{o},S_{3}}$, and thus,
since we treat $\Omega_{L_{o},S_{3}}$ as small in this work, we do the same
for $\Omega_{L_{i},S_{3}}$.

As in [24], we also neglect the precession of $\hat{L}_{i}$ around the spins
of the two stellar mass BH. The opening angle of this precession will be of
order $1^{\circ}$, much less than a typical value of $\lambda_{L}$ [61]. Also,
the effects of this spin-induced precession should be easily distinguishable
from the Doppler shift or de Sitter precession because the spin-induced
precession will occur over just days, rather than years, for GW frequencies in
the bands of space-based observatories.

We also consider Lidov-Kozai oscillations, the Newtonian tidal effect which
exchanges inner orbit eccentricity with inclination between $\hat{L}_{o}$ and
$\hat{L}_{i}$ [62]. These oscillations have a characteristic frequency of [23]

$\Omega_{\text{LK}}=\Omega_{i}\frac{M_{3}}{M_{1}+M_{2}}\Bigr{(}\frac{a_{i}}{a_{o}\sqrt{1-e_{o}^{2}}}\Bigr{)}^{3},$
(32)

where $\Omega_{i}=\sqrt{(M_{1}+M_{2})/a_{i}^{3}}$. The LK timescale is

$\frac{1}{t_{\text{LK}}}=\frac{1}{67\text{
yr}}\bigg{(}\frac{10^{8}M_{\odot}}{M_{3}}\bigg{)}^{2}\bigg{(}\frac{100}{a_{o}/M_{3}}\bigg{)}^{3}\\\
\times\bigg{(}\frac{1-0.3^{2}}{1-e^{2}}\bigg{)}^{3/2}\bigg{(}\frac{10^{-2}\text{
Hz}}{f}\bigg{)}.$ (33)

In our frequency band of interest, this effect occurs over much longer time
scales than the de Sitter precession, and since both de Sitter precession and
Lidov-Kozai oscillations modulate $L_{i}$, we neglect the slower of the two
processes.

We furthermore assume that the eccentricity of the inner binary $e_{i}$ is
zero. As explained in [24], the inner eccentricity does not affect any
component of the measured strain outside of the carrier waveform
$\tilde{h}_{C}(f)$, and thus should influence parameter estimation
uncertainties primarily through the SNR. Furthermore, the eccentric Kozai-
Lidov mechanism can drive periodic modulation of $e_{i}$ between moderate and
very high values. The GW signal frequency from a stellar-mass BBH can be
pushed into the sensitivity range of space-based observatories when the inner
eccentricity is high, so the eccentric Kozai-Lidov mechanism can produce
periodic high SNR bursts in these detectors, driving up the total SNR measured
for that particular binary [63, 64]. However, the time scale of this periodic
burst behavior scales roughly as [65] $\Omega_{o}^{-2}f(1-e_{o}^{2})^{3/2}$.
These effects therefore occur much more slowly than de Sitter and pericenter
precession, and thus are left for implementation into future analyses.

A higher $e_{i}$ also leads to faster merger times; however, high eccentricity
BBHs can still remain in the millihertz and decihertz frequency bands
throughout the entire observation period with just a larger initial separation
between the two stellar mass BHs. So, it is expected that even for $e_{i}$
approaching 1, such BBHs will offer long enough integration times to generate
a moderate SNR, and therefore the inner eccentricity should not significantly
alter the results of the simplified Fisher matrix analysis (see [24] for a
more detailed discussion).

##  III Parameter Estimation with the Fisher Information Matrix

In this analysis, we implement the Fisher information matrix method (as done
in [24]) as a simple estimator for how well properties of a BBH+SMBH triple
system can be measured. We make a number of well-supported assumptions to
reduce the complexity of the numerical methods used to estimate parameter
uncertainties.

###  III.1 Parameter Uncertainties from the Fisher Information Matrix

We first outline how the Fisher information matrix (from now on, Fisher
matrix) is used to estimate parameter measurement uncertainties. The elements
of the Fisher matrix are defined as

$\Gamma_{ab}\equiv\left(\frac{\partial\tilde{h}(f)}{\partial\theta_{a}}\Big{|}\frac{\partial\tilde{h}(f)}{\partial\theta_{b}}\right)\,,$
(34)

where

$\left(\tilde{g}\big{|}\tilde{h}\right)=4\,\text{Re}\int_{0}^{\infty}\frac{\tilde{g}^{\ast}(f)\tilde{h}(f)}{S_{n}(f)}df\,,$
(35)

$\tilde{h}$ is the frequency-domain waveform, $S_{n}(f)$ is the PSD of the
detector noise, and $\theta_{a}$ are the various parameters of the system. In
practice, we limit the frequency bounds of integration to $[f_{min},f_{max}]$,
where $f_{max}$ is at the upper edge of the detector sensitivity range and
$t(f_{max})-t(f_{min})=5$ years (via Eq. (18)) – see Sec. IV.

We note that we use a finite difference method to compute
$\partial\tilde{h}/\partial\theta_{a}$. To choose a finite parameter
difference $\Delta\theta_{a}$ from which to estimate
$\partial\tilde{h}/\partial\theta_{a}$, we minimize the quantity $\epsilon$,
analogous to waveform mismatch,

$\epsilon=1-\frac{(\partial_{[\Delta\theta_{a}]}\tilde{h}|\partial_{[4\Delta\theta_{a}]}\tilde{h})}{\sqrt{(\partial_{[\Delta\theta_{a}]}\tilde{h}|\partial_{[\Delta\theta_{a}]}\tilde{h})(\partial_{[4\Delta\theta_{a}]}\tilde{h}|\partial_{[4\Delta\theta_{a}]}\tilde{h})}},$
(36)

where

$\partial_{[\Delta\theta]}\tilde{h}=\frac{\tilde{h}(\theta+\Delta\theta)-\tilde{h}(\theta-\Delta\theta)}{2\Delta\theta}.$
(37)

Empirically choosing $\Delta\theta_{a}$ to make $\epsilon$ small gives us the
best accuracy in computing the numerical derivative, as $\epsilon$ begins to
increase once $\Delta\theta_{a}$ becomes so small that the changes in
$\tilde{h}$ are smaller than computer precision. The choice of
$4\Delta\theta_{a}$ to compare to $\Delta\theta_{a}$ is arbitrary.

The Fisher information matrix is related to the covariance matrix roughly by

$\Sigma_{ab}=[\Gamma^{-1}]_{ab}+\mathcal{O}(\rho^{-4}),$ (38)

where $\rho$ is the signal-to-noise ratio (SNR). So, in the limit of large
SNR, the covariance between two parameters $\Delta\theta_{i}\Delta\theta_{j}$
is approximately equal to the corresponding element of the inverse of the
Fisher information matrix. As such, the parameter estimation uncertainty is
given by $\Delta\theta_{i}=(\Sigma_{ii})^{0.5}$. If a network of GW detectors
were to observe the same system, the Fisher information matrix would scale as
the sum of the matrix elements for each detector, or

$\left(\Gamma_{ab}\right)^{\text{network}}=\sum_{\text{det}}\Gamma_{ab}^{\text{det}}.$
(39)

This also applies to a triangular observatory, wherein three arms compose two
interferometric detectors.

###  III.2 Reduced Fisher Matrix Dimensions

We can reduce the dimensions of the Fisher matrix by removing certain physical
parameters from the analysis. Doing so reduces the total computation time as
well as the condition number, leading to improved numerical accuracy in the
Fisher matrix inversion [48]. From the parameters listed in Table 1, our
Fisher matrices include the following 12 parameters:

$\theta_{a}=(\log
D_{L},\overline{\theta}_{S},\overline{\phi}_{S},\overline{\theta}_{J},\overline{\phi}_{J},\lambda_{L},\alpha_{0},\\\
\log M_{3},\log\Omega_{o},\gamma_{o},e_{o},\varphi_{0}).$ (40)

We can remove parameters which we expect will have strong priors obtained from
other GW measurements, or which contribute only weakly to the gravitational
waveform. For example, we assume that space-based detectors like LISA or
TianGO will act in conjunction with ground-based observatories, which are far
more sensitive to the chirp mass $\mathcal{M}$, the mass ratio $q$, and the
time and phase of coalescence $t_{c}$ and $\phi_{c}$ [36], and thus treat
these four parameters as perfectly known in our analysis. Removing the chirp
mass from the Fisher matrix also improves the numerical stability of our
analysis. Furthermore, we neglect the spins of the three black holes because
the precessional effects they induce accumulate much more slowly than the
outer orbital motion and de Sitter precession, as described in Sec. II.4.

##  IV Results and Discussion

We examine a BBH+SMBH triple system with fixed parameters
$M_{1}=M_{2}=50M_{\odot}$, $t_{c}=0$, $\phi_{c}=0$, $D_{L}=1$Gpc,
$(\overline{\theta}_{S},\overline{\phi}_{S})=(33^{\circ},147^{\circ})$,
$(\overline{\theta}_{J},\overline{\phi}_{J})=(75^{\circ},150^{\circ})$, and
$\lambda_{L}=45^{\circ}$. For B-DECIGO, TianGO, and LISA, we compute the
Fisher matrix where the integration is taken over a frequency window
corresponding to an observation time of five years and the highest frequency
is $f_{\text{max}}=$ 12 Hz – this roughly corresponds to a lowest frequency of
$f_{\text{min}}\sim 12$ mHz. In Fig. 2, we plot an example frequency-domain
waveform along with the B-DECIGO, TianGO, and LISA sensitivity curves used in
computing Fisher matrix elements.

Figure 2: An example waveform $\tilde{h}(f)$ with $M_{3}=10^{8}M_{\odot}$,
$a_{o}=100M_{3}$, and $e_{o}=0.3$, along with approximate sensitivity curves
for B-DECIGO, TianGO, and LISA used in the Fisher matrix calculations done in
this work. The red dashed curve gives the same waveform but with the effects
of de Sitter precession removed.

In Fig. 3, we plot the fractional uncertainty in the SMBH mass $M_{3}$,
measured by B-DECIGO, as we vary $M_{3}$ and $a_{o}$. The Fisher matrix breaks
down if $e_{o}$ is identically zero, so in order to facilitate comparisons to
the circular orbits used in [24], we use $e_{o}=0.001$. At each point, we
sample the covariance found with the Fisher matrix over combinations of the
three geometrical phases – that is, 6 choices of $\gamma_{o}$, $\varphi_{o}$
and $\alpha_{0}$, or 216 sets of $(\gamma_{o},\alpha_{0},\varphi_{o})$ – and
find the median.

The purple regions denote where the outer binary merges in less than the
proposed observation length of five years. We expect systems in this region to
be exceedingly rare, as there is only a short window for such systems to form
in order to be detected by B-DECIGO. We also shade out the region where the
outer orbital period $P_{\text{outer}}$ exceeds twice the observation
duration. In this region, the most dominant source of waveform modulation –
namely, the Doppler phase shift – is difficult to measure because the BBH only
passes through a small range of angles over the observation period.
Furthermore, when the Doppler phase shift varies slowly, remaining roughly
constant over the observation run, it becomes degenerate with $t_{c}$, which
itself can be changed by a simple redefinition of when $t=0$. So, in this
shaded region, our assumption that $t_{c}$ can be safely removed from the list
of parameters in the Fisher matrix does not hold well. Indeed, we encounter
problems with numerical instability when computing the Fisher matrix in this
region of the contour plots.

Figure 4 gives the same results, but using the LISA detector response and
noise curve instead of that of B-DECIGO. The contour plots using the TianGO
observatory have a similar structure to those using B-DECIGO, as the two
detectors have similar sensitivity curves. Across the majority of the
parameter space studied, the two sets of contours differ only in magnitude and
not in shape, so for the sake of brevity, they are omitted here.

Figure 3: Fractional uncertainty in $M_{3}$ as measured by B-DECIGO for three
different eccentricities $e_{o}=\\{0.001,0.3,0.6\\}$. At each point in the
contour plot, we take the median uncertainty over a set of combinations of
$(\gamma_{o},\alpha_{0},\phi_{o})$. The purple region corresponds to where the
outer binary merges in less time than the observation duration. We lightly
shade out the region with an outer orbital period greater than 10 years, where
the cumulative effect of the Doppler shift becomes small.  Figure 4: Same as
Fig. 3, but measured by LISA instead.

We note that the fractional uncertainty in the outer orbit semimajor axis
$\Delta a_{o}/a_{o}$ follows a similar contour structure to that of $\Delta
M_{3}/M_{3}$. For the outer orbit,

$3\frac{a_{o}^{3}}{M_{3}}\frac{\Delta
a_{o}}{a_{o}}\approx\frac{1}{\Omega_{o}^{2}}\frac{\Delta
M_{3}}{M_{3}}-2\frac{1}{\Omega_{o}^{2}}\frac{\Delta\Omega_{o}}{\Omega_{o}}.$
(41)

Our calculations determined that across the $(M_{3},a_{o}/M_{3})$ parameters
space, $\Delta\Omega_{o}/\Omega_{o}$ is much smaller in magnitude than $\Delta
M_{3}/M_{3}$, so

$\frac{\Delta a_{o}}{a_{o}}\approx\frac{1}{3}\frac{\Delta M_{3}}{M_{3}}.$ (42)

This result is verified in the structures of Figs. 5 and 6, and we observe
that both B-DECIGO and LISA have the potential to realize fractional
uncertainties in $M_{3}$ and $a_{o}$ significantly below the 0.1% level across
a wide range of parameters of the triple systems.

Figure 5: Uncertainty in $a_{o}$ as measured by B-DECIGO for three different
eccentricities $e_{o}=\\{0.001,0.3,0.6\\}$. The same sampling procedure as
used in Fig. 3 is applied here. Figure 6: Same as Fig. 5, but measured by LISA
instead.

To understand the structure of the contour plots, we examine the contour plot
in Fig. 7. For small $M_{3}$ and $a_{o}/M_{3}$, the shape of the contours are
roughly separated by lines of constant $a_{o}^{5}/M_{3}^{3}$. We correlate
these trends to evolving components of the waveform. First, the de Sitter
precession frequency is proportional to
$\Omega_{\text{dS}}\propto\sqrt{M_{3}^{3}/a_{o}^{5}}$. As discussed in App. B,
the Thomas phase and polarization phase scale as
$\Phi_{T}\sim\Omega_{\text{dS}}t$. Thus, measurement accuracy scales with the
number of de Sitter cycles within the five-year window. In this region of
parameter space, the modulations of de Sitter precession are the dominant
effect for how well we can measure $M_{3},a_{o}$.

For larger $M_{3}$ and a wide range of $a_{o}/M_{3}$, the shape of the
contours are roughly separated by lines of constant $a_{o}/M_{3}$. In this
region, the Doppler phase is the dominant term in the frequency domain
waveform phase. The Doppler phase magnitude features a degeneracy between
$a_{o}$ and $\sin\iota_{J}$ (with $\sin\iota_{J}$ being a function of the
angles $\overline{\theta}_{S}$, $\overline{\phi}_{S}$,
$\overline{\theta}_{J}$, and $\overline{\phi}_{J}$), as these quantities
appear in the magnitude only as the product
$a_{o}\sin\iota_{J}=M_{3}^{1/3}\Omega_{o}^{-2/3}\sin\iota_{J}$. This
degeneracy is broken by the inclusion of relativistic pericenter precession,
as this produces different periods in the radial and azimuthal motion of the
BBH in the outer orbit (cf. Sec. II.2). The inclusion of this precession
produces lines of constant $\Delta M_{3}/M_{3}$ that scale roughly with
$(a_{o}/M_{3})^{3/2}$. See App. B for more detailed discussion.

Studying Fig. 3, we see that for $e_{o}\approx 0$, these flat contours do not
appear, as for a circular orbit, pericenter precession is essentially
consistent with an increase in $\Omega_{o}$. The resulting contour plot shape
is similar to the results seen in Fig. 5 of [24], where $e_{o}$ is assumed to
be zero – over a wide range of the parameter space, de Sitter precession is
the dominant effect in determining $\Delta M_{3}/M_{3}$. However, once
$e_{o}>0$, pericenter precession, rather than de Sitter precession, becomes
the leading contribution to $\Delta M_{3}/M_{3}$ over a significant portion of
the parameter space. The importance of pericenter precession is further
emphasized by comparing the magnitudes of $\Delta M_{3}/M_{3}$ in our plots to
Fig. 5 of [24], which sets $e_{o}=0$ and therefore does not include pericenter
precession (though it does include all other effects used in this work). With
pericenter precession included, the parameter uncertainties across a wide
region of the overall parameter space can drop by multiple orders of
magnitude.

Figure 7: Contour plot for the fractional uncertainty in $M_{3}$ as measured
by B-DECIGO, taken from Fig. 3. Plotted on top of the contours are lines of
constant $a_{o}^{5/2}/M_{3}^{3/2}$ and $M_{3}^{1/2}/a_{o}^{1/2}$ to indicate
the structure of the contours.

We also estimate how well the eccentricity can be measured with B-DECIGO and
LISA as shown in Fig. 8 and Fig. 9. These results suggest that the
eccentricity can be constrained to high precision, with B-DECIGO able to
achieve a lower bound of $\Delta e_{o}\sim 10^{-6}-10^{-5}$ and LISA able to
achieve $\Delta e_{o}\sim 10^{-5}-10^{-4}$ across a substantial portion of the
parameter space where precession is detectable. Once again, we see the
importance of de Sitter precession in the measurability of this parameter – in
the portion of the parameter space where de Sitter precession is rapid,
equivalent estimation uncertainties match contours of equal de Sitter
precession period. Unlike the contour plots for $\Delta M_{3}/M_{3}$, the
shape of these contours is not heavily dictated by power laws related to
pericenter precession. Indeed, there are no degeneracies between $e_{o}$ and
other waveform parameters which are broken by pericenter precession.

Figure 8: Uncertainty in $e_{o}$ as measured by B-DECIGO for three different
eccentricities $e_{o}=\\{0.001,0.3,0.6\\}$. The same sampling procedure as
used in Fig. 3 is applied here. Figure 9: Same as Fig. 8, but measured by LISA
instead.

An important question is the impact of increasing outer orbit eccentricity on
the ability to measure parameters like $M_{3}$, $a_{o}$, and $e_{o}$ itself.
In Figs. 10 and 11, we consider B-DECIGO, LISA, and the TianGO concept and
three different combinations of $(M_{3},a_{o}/M_{3})$ across our chosen
parameter space. We study the effect of increasing eccentricity on the
estimation uncertainties in $M_{3}$ and $e_{o}$ (still averaging over initial
orbital angles) and find that increasing eccentricity can produce marginal
improvements in the measurement of $M_{3}$ and $e_{o}$ – a factor of $\sim$ a
few – though such improvement is not universal across $(M_{3},a_{o}/M_{3})$
parameter space.

Considering the arguments given in App. B, we see that the leading
contributions to the Fisher matrix elements come from the derivatives of
$\Phi_{D}$, $\Phi_{P}$, and $\Phi_{T}$. Noting that these phases evolve at
secular rates of $\Omega_{\text{dS}}$ (for $\Phi_{P}$ and $\Phi_{T}$) or
$\Omega_{\text{pericenter}}=\Omega_{o}\frac{3}{p}$ (for $\Phi_{D}$ –
specifically, this is the rate at which the degeneracy between $a_{o}$ and
$\sin\iota_{J}$ is broken), and recalling that these rates scale with
$(1-e_{o}^{2})^{-1}$, it follows that larger eccentricities produce more rapid
evolution, larger Fisher matrix entries, and ultimately smaller parameter
uncertainties.

The relative sensitivities between the three detectors are responsible for the
clear hierarchy in the parameter uncertainties they produce. For example, the
rates of precession and orbital velocity are sensitive to both $M_{3}$ and
$e_{o}$ but with different dependencies, so there exist degeneracies between
these two parameters. These degeneracies can be lifted by observing the system
over long periods of time so that these effects can accumulate, enabling
tighter constraints on their respective individual rates. Examining Fig. 2, we
see that LISA effectively measures the BBH signal over a smaller frequency
band than the other two detectors in the five years prior to merger. Since the
LISA sensitivity is poorer than the other two detectors in the frequencies
sampled in the five year observation run, the SNR of the waveform is reduced
and it becomes more difficult to extract the waveform modulations driven by
orbital and precessional effects over that period of time. Therefore, the
degeneracies are not as cleanly lifted in LISA measurements, especially when
these rates are slow (i.e., low $M_{3}$, high $a_{o}/M_{3}$), producing less
precise parameter estimates.

The primary effect of eccentricity then is to increase the strength of
waveform modulations by increasing the magnitude of the precessional effects
(pericenter, de Sitter); however, we see that for the LISA observatory, the
improvement in parameter estimation uncertainty with rising $e_{o}$ is not as
significant as in B-DECIGO and TianGO, and in some cases, a larger $e_{o}$
produces larger uncertainties. While increasing the eccentricity boosts the
orbit averaged rate of de Sitter and pericenter precession (Cf. Eqs. 7 and
15), the majority of this evolution occurs when the BBH is near the outer
orbit pericenter and the instantaneous precession rate is largest. So, for
systems with slow outer orbits (once again, low $M_{3}$ and high
$a_{o}/M_{3}$), an increasing eccentricity constrains the majority of the
waveform modulation effects to a shorter time window, as the BBH passes
through the region near the pericenter at a faster rate. The GW radiation from
the BBH then evolves through a smaller range of frequencies while the waveform
is significantly modulated.

Figure 10: The fractional uncertainties in $M_{3}$ obtainable by B-DECIGO
(blue), TianGO (orange), and LISA (green) as the eccentricity is varied. The
solid, dashed, and dotted lines correspond to different choices of
$(M_{3},a_{o}/M_{3})$.  Figure 11: Same as Fig. 10, but estimating the outer
eccentricity variance $\Delta e_{o}$.

##  V Conclusion and Future Directions

Using the Fisher information matrix, we have shown that future space-based GW
observatories may be able to precisely constrain the properties of BBH+SMBH
triple systems, like the SMBH mass and outer orbit semimajor axis and
eccentricity, through the GW signal observed from the BBH. We have
demonstrated that the rate of change of the Doppler phase shift and the de
Sitter precession rate are the dominant factors determining of the
measurability of triple system parameters and that an increasing outer orbit
eccentricity leads to improved measurement uncertainties through greater
Doppler phase shift modulation and faster de Sitter precession. We have also
shown that the planned LISA detector is capable of measuring these systems,
though decihertz detector concepts such as TianGO or B-DECIGO would possess a
competitive advantage over LISA in measuring such quantities.

There are some important limitations of the Fisher information method
implemented in this work. As described in [48], a high SNR is required for the
inverse Fisher matrix to give the covariance of the posterior probability
distribution for the true source parameters $\vec{\theta}_{0}$. While the SNR
we compute for our waveform is generally $\sim 40$ for TianGO, it is only
$\sim 4$ for LISA, suggesting that the true parameter estimation uncertainties
may be significantly different than those calculated here. However, the
inverse Fisher matrix is also a _lower bound_ for the uncertainty of an
unbiased estimator of $\vec{\theta}_{0}$ [48], so our results essentially
offer a best-case scenario for the parameter estimation precision obtainable
by future space-based observatories. A more thorough approach to this analysis
will implement a full Bayesian methodology.

We can further develop this work by inclusion of additional effects into the
waveform. One can implement the spin-precession effects that we chose to
neglect in II.4 due to their significantly slower time scales. Furthermore,
for triple systems with lower outer binary merger times (i.e., with $M_{3}$
and $a_{o}/M_{3}$ near the purple regions shown in the contour plots such as
Fig. 3), the semimajor axis and outer eccentricity can evolve significantly in
time due to radiation reaction [37]. Considering the frequency integral that
composes the Fisher matrix elements, we can include the effects of
gravitational redshift and Doppler frequency shift, which would require the
waveform and detector sensitivity to be evaluated at different frequencies in
the integrand. Also, the stationary phase approximation used in the frequency
domain waveform (outlined in App. A) may not hold well for highly eccentric
outer orbits, as the outer orbital angle varies quite rapidly near the
pericenter for such orbits.

In Ref. [44], it is discussed how gravitational lensing of GWs by the SMBH
combined with the de Sitter precession of $\boldsymbol{L}_{\mathrm{i}}$ can
further constrain the parameters of a triple system as estimated by a space-
based GW observatory, even in the case of a circular outer orbit. It would be
interesting to examine the combined effects of an eccentric outer orbit and
repeated GW lensing in parameter estimation problems.

Finally, measurements of the motion of a BBH through space through its
modulated waveform may prove useful for understanding phenomena besides
BBH+SMBH hierarchical triples. For example, measuring the evolving Doppler
shift and aberrations induced by the evolving position and velocity of an
isolated BBH might enable estimates of BBH kicks that occur shortly before
merger or improve the precision of estimates of the Hubble constant by further
constraining the redshifts of GW standard sirens [66, 67].

###### Acknowledgements.

This work was supported in part by the Caltech LIGO Summer Undergraduate
Research Fellowship and by the REU program of the NSF. B.S. acknowledges
support by the National Science Foundation Graduate Research Fellowship under
Grant No. DGE-1745301. H.Y. is supported by NSF PHY-2308415.

## References

  * Abbott  _et al._ [2019] 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 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2021b] 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, (2021b), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ [2020a] R. Abbott _et al._ (LIGO Scientific, Virgo), GW190521: A Binary Black Hole Merger with a Total Mass of $150M_{\odot}$, Phys. Rev. Lett. 125, 101102 (2020a), arXiv:2009.01075 [gr-qc] .
  * Woosley [2017] S. E. Woosley, Pulsational Pair-Instability Supernovae, Astrophys. J. 836, 244 (2017), arXiv:1608.08939 [astro-ph.HE] .
  * Belczynski  _et al._ [2016] K. Belczynski _et al._ , The Effect of Pair-Instability Mass Loss on Black Hole Mergers, Astron. Astrophys. 594, A97 (2016), arXiv:1607.03116 [astro-ph.HE] .
  * Spera and Mapelli [2017] M. Spera and M. Mapelli, Very massive stars, pair-instability supernovae and intermediate-mass black holes with the SEVN code, Mon. Not. Roy. Astron. Soc. 470, 4739 (2017), arXiv:1706.06109 [astro-ph.SR] .
  * Abbott  _et al._ [2020b] R. Abbott _et al._ (LIGO Scientific, Virgo), Properties and Astrophysical Implications of the 150 M⊙ Binary Black Hole Merger GW190521, Astrophys. J. Lett. 900, L13 (2020b), arXiv:2009.01190 [astro-ph.HE] .
  * Gerosa and Fishbach [2021] D. Gerosa and M. Fishbach, Hierarchical mergers of stellar-mass black holes and their gravitational-wave signatures, Nature Astron. 5, 749 (2021), arXiv:2105.03439 [astro-ph.HE] .
  * Antonini  _et al._ [2019] F. Antonini, M. Gieles, and A. Gualandris, Black hole growth through hierarchical black hole mergers in dense star clusters: implications for gravitational wave detections, Mon. Not. Roy. Astron. Soc. 486, 5008 (2019), arXiv:1811.03640 [astro-ph.HE] .
  * Martinez  _et al._ [2020] M. A. S. Martinez, G. Fragione, K. Kremer, S. Chatterjee, C. L. Rodriguez, J. Samsing, C. S. Ye, N. C. Weatherford, M. Zevin, S. Naoz, and F. A. Rasio, Black Hole Mergers from Hierarchical Triples in Dense Star Clusters, Astrophys. J. 903, 67 (2020), arXiv:2009.08468 [astro-ph.GA] .
  * Yang  _et al._ [2019] Y. Yang, I. Bartos, V. Gayathri, K. E. S. Ford, Z. Haiman, S. Klimenko, B. Kocsis, S. Márka, Z. Márka, B. McKernan, and R. O’Shaughnessy, Hierarchical black hole mergers in active galactic nuclei, Phys. Rev. Lett. 123, 181101 (2019).
  * Yang  _et al._ [2020] Y. Yang, I. Bartos, Z. Haiman, B. Kocsis, S. Márka, and H. Tagawa, Cosmic Evolution of Stellar-mass Black Hole Merger Rate in Active Galactic Nuclei, Astrophys. J. 896, 138 (2020), arXiv:2003.08564 [astro-ph.HE] .
  * Bartos  _et al._ [2017] I. Bartos, B. Kocsis, Z. Haiman, and S. Márka, Rapid and Bright Stellar-mass Binary Black Hole Mergers in Active Galactic Nuclei, Astrophys. J. 835, 165 (2017), arXiv:1602.03831 [astro-ph.HE] .
  * O’Leary  _et al._ [2009] R. M. O’Leary, B. Kocsis, and A. Loeb, Gravitational waves from scattering of stellar-mass black holes in galactic nuclei, Mon. Not. Roy. Astron. Soc. 395, 2127 (2009), arXiv:0807.2638 [astro-ph] .
  * Abbott  _et al._ [2023] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Phys. Rev. X 13, 011048 (2023), arXiv:2111.03634 [astro-ph.HE] .
  * Randall and Xianyu [2019] L. Randall and Z.-Z. Xianyu, A Direct Probe of Mass Density Near Inspiraling Binary Black Holes, Astrophys. J. 878, 75 (2019), arXiv:1805.05335 [gr-qc] .
  * Inayoshi  _et al._ [2017] K. Inayoshi, N. Tamanini, C. Caprini, and Z. Haiman, Probing stellar binary black hole formation in galactic nuclei via the imprint of their center of mass acceleration on their gravitational wave signal, Phys. Rev. D 96, 063014 (2017), arXiv:1702.06529 [astro-ph.HE] .
  * Cardoso  _et al._ [2021] V. Cardoso, F. Duque, and G. Khanna, Gravitational tuning forks and hierarchical triple systems, Phys. Rev. D 103, L081501 (2021), arXiv:2101.01186 [gr-qc] .
  * Lim and Rodriguez [2020] H. Lim and C. L. Rodriguez, Relativistic three-body effects in hierarchical triples, Phys. Rev. D 102, 064033 (2020), arXiv:2001.03654 [astro-ph.HE] .
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin induced orbital precession and its modulation of the gravitational wave forms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * Liu  _et al._ [2019] B. Liu, D. Lai, and Y.-H. Wang, Binary Mergers near a Supermassive Black Hole: Relativistic Effects in Triples, Astrophys. J. Lett. 883, L7 (2019), arXiv:1906.07726 [astro-ph.HE] .
  * Liu and Lai [2022] B. Liu and D. Lai, Probing the Spins of Supermassive Black Holes with Gravitational Waves from Surrounding Compact Binaries, Astrophys. J. 924, 127 (2022), arXiv:2105.02230 [astro-ph.HE] .
  * Yu and Chen [2021] H. Yu and Y. Chen, Direct determination of supermassive black hole properties with gravitational-wave radiation from surrounding stellar-mass black hole binaries, Phys. Rev. Lett. 126, 021101 (2021), arXiv:2009.02579 [gr-qc] .
  * Romero-Shaw  _et al._ [2023] I. M. Romero-Shaw, N. Loutrel, and M. Zevin, Inferring interference: Identifying a perturbing tertiary with eccentric gravitational wave burst timing, Phys. Rev. D 107, 122001 (2023), arXiv:2211.07278 [gr-qc] .
  * Meiron  _et al._ [2017] Y. Meiron, B. Kocsis, and A. Loeb, Detecting triple systems with gravitational wave observations, Astrophys. J. 834, 200 (2017), arXiv:1604.02148 [astro-ph.HE] .
  * Hall  _et al._ [2021] E. D. Hall _et al._ , Gravitational-wave physics with Cosmic Explorer: Limits to low-frequency sensitivity, Phys. Rev. D 103, 122004 (2021), arXiv:2012.03608 [gr-qc] .
  * Harms  _et al._ [2013] J. Harms, B. J. J. Slagmolen, R. X. Adhikari, M. C. Miller, M. Evans, Y. Chen, H. Müller, and M. Ando, Low-Frequency Terrestrial Gravitational-Wave Detectors, Phys. Rev. D 88, 122003 (2013), arXiv:1308.2074 [gr-qc] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Luo  _et al._ [2016] J. Luo _et al._ (TianQin), TianQin: a space-borne gravitational wave detector, Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Hu and Wu [2017] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity, Natl. Sci. Rev. 4, 685 (2017).
  * Ruan  _et al._ [2020] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Taiji program: Gravitational-wave sources, Int. J. Mod. Phys. A 35, 2050075 (2020), arXiv:1807.09495 [gr-qc] .
  * Kawamura  _et al._ [2021] S. Kawamura _et al._ , Current status of space gravitational wave antenna DECIGO and B-DECIGO, PTEP 2021, 05A105 (2021), arXiv:2006.13545 [gr-qc] .
  * Sato  _et al._ [2017] S. Sato _et al._ , The status of DECIGO, J. Phys. Conf. Ser. 840, 012010 (2017).
  * Kuns [2019] K. A. Kuns,  _Future Networks of Gravitational Wave Detectors: Quantum Noise and Space Detectors_ , Ph.D. thesis, University of California, Santa Barbara (2019).
  * Kuns  _et al._ [2020] K. A. Kuns, H. Yu, Y. Chen, and R. X. Adhikari, Astrophysics and cosmology with a decihertz gravitational-wave detector: TianGO, Phys. Rev. D 102, 043001 (2020), arXiv:1908.06004 [gr-qc] .
  * Maggiore [2007] M. Maggiore,  _Gravitational Waves. Vol. 1: Theory and Experiments_ , Oxford Master Series in Physics (Oxford University Press, 2007).
  * B.M. Peterson [2014] B.M. Peterson, Measuring the Masses of Supermassive Black Holes, Space Sci. Rev. 183 (2014).
  * Barth  _et al._ [2001] A. J. Barth, M. Sarzi, H.-W. Rix, L. C. Ho, A. V. Filippenko, and W. L. W. Sargent, Evidence for a Supermassive Black Hole in the S0 Galaxy NGC 3245, Astrophys. J. 555, 685 (2001), arXiv:astro-ph/0012213 [astro-ph] .
  * Onken  _et al._ [2014] C. A. Onken, M. Valluri, J. S. Brown, P. J. McGregor, B. M. Peterson, M. C. Bentz, L. Ferrarese, R. W. Pogge, M. Vestergaard, T. Storchi-Bergmann, and R. A. Riffel, The Black Hole Mass of NGC 4151. II. Stellar Dynamical Measurement from Near-infrared Integral Field Spectroscopy, Astrophys. J. 791, 37 (2014), arXiv:1406.6735 [astro-ph.GA] .
  * Bentz  _et al._ [2021] M. C. Bentz, P. R. Williams, R. Street, C. A. Onken, M. Valluri, and T. Treu, A Detailed View of the Broad-line Region in NGC 3783 from Velocity-resolved Reverberation Mapping, Astrophys. J. 920, 112 (2021), arXiv:2108.00482 [astro-ph.GA] .
  * Bentz  _et al._ [2022] M. C. Bentz, P. R. Williams, and T. Treu, The Broad Line Region and Black Hole Mass of NGC 4151, Astrophys. J. 934, 168 (2022), arXiv:2206.03513 [astro-ph.GA] .
  * Akiyama  _et al._ [2019] K. Akiyama _et al._ (Event Horizon Telescope), First M87 Event Horizon Telescope Results. VI. The Shadow and Mass of the Central Black Hole, Astrophys. J. Lett. 875, L6 (2019), arXiv:1906.11243 [astro-ph.GA] .
  * Yu  _et al._ [2021] H. Yu, Y. Wang, B. Seymour, and Y. Chen, Detecting gravitational lensing in hierarchical triples in galactic nuclei with space-borne gravitational-wave observatories, Phys. Rev. D 104, 103011 (2021), arXiv:2107.14318 [gr-qc] .
  * Sberna  _et al._ [2022] L. Sberna _et al._ , Observing GW190521-like binary black holes and their environment with LISA, Phys. Rev. D 106, 064056 (2022), arXiv:2205.08550 [gr-qc] .
  * Kuntz and Leyde [2023] A. Kuntz and K. Leyde, Transverse Doppler effect and parameter estimation of LISA three-body systems, Phys. Rev. D 108, 024002 (2023), arXiv:2212.09753 [gr-qc] .
  * Heggie [1975] D. C. Heggie, Binary evolution in stellar dynamics, Mon. Not. Roy. Astron. Soc. 173, 729 (1975).
  * 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 .
  * Martens and Joffre [2021] W. Martens and E. Joffre, Trajectory Design for the ESA LISA Mission 10.1007/s40295-021-00263-2 (2021), arXiv:2101.03040 [gr-qc] .
  * Cutler  _et al._ [1994] C. Cutler, D. Kennefick, and E. Poisson, Gravitational radiation reaction for bound motion around a schwarzschild black hole, Physical Review D 50, 3816–3835 (1994).
  * Poisson and Will [2014] E. Poisson and C. M. Will,  _Gravity: Newtonian, Post-Newtonian, Relativistic_ (Cambridge Univ. Press, 2014).
  * Hagihara [1930] Y. Hagihara, Theory of the Relativistic Trajeetories in a Gravitational Field of Schwarzschild, Japanese Journal of Astronomy and Geophysics 8, 67 (1930).
  * Darwin [1959] C. Darwin, The Gravity Field of a Particle, Proceedings of the Royal Society of London Series A 249, 180 (1959).
  * Darwin [1961] C. Darwin, The Gravity Field of a Particle. II, Proceedings of the Royal Society of London Series A 263, 39 (1961).
  * Chandrasekhar [1985] S. Chandrasekhar,  _The mathematical theory of black holes_ (1985).
  * Will [2018] C. M. Will, New General Relativistic Contribution to Mercury’s Perihelion Advance, Phys. Rev. Lett. 120, 191101 (2018), arXiv:1802.05304 [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 .
  * Cutler [1998] C. Cutler, Angular resolution of the LISA gravitational wave detector, Phys. Rev. D 57, 7089 (1998), arXiv:gr-qc/9703068 .
  * Droz  _et al._ [1999] S. Droz, D. J. Knapp, E. Poisson, and B. J. Owen, Gravitational waves from inspiraling compact binaries: Validity of the stationary phase approximation to the Fourier transform, Phys. Rev. D 59, 124016 (1999), arXiv:gr-qc/9901076 .
  * Christian  _et al._ [2018] P. Christian, S. Vitale, and A. Loeb, Detecting Stellar Lensing of Gravitational Waves with Ground-Based Observatories, Phys. Rev. D 98, 103022 (2018), arXiv:1802.02586 [astro-ph.HE] .
  * Tagawa  _et al._ [2020] H. Tagawa, Z. Haiman, I. Bartos, and B. Kocsis, Spin Evolution of Stellar-mass Black Hole Binaries in Active Galactic Nuclei, Astrophys. J. 899, 26 (2020), arXiv:2004.11914 [astro-ph.HE] .
  * Suzuki  _et al._ [2020] H. Suzuki, P. Gupta, H. Okawa, and K.-i. Maeda, Post-Newtonian Kozai–Lidov mechanism and its effect on cumulative shift of periastron time of binary pulsar, Mon. Not. Roy. Astron. Soc. 500, 1645 (2020), arXiv:2006.11545 [gr-qc] .
  * Knee  _et al._ [2023a] A. Knee, J. McIver, S. Naoz, B.-M. Hoang, and I. Romero-Shaw, Characterizing LISA-band gravitational-wave signals from eccentric compact binaries orbiting a supermassive black hole (2023a), Presented at APS April Meeting.
  * Knee  _et al._ [2023b] A. Knee _et al._ , in prep., (2023b).
  * Naoz [2016] S. Naoz, The eccentric kozai-lidov effect and its applications, Annual Review of Astronomy and Astrophysics 54, 441 (2016).
  * Torres-Orjuela  _et al._ [2020] A. Torres-Orjuela, X. Chen, and P. Amaro-Seoane, Phase shift of gravitational waves induced by aberration, Phys. Rev. D 101, 083028 (2020), arXiv:2001.00721 [astro-ph.HE] .
  * Torres-Orjuela and Chen [2023] A. Torres-Orjuela and X. Chen, Moving gravitational wave sources at cosmological distances: Impact on the measurement of the Hubble constant, Phys. Rev. D 107, 043027 (2023), arXiv:2210.09737 [astro-ph.CO] .
  * Bender and Orszag [1999] C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers I_ (Springer New York, 1999).
  * Chamberlain  _et al._ [2019] K. Chamberlain, C. J. Moore, D. Gerosa, and N. Yunes, Frequency-domain waveform approximants capturing Doppler shifts, Phys. Rev. D 99, 024025 (2019), arXiv:1809.04799 [gr-qc] .

##  Appendix A Validity of Waveform Approximations

We consider some approximations that are made in the formulation of the
frequency domain waveform. For the source frame waveform
$h_{s}(t_{s})=A_{s}(t_{s})e^{-i\Phi_{s}(t_{s})}$, we define the connection
between time and gravitational wave frequency in the source frame by

$\dot{\Phi}_{s}(t_{s})=2\pi f_{s}.$ (43)

The outer orbit of the BBH induces a change in the time at which a GW of a
particular frequency reaches a fixed observer, which we denote $t_{o}$.
Clearly marking the dependencies of various times on one another,

$t_{o}(t_{s})=t_{s}+t_{\parallel}(t_{s})+D_{L}/c,$ (44)

where $t_{\parallel}=a_{o}\sin\iota_{J}\sin\varphi$ is the time is takes for a
radiated GW to propagate from the BBH to the SMBH along the direction of
$\hat{N}$, assuming a circular outer orbit.

With $h_{s}$ and $h_{o}$ being the strain in the source and observer frames,
we then have

$h_{o}(t_{o})=h_{s}(t_{s})=h_{s}(t_{o}-t_{\parallel}-D_{L}/c),$ (45)

which in the Fourier transform (as a function of the observed frequency
$f_{o}$) becomes

$\tilde{h}_{o}(f_{o})=\int h_{s}(t_{o}-t_{\parallel}-D_{L}/c)e^{2\pi
if_{o}t_{o}}dt_{o}.$ (46)

Inputting the form of $h_{s}(t_{s})$ gives

$\tilde{h}_{o}(f_{o})=\int
A_{s}(t_{o}-t_{\parallel}-D_{L}/c)e^{-i\Phi_{s}(t_{o}-t_{\parallel}-D_{L}/c)}e^{2\pi
if_{o}t_{o}}dt_{o},$ (47)

and assigning $t_{o}-D_{L}/c=t$ so that $t_{s}=t-t_{\parallel}$ produces

$\tilde{h}_{o}(f_{o})=e^{2\pi if_{o}D_{L}/c}\int
A_{s}(t-t_{\parallel})e^{-i\Phi_{s}(t-t_{\parallel})}e^{2\pi if_{o}t}dt.$ (48)

For a typical system we study (e.g., $M_{3}=10^{8}M_{\odot}$,
$a_{o}=100M_{3}$), $t_{\parallel}\sim 10^{4}-10^{5}$ seconds (depending on the
orbital angle) and $\dot{t}_{\parallel}\sim t_{\parallel}\Omega_{o}\sim
0.01-0.1$.

We can make a number of simplifications to this expression. First,
$\Phi_{s}(t-t_{\parallel})\approx\Phi_{s}(t)-2\pi f_{s}t_{\parallel}$ as long
as $\dot{f}_{s}t_{\parallel}\ll f_{s}$. For an inspiral regime BBH with two
50$M_{\odot}$ BH, $f_{s}/\dot{f}_{s}\sim 3\times 10^{3}f_{s}^{-8/3}$s [37],
and with the majority of the time-integration taking place with $f_{s}\lesssim
0.1$Hz, the approximation using $\dot{f}_{s}t_{\parallel}\ll f_{s}$ holds
well.

We apply a similar approximation to $A_{s}(t-t_{\parallel})$. The time scale
for the evolution of the GW amplitude is roughly [57, 37]
$A_{s}/\dot{A}_{s}\sim\frac{3}{2}\frac{f_{s}}{\dot{f}_{s}}\gg t_{\parallel}$,
so we can reasonably approximate $A_{s}(t-t_{\parallel})\approx A_{s}(t)$.
This simplifies Eq. (48) to

$\tilde{h}_{o}(f_{o})\approx e^{2\pi if_{o}D_{L}/c}\int
A_{s}(t)e^{-i[\Phi_{s}(t)-2\pi f_{s}t_{\parallel}(t)]}e^{2\pi if_{o}t}dt.$
(49)

Consider an expansion of $F(t)=\Phi_{s}(t)-2\pi f_{s}t_{\parallel}(t)$ around
some time $t^{\prime}$. Noting that $t_{\parallel}$ is a function of $t$, we
find

$\Phi_{s}(t)-2\pi f_{s}t_{\parallel}(t)\approx\Phi_{s}(t^{\prime})-2\pi
f_{s}t_{\parallel}(t^{\prime})\\\
+2\pi(f_{s}-f_{s}\dot{t}_{\parallel}-\dot{f}_{s}t_{\parallel})(t-t^{\prime})+\frac{1}{2}\ddot{F}(t^{\prime})(t-t^{\prime})^{2}+...$
(50)

Since $\dot{f}_{s}/f_{s}\ll\dot{t}_{\parallel}/t_{\parallel}$, the linear term
approximates to $2\pi f_{s}(1-\dot{t}_{\parallel})$. Now, to complete the
Fourier transform, we turn to the stationary phase approximation. Namely, the
majority of the integral comes from the region where the argument of the
oscillating term is stationary, which occurs at a time $\tau$ when
$f_{o}-f_{s}(\tau)(1-\dot{t}_{\parallel}(\tau))=0$.

Inserting the second order expansion found above into the Fourier transform
gives

$\tilde{h}_{o}(f_{o})\approx e^{2\pi if_{o}D_{L}/c}\int A_{s}(\tau)\\\
\exp\Bigr{\\{}-i[-2\pi f_{o}t+\Phi_{s}(\tau)-2\pi
f_{s}(\tau)t_{\parallel}(\tau)\\\ +2\pi
f_{s}(\tau)(1-\dot{t}_{\parallel}(\tau))(t-\tau)\\\
+\frac{1}{2}\ddot{F}(\tau)(t-\tau)^{2}+...]\Bigr{\\}}dt,$ (51)

which simplifies to

$\tilde{h}_{o}(f_{o})\approx e^{2\pi
i[f_{o}D_{L}/c+f_{s}(\tau)t_{\parallel}(\tau)]}e^{2\pi
if_{s}(\tau)(1-\dot{t}_{\parallel}(\tau))\tau}\\\ \int
A_{s}(\tau)e^{-i[\Phi_{s}(\tau)+\frac{1}{2}\ddot{F}(\tau)(t-\tau)^{2}+...]}dt,$
(52)

recalling that $f_{s}(\tau)(1-\dot{t}_{\parallel}(\tau))=f_{o}$.

The expression $\ddot{F}$ has terms proportional to $\dot{f}_{s}$,
$\ddot{f}_{s}t_{\parallel}$, $\dot{f}_{s}\dot{t}_{\parallel}$, and
$f_{s}\ddot{t}_{\parallel}$. For our typical system, the three latter terms
are generally much smaller than the first, allowing a reasonable approximation
of $\ddot{F}\sim\ddot{\Phi}_{s}$.

Carrying out this integral is a standard exercise as in [68], and we see that
Eq. (52) evaluates to

$\tilde{h}_{o}(f_{o})\approx\frac{\sqrt{\pi}}{2}e^{2\pi
if_{o}D_{L}/c+i\dot{F}(\tau(f))\tau-iF(\tau(f))+i\pi/4}\\\ \times
A_{s}(\tau(f))\Bigr{(}\frac{2}{\ddot{F}_{s}(\tau(f))}\Bigr{)}^{1/2},$ (53)

recalling $F=\Phi_{s}-2\pi f_{s}t_{\parallel}$ and noting that Eq. (18) can be
used to convert time-dependent quantities into frequency dependent ones. It is
clear that the Doppler shift is built directly into this definition – the
components of the source radiation (e.g. $A_{s}$,
$\ddot{F}\approx\ddot{\Phi}_{s}$) which appear in the observed strain are
evaluated at a time $\tau$ when the source radiation is emitted at frequency
$f_{s}=f/(1-\dot{t}_{\parallel})$. The ratio of the observed frequency to
source frequency matches the expansion of the exact form of the Doppler shift
given in Eq. (55) below (in the low $\dot{t}_{\parallel}$ limit).

This result poses a problem for carrying out Fisher matrix calculations. The
Fisher matrix formalism relies on integration over the observed frequencies;
while there exists a monotonic relation between the source GW frequency and
time, the inclusion of the Doppler frequency shift results in the same
observed frequency originating from multiple distinct source frequencies.
Furthermore, the PSD term must be evaluated at the observed frequency, while
standard results for the GW signal in the frequency domain are parameterized
by the source frequency. Without a one-to-one relationship between the
observed frequency and source frequency, carrying out the Fisher matrix
calculations requires careful attention to these subtleties when evaluating
the frequency domain integrand. Future iterations of this analysis will more
carefully implement the Doppler shift in comparing the GW signal in the source
and observer frames.

The best path forward we can take to produce a reasonable approach for the
Fisher matrix computations is to make the approximation that
$\dot{t}_{\parallel}\ll 1$. With this simplification (which holds fairly well
across the majority of our parameter space, as we describe below),
$f_{obs}\approx f_{src}$ and thus the following simple relation emerges:

$\tilde{h}_{o}(f)\approx e^{2\pi if(D_{L}/c+t_{\parallel})}\tilde{h}_{s}(f).$
(54)

Let us more completely explain why we can reasonably neglect the change in GW
frequency due to the Doppler shift induced by the BBH orbital velocity around
the SMBH. The (exact) longitudinal Doppler shifted frequency is given by

$f_{obs}=f_{src}\frac{\sqrt{1+\beta_{\parallel}}}{\sqrt{1-\beta_{\parallel}}},$
(55)

where $f_{obs}$ and $f_{src}$ are the observed and source frequencies, and
$\beta_{\parallel}=v_{\parallel}/c=\dot{t}_{\parallel}$ is the source velocity
along the line of sight. Using the methods of Sec. II.2, the maximum orbital
velocity occurs at pericenter, with

$\beta=\frac{p-2-2e_{o}}{p\sqrt{(p-2)^{2}-4e_{o}^{2}}}(1+e_{o})\sqrt{p-6+2e_{o}}.$
(56)

Depending on the argument of pericenter, the magnitude of the source’s line of
sight velocity can reach up to this value. In the parameter space we study,
this velocity is maximized over eccentricity when $e_{o}=0.9$. At this
eccentricity, the maximum velocity over semimajor axes in our parameter space
occurs when $p\approx 12$ with $\beta_{max}\approx 0.37$. Then, the largest
increase in GW frequency due to the Doppler shift is roughly 50%.

Making the approximation that $f_{obs}\approx f_{src}$ can produce some
inaccuracies in the distance-accumulated phase terms in Eq. (53); however,
because the Fisher matrix formalism includes integration of the product of the
waveform and its complex conjugate, the value of the accumulated phase does
not have any effect on the results. Furthermore, the sensitivity curves in
Fig. 2 vary slowly in frequency, so therefore, our approximation that the
source frame and observer frame GW frequencies are roughly equal does not
significantly affect our Fisher matrix calculations and resulting parameter
uncertainties.

Future work may implement more rigorous treatment of the waveform in the
frequency domain, including corrections suggested in [69], for example.
However, we expect that the information provided by examining the shifts in
observed GW frequencies, which track the BBH orbital velocity, is essentially
already provided by the Doppler phase, which tracks the BBH orbital position.
As such, we anticipate little improvement in parameter measurability by
including this additional effect.

##  Appendix B Description of Measurement Accuracy

In Fig. 7, we note two power laws for the measurement accuracy of $M_{3}$. We
say that the pericenter precession gives $(a_{o}/M_{3})^{3/2}$, while de
Sitter precession scales like $a_{o}^{5/2}/M_{3}^{3/2}$. Below, we give
scaling arguments to understand this plot. As described in Sec. IV, we measure
$a_{o}\times\sin\iota_{J}$ and $\Omega_{o}$ very well via the Doppler shift

$\Phi_{D}=2\pi f\hat{N}\cdot\vec{r}=2\pi fr\sin\iota_{J}\sin\phi\,,$ (57)

however we need an additional effect to break the degeneracy between $a_{o}$
and $\sin\iota_{J}$. In this appendix, we will show how the effect of
precession can be understood as separating the radial and azimuthal periods
$P_{r}\approx P_{\phi}-\\#\frac{M}{a}$, and breaks the degeneracy. We will
also explain how the Thomas phase and polarization phase terms allow us to
also break the degeneracy.

Let us study a simple analytic toy model for precession where the radial
period is shortened by a 1PN term. We will consider only a waveform with
Doppler phase here to extract the physical reason that pericenter precession
helps us measure the SMBH’s orbit. We set

$\tilde{h}(f)=\exp\Bigr{[}i(\Phi_{D}+2\pi ft_{c})\Bigr{]},$ (58)

and provide the following simple dynamics to the outer orbit:

$\Phi_{D}(t)=2\pi
f\sin\iota_{J}\frac{a_{o}(1-e_{o})^{2}}{1+e_{o}\cos\xi_{r}}\sin(\xi_{\varphi}),$
(59)
$\xi_{r}=\xi_{r,0}+\Omega_{o}\left(1-\epsilon\frac{3M_{3}}{a_{o}(1-e_{o}^{2})}\right)t,$
(60) $\xi_{\varphi}=\xi_{r,0}+\xi_{\varphi,0}+\Omega_{o}t\,,$ (61)

where the variable $\epsilon$ is a counting parameter for precession, and set
to 0 or 1 at the end of the calculation. This system gives a rough
approximation for a low eccentricity orbit that includes relativistic
pericenter precession, assuming $1/\Omega_{o}\ll t_{obs}$ so that many orbits
are completed during $t_{obs}$ and thus the angular velocities for the orbital
and precessional motion average out to their secular values. The quantities
$\xi_{r,0}$ and $\xi_{\varphi,0}$ are analogous to $\phi_{0}$ and
$\gamma_{o}$, respectively.

Since this is such a simple model, we can compute the Fisher matrix
analytically under certain assumptions. We first compute the Fisher matrix
with elements $\\{M_{3},\Omega_{o},\sin\iota_{J}\\}$. After computing the
derivatives of $\tilde{h}$, we expand in $e_{o}\ll 1$ and consider only the
secular effect of the trigonometric functions in $\xi_{r}$ and
$\xi_{\varphi}$, as the contributions from oscillating terms will be minimal
after integration over $f$. We substitute $t\propto\mathcal{M}^{-5/3}f^{-8/3}$
(cf. Eq. (18)) and after integrating over $f$ with a flat PSD, we find that
the resulting matrix is invertible only if we have some pericenter precession
$\epsilon\neq 0$. The fractional error in mass then scales like

$\frac{\Delta
M_{3}}{M_{3}}\propto\frac{1}{\epsilon}\Bigr{(}\frac{a_{o}}{M_{3}}\Bigr{)}^{3/2}.$
(62)

This result shows that pericenter precession will produce contours of constant
$\Delta M_{3}/M_{3}$ which scale with $(a_{o}/M_{3})^{3/2}$, which is in good
agreement with the lines in Fig. 3, for example. These lines do not appear
when $e_{o}=0$, however, because the difference in the radial and azimuthal
frequencies does not produce any change in the BBH’s actual orbit when that
orbit is circular.

We also expanded the Fisher matrix dimensions to include
$\\{M_{3},\Omega_{o},\sin\iota_{J},\xi_{r,0},\xi_{\varphi,0},t_{c}\\}$ and
found identical scaling in $\Delta M_{3}/M_{3}$.

In contrast, let us examine the behavior of dS precession. We will now show
how mass and orbital frequency can be independently measured by the the Thomas
phase $\Phi_{T}$ and the polarization phase $\Phi_{P}$. As discussed in [21],
the Thomas phase for a circular orbit is approximately 444This comes from Eq.
(65) of [21]. There are multiple expressions for the average rate depending on
the sign and dot products of $(\hat{L},\hat{N})$, which only differ by factors
of of order unity.

$\Phi_{T}\sim 2\pi(1+\cos\lambda_{L})\Omega_{\text{dS}}t\,,$ (63)

where
$\Omega_{\text{dS}}=\frac{3}{2}\frac{M_{3}\Omega_{o}}{a_{o}}\approx\frac{3}{2}M_{3}^{2/3}\Omega_{o}^{5/3}$.

The Thomas phase breaks the degeneracy between $M_{3}$, $\Omega_{o}$, and
$\sin\iota_{J}$. If we consider a waveform with perfectly measured
$\Omega_{o}$, with a simple waveform $h\propto e^{i\Phi_{T}}$ one can show the
measurement scales as

$\frac{\Delta M_{3}}{M_{3}}\approx\sqrt{\frac{1}{\Gamma_{\log M_{3}\log
M_{3}}}}\propto\frac{1}{\partial_{\log M_{3}}\Omega_{\text{dS}}}\\\
=\frac{1}{M_{3}^{2/3}\Omega_{o}^{2}}\propto\frac{a_{o}^{5/2}}{M_{3}^{3/2}}$
(64)

For contour plots such as Fig. 3, this result matches the power law for lines
of constant $\Delta M_{3}/M_{3}$ found in the region where de Sitter
precession is most rapid.

The polarization phase breaks the degeneracy between $M_{3}$ and $\Omega_{o}$
in the same manner as the Thomas phase. The power law can be seen by
considering a waveform $h\propto e^{i\Phi_{P}}$. The polarization phase is
defined as

$\Phi_{P}(t)=-\arctan\left[\frac{A_{\times}(t)F_{\times}(t)}{A_{+}(t)F_{+}(t)}\right].$
(65)

In the context of the Fisher matrix, it is useful to compute the derivative of
this phase,

$\partial_{\theta}\Phi_{P}=-\frac{1}{(1+\tan^{2}\Phi_{P})}\partial_{\theta}\Bigr{(}\frac{A_{\times}(t)F_{\times}(t)}{A_{+}(t)F_{+}(t)}\Bigr{)}.$
(66)

The fraction of amplitude factors and antenna patterns depends on the angles
$\theta_{S},\phi_{S},\theta_{J},\phi_{J},\lambda_{L},\text{ and }\alpha_{0}$,
as well as the integrated de Sitter precession rate, all of which appear
exclusively in trigonometric functions (cf. Eqs. (8), (19), (20), (21), (22),
and (23)). We assume the instantaneous de Sitter precession rate does not vary
significantly, assigning the rate to its secular value
$\Omega_{dS}\approx\frac{3}{2}\frac{M_{3}\Omega_{o}}{a_{o}}$.

Taking $\theta=\log M_{3}$, then, the only dependence on $M_{3}$ in this
fraction is through $\Omega_{dS}$. Since $\Omega_{dS}$ appears only in the
argument of sines and cosines, we expect that after applying the chain rule
$\partial_{\log M_{3}}\rightarrow(\partial\Omega_{dS}/\partial\log
M_{3})\partial_{\Omega_{dS}}$, the magnitude of $\partial_{\log
M_{3}}\Phi_{P}$ is primarily influenced by $\partial\Omega_{dS}/\partial\log
M_{3}$. This expression evaluates to (cf. Eq. (64))

$\frac{\partial\Omega_{dS}}{\partial\log
M_{3}}=\frac{M_{3}^{3/2}}{a_{o}^{5/2}},$ (67)

Thus, assuming $\Omega_{o}$ is known, the Fisher matrix will scale like

$\Gamma_{\log M_{3}\log M_{3}}=\int df\left(\partial_{\log
M_{3}}\Phi_{P}(f)\right)^{2}\propto(\partial_{\log
M_{3}}\Omega_{\text{dS}})^{2}.$ (68)

Therefore, considering both $\Phi_{T}$ and $\Phi_{P}$, de Sitter precession
produces an uncertainty in $\log M_{3}$ that scales roughly as
$a_{o}^{5/2}/M_{3}^{3/2}$, following lines of constant $\Omega_{dS}$.

It may seem counter-intuitive that two processes with identically-scaling
rates (i.e., pericenter precession and de Sitter precession) produce different
power laws in the shape of the contour lines. However, even though the rate of
pericenter precession $\Omega_{\text{peri prec}}\propto
M_{3}^{2/3}\Omega_{o}^{5/3}$ occurs at the same PN order as de Sitter
precession $\Omega_{\text{dS}}\sim M_{3}^{2/3}\Omega_{o}^{5/3}$, there is an
additional factor of $a_{o}$ in the pericenter precession contribution due to
the Doppler shift term being proportional to $r$. This causes the power law
dependence for the measurement accuracy of $M_{3}$ to scale differently by a
factor of $a_{o}$.