# Detecting environmental effects in gravitational waves from binaries
perturbed by periodic forces

Lorenz Zwick1<EMAIL_ADDRESS>Christopher Tiede1 Alessandro A. Trani1
Andrea Derdzinski2,3 Zoltan Haiman4,5 Daniel J. D’Orazio1 Johan Samsing1
1Niels Bohr International Academy, The Niels Bohr Institute, Blegdamsvej 17,
DK-2100, Copenhagen, Denmark
2Department of Life and Physical Sciences, Fisk University, 1000 17th Avenue
N., Nashville, TN 37208, USA
3Department of Physics & Astronomy, Vanderbilt University, 2301 Vanderbilt
Place, Nashville, TN 37235, USA
4Department of Astronomy, Columbia University, 550 West 120th Street, New
York, NY 10027, USA
5Department of Physics, Columbia University, 550 West 120th Street, New York,
NY 10027, USA

###### Abstract

We study the gravitational wave (GW) emission of sources perturbed by periodic
dynamical forces which do not cause secular evolution in the orbital elements.
We construct a corresponding post-Newtonian waveform model and provide
estimates for the detectability of the resulting GW phase perturbations, for
both space-based and future ground-based detectors. We validate our results by
performing a set of Bayesian parameter recovery experiments with post-
Newtonian waveforms. We find that, in stark contrast to the more commonly
studied secular dephasing, periodic phase perturbations do not suffer from
degeneracies with any of the tested vacuum binary parameters. We discuss the
applications of our findings to a range of possible astrophysical scenarios,
finding that such periodic perturbations may be detectable for massive black
hole binaries embedded in circum-binary discs, extreme mass-ratio inspirals in
accretion discs, as well as stellar-mass compact objects perturbed by tidal
fields. We argue that modelling conservative sub-orbital dynamics opens up a
promising new avenue to detect environmental effects in binary sources of GWs
that should be included in state-of-the-art waveform templates.

##  I Introduction

Astrophysical sources of gravitational waves (GW) are likely to form in high
density environments, such as star clusters, nuclear clusters, and accretion
or circumbinary discs. Understanding and modelling the physical mechanisms at
play in these environments is crucial, as they are responsible for determining
the efficiency of the various formation channels for black hole binaries
(BHB), and ultimately result in the expected event rates of detectable signals
for current and future GW detectors [see the following reviews and white
papers, among the innumerable relevant works, 1, 46, 85, 4, 47, 33]. In
addition to determining the rates, different astrophysical environments set
characteristic ranges and correlations in the expected parameters of observed
BHBs, as exemplified by the numerous efforts to link LIGO/Virgo/Kagra
observations to specific formation channels [16, 91, 110, 6, 142, 71, 154, 73,
117].

In almost all astrophysically motivated scenarios, the direct influence of the
source’s environment on the emitted GWs is expected to be extremely small once
BHBs reach the separations required to produce detectable GWs [with the
possible exception of Pulsar Timing Array sources, see, e.g., 65, 119].
However, extremely small need not necessarily mean negligible. These
influences are referred to with the umbrella term of ”environmental effects”
(EE), which are defined as small perturbations with respect to the expected
vacuum GWs that arise as a consequence of the source’s astrophysical
environment. Common examples of EEs include dynamical perturbations caused by
third bodies, effects of peculiar velocity and other accelerations, as well as
the study of frictional forces motivated by both gaseous and dark matter
backgrounds [13, 68, 88, 20, 135, 99, 49, 80, 149, 62, 32, 28, 155, 131].

Since the early mentions in the nineties [25, 108] to the seminal work of the
mid two-thousands [78, 74, 12], the study of EE, i.e. astrophysically
motivated perturbations to GW emission, has been gaining traction due to the
ever closer dawn of space-based GW detectors sensitive to the milli-Hz band
[in particular the recently adopted LISA mission, 33]. The scientific
motivations are threefold: Firstly, the detection of even a single source that
is clearly perturbed by EEs would allow to set unprecedented constraints on
the surroundings of BHBs and the physical processes that facilitate their
formation and shape their evolution [41, 38, 156, 124, 32, 42, 62, 114].
Secondly, the presence of unmodelled EEs may produce biases in the recovered
binary parameters (see references above). Thirdly, deviations from the
expected vacuum evolution of sources may be mislabelled as resulting from non-
standard physics, in the case that astrophysical influences are not modelled
correctly [see, e.g., 90, 151, 23, 22, 32, 155]. The latter point is
especially relevant when interpreting signals from extreme mass-ratio
inspirals, which are supposed to provide extremely precise measurement of the
metric in the close vicinity of massive BHs [11, 60, 3, 97]. In short, the
careful study of EE is is required to both preserve and broaden the scientific
goals of future GW detectors in the joint fields of fundamental physics and
astrophysics.

Focusing now on binary sources of GWs, the state of the art in the modelling
of EE currently consists of the implementation of simple secular de-phasing
prescriptions in pre-existing vacuum waveform templates, as well as their
study through Bayesian parameter estimation techniques [136, 32, 124, 61].
Indeed, the observable that is most often associated to EE is that of a phase
difference between the perturbed waveform and a reference vacuum waveform,
which slowly accumulates throughout the binary source’s inspiral phase. The
main motivations for focusing on de-phasing are twofold: Firstly, the process
of matched filtering is known to be extremely sensitive to cumulative
perturbations to the GW phase. Phase shifts $\delta\phi$ as small as:

$\displaystyle\delta\phi\sim\text{few}/\text{SNR},$ (1)

where SNR is the signal to noise ratio of the given GW, can significantly
reduce the match between a waveform template and the data and reduce the
effectiveness of parameter estimation [see, e.g., 34, 35, 74, 70, 92].
Secondly, most environmental effects may indeed be modelled with simple
secular de-phasing prescriptions, at least as a first approximation [see 155,
for an impression on the broad relevance of dephasing in the context of EE].
As an example, any coupling between the environment’s and the binary’s energy
will cause a de-phasing by the following relation:

$\displaystyle\delta\phi\sim\int\int\left(\Delta\dot{f}\,{\rm{d}}t\right)\,\,{\rm{d}}t,$
(2)

where $\Delta\dot{f}$ is a change in the binary’s chirp due to the coupling.
Additionally, any kind of EE that involves a Doppler shift [see, e.g., 26,
where modifications to the GW phase and GW amplitude are discussed in detail],
such as peculiar motion or gravitational red-shift, is also well approximated
via a dephasing [88, 68, 135, 29, 152]:

$\displaystyle\delta\phi\sim\int\Delta f\,{\rm{d}}t,$ (3)

where $\Delta f$ is the Doppler shift in frequency.

Secular de-phasing is often considered the most promising avenue to detect EE
in GWs, and several works have used the reference value of Eq. 1 to claim
detectability111Including, admittedly, the corresponding author’s own work..
While the latter criterion can provide a good estimate in certain limits,
several careful parameter inference studies have highlighted how de-phasing
suffers from many-fold intrinsic limitations which may greatly hinder the
prospective of correctly identifying it in realistic binary sources [24, 72,
92]. Firstly, secular dephasing is preferentially accumulated in the very
early inspiral stages of a binary’s observable evolution, where the binary
only chirps slowly. This corresponds to a rather low available signal-to-noise
(SNR) in which the perturbation may effectively be constrained. Secondly,
tracking an accumulated phase shift over a very long inspiral may bee
exceedingly hard due to the intrinsic challenges of matched-filtering [53].
Thirdly, secular de-phasing is similar to many other post-Newtonian (PN)
corrections to the binary’s vacuum evolution, which also add slowly
accumulating contributions to the total GW phase [see e.g. 18]. Finally, and
perhaps most crucially, it is the case that a hypothetical constant dephasing
term, i.e. a phase offset, has exactly the same form as the source’s initial
phase (as can be seen at a glance in, e.g., Eq. 28), leading to potential
degeneracies. The consequence of this is that only the rate of variation in
the de-phasing is truly informative, rather than the magnitude of the de-
phasing itself.

These intrinsic challenges raise the following questions: Are there any other
features in GWs from perturbed binaries that may carry detailed information
regarding the source’s interaction with its environment? Are any of these
features less prone to the inherent limitations of secular dephasing when it
comes to parameter estimation? And is there any expectation that GW detectors
may be able to detect these in the not too distant future? Motivated by these
questions, we take a step away from modelling forces or perturbations that
would produce any long term evolution on the binary’s orbital elements, since
such forces ultimately appear in the GW as various flavours of secular
dephasing. Instead, we focus entirely on forces that may influence the
binary’s motion on orbital or even sub-orbital timescales without causing any
long term evolution. We will find that such forces leave characteristic and
potentially extremely informative imprints in the GW emission of binaries.

This work is structured as follows: In section II, we solve the equations of
motion for a circular binary that is influenced by an arbitrary, periodic
perturbing force. We then use the latter to construct a ready-to-use post
Newtonian waveform model. In sections III and IV, we analyse the resulting
periodic perturbations in the GW phase and provide some detectability criteria
using both analytical and numerical tools. In section V we discuss a range of
plausible astrophysical scenarios in which such perturbations may be
detectable. Finally, we lay out concluding remarks in section VI.

Figure 1: A simple representation of the type of perturbations analysed in
this work. The sketch showcases a binary (represented by a total mass $M$ and
reduced mass $\mu$) being affected by time-varying forces (decomposed in their
tangential, radial and normal components). Here we only consider forces that
do not cause a secular evolution in the orbital elements. The latter influence
the binary’s instantaneous velocity and lever arm, which in turn modifies the
system’s GW emission. The GW strain will therefore present features associated
to the forces acting on the binary. The magnitude of the effects are greatly
exaggerated for clarity.

##  II Gravitational waves from binaries subject to time-varying forces

###  II.1 Motion of a circular binary under time varying forces

We consider the case of a binary source of GWs with arbitrary total mass $M$
and reduced mass $\mu$. The binary’s orbit is described by the solutions to
the six osculating elements $\left[p,e,\nu,\omega,\Omega,I\right]$222Semi-
latus rectum, eccentricity, true anomaly, argument of peri-apsis, longitude of
ascending node and inclination., which fix the coordinate position of the
reduced mass in space [following here the formalism detailed in 96]:

$\displaystyle\begin{pmatrix}x\\\ y\\\
z\end{pmatrix}=r\begin{pmatrix}\cos\Omega\cos\left(\omega+\nu\right)-\sin\Omega\sin\left(\omega+\nu\right)\cos
I\\\
\sin\Omega\cos\left(\omega+\nu\right)+\cos\Omega\sin\left(\omega+\nu\right)\cos
I\\\ \sin\left(\omega+\nu\right)\sin I\end{pmatrix},$ (4)

where the instantaneous radius is:

$\displaystyle r=\frac{p}{1+e\cos\nu},$ (5)

In the GW driven regime, most sources will be circularised by the effect of
radiation reaction. Therefore, we will operate in the limit of $e\to 0$ as a
first approximation, and comment on the possible effects of eccentricity
below.

We study the effect of a perturbing force ${\rm{d}}\mathbf{F}$ which is acting
on the binary, which represents the influence of the environment with respect
to the vacuum gravitational force. Here ${\rm{d}}\mathbf{F}$ has units of
acceleration, and may be decomposed into its radial, tangential and normal
components:

$\displaystyle{\rm{d}}\mathbf{F}=R\mathbf{\hat{r}}+T\mathbf{\hat{t}}+N\mathbf{\hat{z}},$
(6)

where $\mathbf{\hat{r}}$, $\mathbf{\hat{t}}$ and $\mathbf{\hat{z}}$ denote
unit vectors in the aforementioned directions (see Fig. 1).

In this work, we restrict our analysis to in-plane perturbative forces, as out
of plane displacements will only ever affect the orbital inclination angles.
Note that, if sufficiently strong, variations in the inclination angles may
well produce detectable modulations in the relative amplitude of the GW’s plus
and cross polarisations. We leave this effect for future investigation. We
treat the binary’s orbit in the Newtonian limit. While this assumption is
known to break down as the orbital velocity approaches relativistic values, it
is justifiable by two arguments: First, most EEs will affect the binary more
strongly at large separations, where orbital speeds are low and GW emission is
weaker. Second, the assumption is justified in a strictly perturbative sense,
as additional cross terms333Cross terms between environmental perturbations
and post Newtonian dynamics have been extensively studied in [145, 146], and
have interesting consequences on the secular evolution of triple systems and
orbits in accretion discs [157]. between the environmental perturbations and
relativistic dynamics will be suppressed with respect to the leading order
terms in $1/c^{2}$. Note that under these assumptions, the binary’s motion
will be confined to a unique plane, which we identify as the $(x,y)$ plane.

In the Newtonian limit, the evolution equations for the osculating orbital
elements are known from the seminal Lagrange planetary equations:

$\displaystyle\frac{{\rm d}p}{{\rm
d}t}=2\sqrt{\frac{p^{3}}{GM}}\frac{T}{\left(1+e\cos\nu\right)}$ (7)
$\displaystyle\frac{{\rm d}\nu}{{\rm d}t}$
$\displaystyle=\sqrt{\frac{GM}{p^{3}}}(1+e\cos\nu)^{2}$
$\displaystyle+\frac{1}{e}\sqrt{\frac{p}{GM}}\left[\cos\nu
R-\frac{2+e\cos\nu}{1+e\cos f}\sin\nu T\right]$ (8) $\displaystyle\frac{{\rm
d}e}{{\rm
d}t}=\sqrt{\frac{p}{GM}}\left[R\sin\nu+\frac{2\cos\nu+e(1+\cos^{2}\nu)}{1+e\cos\nu}T\right]$
(9)

$\displaystyle\frac{{\rm d}\omega}{{\rm d}t}$
$\displaystyle=\frac{1}{e}\sqrt{\frac{p}{GM}}\Big{[}-\cos\nu
R+\frac{2+e\cos\nu}{1+e\cos\nu}\sin\nu\,T\Big{]}$ (10)

which we can solve in the perturbative limit via standard techniques. Our goal
is to capture the effect of forces which do not cause a secular evolution in
the orbital elements of the binary. Therefore, we the follow the groundwork
laid out in [156] and describe a generic force as a sum of Fourier components
with a fundamental frequency equal to the binary’s orbital frequency $f_{\rm
K}$:

$\displaystyle{\rm{d}}\mathbf{F}^{\rm p}$
$\displaystyle=\text{Re}\Bigg{(}\sum_{n}\Big{[}B^{\rm{T}}_{n}(r)\exp\left(i2\pi
nf_{\rm K}t\right)\mathbf{\hat{t}}$
$\displaystyle+B^{\rm{R}}_{n}(r)\exp\left(i2\pi nf_{\rm
K}t\right)\mathbf{\hat{r}}\Big{]}\Bigg{)},$ (11)

where the complex coefficients $B^{\rm T}_{n}$ and $B^{\rm R}_{n}$ are only
functions of the binary’s separation and their magnitude and dependence on $r$
encodes the forces that are applied on the binary by the given EE. The form of
Eqs. II.1 is to be understood as a Fourier series in the limit in which the
binary residency time at the separation $r$ is much larger than one orbital
period, i.e. the adiabatic limit. In that case, we can also allow (with a
slight abuse of notation) the parameter $n$ to take non-integer values in
order to describe periodic perturbations that do not exactly match a multiple
of the system’s dynamical time. We can insert the perturbing force into the
Lagrange planetary equations to find analytical solutions to the orbital
elements. As it will be directly related to GW observables, we report the
explicit first order perturbative solution444Here we take the chance to note
an imprecision present in [156], in which the computations did not properly
propagate the perturbation in the Kepler equation and resulted in a mistake in
Eq. (33). However, the latter does not change the qualitative results in
[156], as the physical scaling remains unaffected. of the true longitude
$\theta=\nu+\omega$:

$\displaystyle\theta$ $\displaystyle=2\pi f_{\rm
K}t+\sum_{n}\left[\delta\theta_{n}^{\rm T}(t)+\delta\theta_{n}^{\rm
R}(t)\right]$ (12) $\displaystyle\delta\theta_{n}^{\rm T}$
$\displaystyle=\frac{B_{n}^{\rm T}p_{0}^{2}\exp(i2\pi nf_{\rm
K}t)}{GM\left(n^{6}-5n^{4}+4n^{2}\right)}\Big{(}-12+7n^{2}-n^{4}$
$\displaystyle+2n^{2}(2+n^{2})\mathcal{C}_{2}-i6n^{3}\mathcal{S}_{2})\Big{)}$
(13) $\displaystyle\delta\theta_{n}^{\rm R}$ $\displaystyle=-\frac{B_{n}^{\rm
R}p_{0}^{2}\exp(i2\pi nf_{\rm
K}t)}{GM\left(n^{6}-5n^{4}+4n^{2}\right)}n\Big{(}i(-4+n^{2})4$
$\displaystyle+i3n^{2}\mathcal{C}_{2}-n(2+n^{2})\mathcal{S}_{2})\Big{)}$ (14)

where $p_{0}$ is the semi-latus rectum of the unperturbed orbit,
$\mathcal{S}_{m}\equiv\sin(2\pi mf_{\rm K})$ and
$\mathcal{C}_{m}\equiv\cos(2\pi mf_{\rm K})$, and only the real part of the
perturbation is relevant. Here the results are shown for the case
$n\neq[0,1,2]$, which are discussed below.

The remaining orbital elements take the schematic form:

$\displaystyle p(t)$ $\displaystyle=p_{0}+\delta p(t)$ (15) $\displaystyle
e(t)$ $\displaystyle=\delta e(t)$ (16) $\displaystyle\omega(t)$
$\displaystyle=\delta\omega(t)$ (17) $\displaystyle\nu(t)$ $\displaystyle=2\pi
f_{\rm K}t+\delta\nu(t),$ (18)

where the perturbations (denoted by $\delta$) are purely oscillatory and
indeed do not cause any secular evolution of the binary over timescales longer
than the forcing. We also report the explicit solution for the binary’s
separation, which oscillates about its unperturbed value:

$\displaystyle r(t)$ $\displaystyle=p_{0}+\sum_{n}\Bigg{[}\frac{B^{\rm
R}_{n}p_{0}^{3}\exp(i2\pi nf_{\rm K}t)}{GM(1-n^{2})}$
$\displaystyle-\frac{i2B^{\rm T}_{n}p_{0}^{3}\exp(i2\pi nf_{\rm
K}t)}{GMn(1-n^{2})}\Bigg{]}.$ (19)

We find, as expected, that oscillatory forces cause a corresponding behaviour
in the orbital elements with a periodicity comparable to the forcing. A
visualisation of the resulting binary motion can be seen in Figure 1. In
essence, Eqs. 12 and II.1 state that a periodic, in-plane perturbation with an
amplitude $B$ and with a periodicity of $nf_{\rm K}$ will produce variations
in the orbital phase of typical order $Bp_{0}^{2}/(GM)$ at the different
harmonics $(2-n,n,2+n)$. The amplitude of high-frequency fluctuations is
suppressed as $1/n^{2}$, due to the shorter available time to apply an
acceleration. We can summarise the results schematically as follows:

Force $\to$ Phase $\big{\\{}B,[n]f_{\rm K}\big{\\}}$ $\to$
$\Big{\\{}\sim\frac{Bp_{0}^{2}}{GM},[2-n,n,2+n]f_{\rm K}\Big{\\}}$

In Fig. 2 we show the dependence of the true longitude perturbation
coefficients on the periodicity $n$. We find that, when neglecting some
special values of $n$ as discussed below, the coefficients for tangential and
radial perturbations can be well fit by the simpler expressions $3/n^{2}$ and
$1/(n+n^{2})$, respectively.

The cases of $n=0$, $n=1$, and $n=2$ must be treated separately in order to
avoid divergences. They can be understood as the binary’s orbital frequency
entering into resonance with a corresponding component of the perturbing
force. In general, resonances are an important factor in determining the
evolution of perturbed binaries, and may even determine the dynamical
stability of triple systems [30, 86]. Additional PN contributions to the
binary’s motion, such as the general relativistic apsidal precession of the
binary, may also resonate with the outer forcing, either constructively or
destructively. Finally, we further note that the interaction of periodic
forces with eccentricity will give rise to further resonances at higher epi-
cyclical harmonics. For the remainder of this analysis, we will not model the
effect of resonances, leaving them for future investigation.

Figure 2: We plot the normalised amplitudes of the true longitude perturbation
coefficients (Eqs. 12). Both radial and tangential perturbations are split
into several components, at frequencies of $\left[n-2,n,n+2\right]f_{\rm K}$
and in sinusoidal and co-sinusoidal parts. The dominant perturbations are
tangential at the frequency $nf_{\rm K}$, and their scaling with the parameter
$n$ can be well approximated by $3/n^{2}$ when disregarding resonances. We do
not fit the $n\pm 2$ perturbations (dashed lines) as they are sub-dominant.
Note that tangential perturbations always dominate unless $n>>1$.

###  II.2 Vacuum and perturbed waveform models

The GW emission of a binary depends on the physical parameters describing its
position in space, i.e. its radius $r$, the true longitude $\theta$ as well as
the inclination angles with respect to the observer. In the Newtonian limit,
the relation between the source and the gravitational radiation far away is
given by the quadrupole formula:

$\displaystyle h_{ij}(t)=\frac{2G}{c^{4}D}\ddot{\mathcal{I}}_{ij}(t_{\rm r}),$
(20)

where $\mathcal{I}$ is the mass quadrupole tensor in its traceless form, $D$
is the source’s luminosity distance and $t_{\rm r}$ is the retarded time. In
our case, additional perturbations to the quadrupole tensor arise from from
both the radius and the true anomaly (Eqs. 12 and II.1). We insert the
explicit solutions to the orbital elements in Eq. 20 and proceed to calculate
the components of the GW strain. After some intensive but straightforward
algebraic computations we find that the perturbed gravitational waveform can
be expressed in the following way:

$\displaystyle h\sim h_{\rm N}(\theta^{\prime})+\delta h(\delta
r,\delta\dot{r},\delta\ddot{r},\delta\dot{\theta},\delta\ddot{\theta}),$ (21)

where we defined the perturbed variable $\theta^{\prime}=\theta+\delta\theta$
as in Eq. 12. Here $h_{\rm N}$ is the familiar Newtonian quadrupolar waveform,
while $\delta h$ are additional perturbations to the radiated GW. The latter
have been studied in detail in [156], and only contain terms explicitly
proportional to $\delta r$, $\delta\dot{r}$, $\delta\ddot{r}$,
$\delta\dot{\theta}$ and $\delta\ddot{\theta}$. The interpretation of this
result is that the perturbed waveform consists of two separate parts: First, a
Newtonian waveform evaluated at a perturbed phase. Second, a series of
additional low-amplitude harmonics dependent on the perturbations in the
orbital elements and their derivatives, which arise directly from variations
in the binary’s lever arm $r+\delta r$ or radial and tangential radial
velocities.

All the components in Eq. 21 carry information regarding the forces acting on
the binary and could in principle be included in a fully general perturbed
waveform model. Additionally, the mass distributions responsible for producing
the forces could themselves produce additional GW radiation. For the purposes
of this work, we focus solely on the first part of Eq. 21, led by the standard
intuition that phase perturbations are more easily detectable than amplitude
perturbations of the same order [34]. Furthermore, we will see that
perturbations of this type already contain an incredible amount of information
regarding the forces acting on the binary. Our goal is thus to construct a
waveform model that can account for such periodic phase perturbations over the
full observation of a chirping source. We base our methodology on the seminal
work in [34] and construct a Fourier space GW template using the stationary
phase approximation, which is justified as long as the perturbations are
sufficiently small. In the stationary phase approximation, a Newtonian
waveform with a phase $\phi_{\rm GW}$ is given by:

$\displaystyle\tilde{h}(f)=\frac{Q}{D}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{5/6}f^{-7/6}\exp\left[i(-\phi_{\rm
GW}-\pi/4)\right],$ (22)

where $f$ is the observed GW frequency, $\mathcal{M}$ is the red-shifted chirp
mass, $D$ is the luminosity distance, $Q$ is a geometric pre-factor that
accounts for projections of the wave onto a given detector and the tilde
denotes a Fourier transform. In our case then, the GW phase is simply the sum
of the vacuum phase and the phase perturbations:

$\displaystyle\phi_{\rm GW}$ $\displaystyle=\phi_{\rm
vac}+2\sum_{n}\left[\delta\theta_{n}^{\rm T}(t)+\delta\theta_{n}^{\rm
R}(t)\right]$ (23) $\displaystyle\phi_{\rm vac}$ $\displaystyle=\phi_{\rm
c}-2\pi ft+2\left(8\pi\mathcal{M}f\right)^{-5/3}(+\text{ PN corrections}),$
(24)

where $\phi_{c}$ is the initial phase and the factor 2 multiplying the sum
arises from the relation $\phi_{\rm GW}=2\theta$. Note that now it is
important to make explicit the radial dependence of the force Fourier
coefficients also in terms of the GW frequency, due to the frequency
dependence of Eqs. 12 and II.1:

$\displaystyle B^{\rm{T}}_{n}$ $\displaystyle\to B^{\rm{T}}_{n}(p_{0})\text{
or }B^{\rm{T}}_{n}(f)$ (25) $\displaystyle B^{\rm{R}}_{n}$ $\displaystyle\to
B^{\rm{R}}_{n}(p_{0})\text{ or }B^{\rm{R}}_{n}(f).$ (26)

Higher accuracy in terms of post-Newtonian terms may also easily be achieved
by applying the perturbed phase appropriately to the required higher order GW
modes 555One may also easily add PN corrections to the binary’s phase without
adding the corresponding higher-order GW multipole terms, as is commonly done
in analytical waveform models. Adapting the notation in [19] we have:

$\displaystyle
h(f)=\frac{Q}{D}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{5/6}f^{-7/6}\times\sum_{lm}\mathcal{H}^{lm}\exp\left[-im\psi\right],$
(27)

where the mode coefficients $\mathcal{H}^{lm}$ are known to very high post-
Newtonian orders, and:

$\displaystyle\psi$ $\displaystyle=-\pi/4+\phi_{\rm c}-2\pi
ft+2\left(8\pi\mathcal{M}f\right)^{-5/3}$ (28) $\displaystyle+\text{ PN
corrections}+2\sum_{n}\left[\delta\theta_{n}^{\rm T}(t)+\delta\theta_{n}^{\rm
R}(t)\right].$ (29)

Eqs. 27 and 28 represent a post-Newtonian, quasi-circular waveform model that
can describe the effect of arbitrary time-varying forces on the binary’s
phase. The number of additional free parameters with respect to a vacuum
waveform is limited to the number of free parameters required to specify the
perturbations on the binary. For clarity’s sake, let us summarise again the
assumptions that went into the construction of the waveform model:

  * •

We focus on circular binaries subject to in-plane perturbative forces.

  * •

The interaction between the relativistic binary and the perturbative forces is
treated at lowest order, i.e. Newtonian.

  * •

Any secular evolution of the orbital elements caused by perturbative forces is
removed by construction, including that arising from resonances.

  * •

Perturbations to the GW emission are similarly treated at lowest order and
only phase perturbations are considered. The additional low amplitude
harmonics were studied in [156].

  * •

GW emission resulting from the stress energy tensor of the environment causing
the perturbation is neglected.

Given these assumptions, we propose our perturbed waveform model as a proxy
for more sophisticated templates that can account for both strong gravity
effects and environmental perturbations in a fully relativistic way. We
highlight the exciting prospect of measuring forces directly from, e.g.,
hydrodynamical and GRMHD simulations of gas-embedded binaries or precise
integrations of few body dynamical scatterings, including them in Eqs. 27 and
28, and performing parameter estimation directly on such numerical results.

##  III SNR criteria

###  III.1 Baseline model and Newtonian SNR estimates.

From this point on, we will omit to specify both the summation and the
distinction between radial and tangential forces, and simply refer to a
generic phase perturbation caused by a force of magnitude $B$ with a
periodicity of frequency $nf_{\rm K}$. From Eqs. 12, we know that the phase
perturbation will consist of three distinct oscillations at the characteristic
frequencies $[2-n,n,2+n]f_{\rm K}$, and a characteristic dimensionless
amplitude$\sim B_{n}p_{0}^{2}/(GM)$. We can further simplify the waveform
model presented in Section II by assuming that the central component of the
phase perturbation triplet $[2-n,n,2+n]f_{\rm K}$ dominates, as suggested by
Fig. 2. Furthermore, we assume that the radial dependence of the Fourier
coefficients follows a power law.

Under these assumptions, we can then construct a simplified version of the
resulting GW-phase perturbation:

$\displaystyle\phi^{\rm ess}(f)$
$\displaystyle=\phi_{\rm{vac}}(f)+B^{\prime}\left(\frac{f}{f_{\rm{in}}}\right)^{\beta^{\prime}}\cos(nft(f)+\phi_{0}^{\rm
B}),$ (30)

where $B^{\prime}$ is a dimensionless amplitude which we scale at the
frequency $f_{\rm in}$ at which the source becomes observable to a particular
detector. Eq. 30 still captures the qualitative effect of the complete true
longitude perturbations on the actual GW phase, and is more easily treated for
order of magnitude estimates. We therefore refer to it as a baseline model for
the perturbations treated in this work. Note that, for clarity’s sake, here
and from now on, the primed quantities are always specifically related to GWs,
while the non-primed quantities relate to the forces acting on the binary. If
the radial dependence of the perturbing force follows a power law with an
exponent $\beta$, we obtain the following relations between the forces and the
phase perturbations:

$\displaystyle\beta^{\prime}$ $\displaystyle=\beta-4/3,$ (31) $\displaystyle
B^{\prime}$
$\displaystyle\approx\frac{6B_{n}}{n^{2}(GM)^{1/3}}\pi^{-4/3}f_{\rm
in}^{-4/3},$ (32)

where in the last step we replace the polynomial in $n$ from Eqs. 12 with the
simplified fit of $3/n^{2}$ (see also Fig. 2), and here $\beta^{\prime}$
describes the frequency dependance.

For the purposes of this work, we will compute simple SNR estimates using the
standard definition of the noise-weighted scalar product for two waveforms
$h_{1}$ and $h_{2}$:

$\displaystyle(h_{1},h_{2})=2\int_{0}^{\infty}\frac{\tilde{h}_{1}^{*}\tilde{h}_{2}+\tilde{h}_{1}\tilde{h}_{2}^{*}}{S_{n}}{\rm{d}}f^{\prime},$
(33)

where $S_{n}$ is the detector’s noise power spectral density. The signal to
noise of a GW is then defined simply as:

$\displaystyle\text{SNR}^{2}=(h,h).$ (34)

We also use a simple prescription to assess the detectability of GW phase
perturbations, which is commonly used in works that treat environmental
effects. We require the noise-weighted scalar product of the GW perturbation
to exceed a certain threshold which we leave unspecified, but is commonly
chosen to be $\sim\mathcal{O}(10)$:

$\displaystyle\delta\text{SNR}^{2}=(\delta h,\delta
h)=8\int\frac{\tilde{h}^{2}}{S_{n}}\left[1-\cos\left(\delta\phi\right)\right]\,{\rm{d}}f^{\prime}.$
(35)

where the prefactor 8 arises from the normalisation of the one sided power
spectrum [101]. Note that the latter equation is often only a rough estimate
of the detectability of a waveform perturbation. In particular, the commonly
chosen criterion $\delta$SNR ¿8 does not account for degeneracies and
correlations between the binary vacuum parameters and the perturbations under
scrutiny. Thus, Eq. 35 often overestimates the detectability of secular de-
phasing, which is strongly degenerate with both the binary’s initial phase as
well as other parameters such as effective spin666This particular degeneracy
is in fact responsible for the mismatch between the Bayesian analysis
performed in this work and recent works [see, e.g., 61] on gas de-phasing
based on the pycbc package, as the latter treats the initial phase as an
extrinsic parameter., as seen in Fig. 7. For this reason, we will test the
simple detectability criteria with several MCMC parameter recovery tests. From
here on, we will refer to the results of Eq. 35 as a ”$\delta$SNR” estimate.

Let us also briefly discuss secular de-phasing. In binaries perturbed by
external forces, secular de-phasing is caused by the exchange of energy and
angular momentum between the binary and a separate reservoir. It is mediated
by a residual tangential force $A_{0}^{\rm T}(r)$, from which we remove any
time variation by performing a time average. We can find the resulting secular
drift in binary separation by using the Lagrange planetary equations:

$\displaystyle\dot{p}=2A_{0}^{\rm T}(p_{0})\sqrt{\frac{p_{0}^{3}}{GM}}.$ (36)

When it comes to GW observables, the additional drift in orbital separation
modifies the source’s chirp evolution:

$\displaystyle\frac{\dot{p}}{p}=-\frac{3}{2}\frac{\dot{f}}{f}=-\frac{3}{2}\frac{\dot{f}_{\rm{GW}}+\dot{f}_{\rm
A}}{f},$ (37)

where:

$\displaystyle\dot{f}_{\rm GW}$
$\displaystyle=\frac{96}{5}\pi^{8/3}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{5/3}f^{11/3},$
(38) $\displaystyle\dot{f}_{\rm A}$ $\displaystyle=-\frac{4}{3\pi}A_{0}^{\rm
T}(f)\left(\frac{f^{2}}{GM\pi^{2}}\right)^{1/3}.$ (39)

The accumulated GW phase over the binary chirp is then given by:

$\displaystyle\phi=2\pi\int\frac{f}{\dot{f}_{\rm{GW}}+\dot{f}_{\rm
gas}}\,{\rm{d}}f\approx\phi_{\rm vac}-\int\frac{f\dot{f}_{\rm
A}}{\dot{f}_{\rm{GW}}^{2}}{\rm{d}}f,$ (40)

where:

$\displaystyle\phi_{\rm
vac}\approx\frac{1}{32\pi^{8/3}}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{-5/3}f_{\rm
in}^{-5/3}$ (41)

The cumulative effect of a tangential force thus causes a secular de-phasing
of the form:

$\displaystyle\phi^{\rm d}(f)$
$\displaystyle=\phi_{\rm{vac}}(f)+A^{\prime}\left(\frac{f}{f_{\rm{in}}}\right)^{\alpha^{\prime}},$
(42)

where $A^{\prime}$ is a dimensionless coefficient and we assumed that the
radial dependence of the secular force follows a power law with an exponent
$\alpha$. The general relation between the force and the secular de-phasing
is:

$\displaystyle\alpha^{\prime}$ $\displaystyle=\alpha-14/3$ (43) $\displaystyle
A^{\prime}$ $\displaystyle=\frac{A_{0}^{\rm T}}{(GM)^{1/3}}\frac{25f_{\rm
in}^{-14/3}}{2304(14-3\alpha)\pi^{7}}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{-10/3}.$
(44)

Comparing this result with Eq. 31, we find that the amplitude of secular de-
phasing will typically accumulate in the early inspiral and decay extremely
rapidly as the source chirps.

###  III.2 Bayesian parameter recovery with post-Newtonian waveforms.

In this paper we perform several numerical parameter recovery tests to
validate the simpler analytical prescriptions of the previous sections. We
adopt a vacuum waveform model of the type described by Eqs. 27 and 28. We
include the phase evolution of a non-spinning binary up to third PN order by
using the formulas from [19]. Similarly, we also account for all higher order
modes beyond quadrupolar GW emission up to order $x^{5/2}$, where $x$ is the
post-Newtonian parameter. Following [34], we also include the leading order
effect of spin through the effective parameter $\mathcal{B}$, which enters the
phase at first-and-a-half post-Newtonian order and represents a projection of
the orbital angular momentum and the individual binary component spins. In
this work, we do not include any spin evolution term nor the effects of spin-
spin coupling. Finally, we include the effect of periodic dephasing by means
of the baseline model detailed above, for a single periodicity. While this
suffices to obtain correct detectability estimates, we stress that the latter
cannot replace the full formulas for true longitude perturbations in actual
applications.

We limit our analysis to the binary’s intrinsic parameters and initial phase,
while performing a sky average of a source’s position on the sky and
orientation. Furthermore, we assume that the binary coalescence time may be
measured with high precision through the merger and ring-down portions of the
GW signal, which we do not model here. Therefore, our post-Newtonian waveform
is completely specified by the following nine independent parameters:

$\displaystyle\text{Binary}=$ $\displaystyle\begin{cases}z;&\text{Redshift}\\\
\mathcal{M};&\text{Chirp mass}\\\ \mu;&\text{Reduced mass}\\\ \phi_{\rm
c};&\text{Initial binary phase}\\\ \mathcal{B};&\text{Effective
spin}\end{cases}$ (45) $\displaystyle\text{Perturbation}=$
$\displaystyle\begin{cases}\log_{10}(B^{\prime});&\text{Initial amplitude}\\\
\beta^{\prime};&\text{Frequency scaling}\\\ n;&\text{Periodicity}\\\
\phi_{0}^{B};&\text{Init. perturbation phase}\end{cases}$ (46)

We use Monte-Carlo-Markhov-Chain (MCMC) methods to perform parameter
estimation studies numerically. The goal of MCMC based parameter estimation is
to sample the posterior distributions of a set of waveform parameters $\Theta$
given a likelihood function $\mathcal{L}$ of the form:

$\displaystyle\mathcal{L}\left(\Theta\right)\propto\exp\left[-\big{(}h(\Theta)-h(\Theta_{\rm
GT}),h(\Theta)-h(\Theta_{\rm GT})\big{)}\right]$ (47)

where $\Theta_{\rm GT}$ is a vector of parameters for the injected signal. The
detector noise is treated as being stationary, representing the average value
for the noise power spectral density over many realisations. We use the affine
invariant sampler emcee [56] to perform the numerical tests, typically running
32 parallel walkers for approximately 30’000 to 50’000 steps, where the
typical auto-correlation time of the walkers is $\sim 200$. We initialise the
walkers in the vicinity of the injected true values in order to speed up
convergence and to avoid having to specify priors for the vacuum binary
parameters (which is appropriate only for high SNR sources). The priors for
the perturbation parameters are simply chosen for convenience and do not
influence our results unless stated otherwise.

##  IV Detectability of Periodic Phase perturbations

###  IV.1 Detectability for monochromatic sources

Here we use the baseline GW phase model constructed in Sec III.1 to extract
some qualitative insight regarding the detectability of periodic
perturbations. We start from the definition of $\delta$SNR:

$\displaystyle\delta\text{SNR}^{2}=(\delta h,\delta
h)=8\int\frac{\tilde{h}_{\rm
vac}^{2}}{S_{n}}\left[1-\cos\left(\delta\phi\right)\right]\,{\rm{d}}f^{\prime},$
(48)

and want to evaluate it for:

$\displaystyle\delta\phi=B^{\prime}\left(\frac{f}{f_{\rm{in}}}\right)^{\beta^{\prime}}\cos(nft(f)+\phi_{0}^{\rm
B}),$ (49)

representing a periodic perturbation with a typical amplitude $B^{\prime}$ at
the GW frequency $f_{\rm in}$, a frequency scaling given by $\beta^{\prime}$,
a typical periodicity described by the frequency $nf_{\rm K}$ and an inital
phase $\phi_{0}^{\rm B}$. Evaluating Eq. 48 in the limit of small $B^{\prime}$
and also for monochromatic sources, results in the following expression:

$\displaystyle\delta\text{SNR}^{2}$ $\displaystyle\approx\frac{\tilde{h}_{\rm
vac}^{2}\dot{f}}{S_{\rm n}}B^{\prime
2}\left(\frac{f}{f_{\rm{in}}}\right)^{2\beta^{\prime}}\int^{T}\cos(nft+\phi_{0}^{\rm
B})^{2}\,{\rm d}t.$ (50)

where we replaced the integration in frequency with an integration over the
observation time $T$. Here we can already see that, as $T$ is much longer than
a wave period, the dependence on the parameters $n$ and $\phi_{0}^{\rm B}$
should be washed away by the integration. This is also true for chirping
sources, as long as the number of cycles at a given frequency is larger than
$\sim n^{-1}$, as shown in Fig. 3. The average of the integral then simply
evaluates to $\sim 1/2\times T$ and reduces to the convenient form:

$\displaystyle\delta\text{SNR}^{2}\approx\frac{1}{2}\text{SNR}(f)^{2}B^{\prime}(f)^{2},$
(51)

where we used the definition of the signal-to-noise ratio for a monochromatic
unperturbed source [see, e.g., 19, 83, 101]. We can also easily extrapolate
the results of performing the integral in Eq. 50 for a full spectrum of
perturbations. Re-introducing the subscript $n$ to denote the various periodic
force components, we find that the total SNR grows in quadrature:

$\displaystyle\delta{\text{SNR}}_{\rm
tot}=\text{SNR}\times\sqrt{\frac{1}{2}\sum_{n}(B^{\prime}_{n})^{2}}\approx\text{SNR}\times
B^{\prime}\sqrt{N/2},$ (52)

where in the latter approximation we define $N$ as an effective number of
periodic perturbations with amplitude $B^{\prime}$. The latter factor can
significantly boost the detectability of periodic force perturbations from EE
that naturally produce perturbations at many frequencies, such as, e.g., the
complex gas hydrodynamics around merging binaries, see Sec. V. Finally, we
report the threshold acceleration required for detection of a single periodic
perturbation:

$\displaystyle B_{n}^{\rm det}$ $\displaystyle>A_{\rm
G}\frac{n^{2}\sqrt{2}}{3\times\text{SNR}}$ (53)
$\displaystyle=\frac{n^{2}}{\text{SNR}}\frac{\pi^{4/3}\sqrt{2}}{6}(GM)^{1/3}f^{4/3},$
(54)

where $A_{\rm G}$ is the binary’s gravitational acceleration. This equation
states that individual periodic force components are detectable if they are of
the order of $A_{\rm G}/n^{2}/$SNR, which strongly prefers sub-orbital
perturbations.

Figure 3: We show the variations in the $\delta$SNR caused by varying the
parameters $n$ and $\phi_{0}^{\rm B}$ in Eq. 49 for a typical chirping source
that completes $N_{\rm c}$ cycles in band. As expected, we find that chirping
sources that complete many cycles tend to wash out the importance of the
perturbation’s initial phase. The latter’s variation is completely negligible
as long as the source typically completes more than $\sim 10/n$ cycles at each
corresponding frequency bin.

####  IV.1.1 Comparison with secular dephasing

Eq. 50 suggests that for monochromatic sources, phase oscillations of
approximately $\delta\phi\sim{\text{few }}\times 1/$SNR should be detectable,
which seems comparable to the analogous requirement for secular de-phasing, in
the limit where degeneracies may be ignored. The comparison is however not
trivial, as in the latter case dephasing arises as a consequence of the slow,
coherent accumulation of the perturbation over the entire observation. For
almost monochromatic sources, the GW phase accumulates as:

$\displaystyle\phi$
$\displaystyle=\int_{0}^{T}\left[f+\dot{f}(f)t\right]\,\rm{d}t.$ (55)

where $T$ is the observation time. In case of a secular perturbation, the
total accumulated dephasing $\Delta\phi$ is given by:

$\displaystyle\Delta\phi=\frac{1}{2}\dot{f}_{\rm A}T^{2}$ (56)

Evaluating Eq. 35, the expected $\delta$SNR produced by secular de-phasing
amounts to:

$\displaystyle\delta\text{SNR}^{2}$ $\displaystyle\approx\frac{\tilde{h}_{\rm
vac}^{2}\dot{f}}{S_{\rm n}}\dot{f}_{\rm gas}^{2}\int^{T}t^{4}\,{\rm
d}t=\frac{4}{5}\text{SNR}(f)^{2}\Delta\phi^{2}.$ (57)

The previous equation speaks to the high potential detectability of secular
dephasing. It suggests that, indeed, cumulative shifts of order few
$1/\text{SNR}$ may be detectable. However, this conclusion ignores the fact
that secular dephasing, by its own nature of being slowly accumulated, can be
degenerate with many binary parameters that act in a similar way. In
particular, it is highly degenerate with the binary’s initial phase, along
with any post-Newtonian modifications to the binary’s chirp behaviour with a
similar magnitude or frequency scaling. A more conservative requirement for
the detectability of secular dephasing is that it should exceed a full cycle,
i.e. $2\pi$, while many works have also highlighted cases where several cycles
of dephasing are necessary to perform parameter estimation [see, e.g., 124,
72]. Here we use $2\pi$ as a detectability criterion, finding that the
threshold acceleration for secular de-phasing amounts to:

$\displaystyle A_{0}^{{}^{\rm det}}=3\pi^{5/3}(GM)^{1/3}f^{-2/3}T^{2}.$ (58)

The ratio between the detection thresholds for periodic and secular dephasing
yields some insight:

$\displaystyle\frac{B_{n}^{{}^{\rm det}}}{A_{0}^{{}^{\rm
det}}}=\frac{1}{{\text{SNR}}}\frac{\sqrt{2}4n^{2}}{9\pi^{1/3}}(fT)^{2}\sim\frac{N_{\rm
c}^{2}}{\text{SNR}},$ (59)

where $N_{\rm{c}}$ is the number of cycles at the frequency $f$.
Interestingly, this suggests that periodic perturbations may be the dominant
EE for high SNR sources that complete comparatively fewer cycles in the
sensitivity band of their respective detectors. As we will see in Sec. V, this
may be realised in several astrophysical scenarios. The euation above also
suggests that, in sources where both secular and periodic forces are present,
the latter must be significantly higher in amplitude to be detected at the
same $\delta$SNR. However a more precise statement on detectability requires a
careful Bayesian analysis. Indeed, as we will see in IV.4, these naive
expectations are often misleading.

###  IV.2 Detectability for chirping LISA sources

When stepping beyond the monochromatic case, it is necessary to account for
the full noise power spectral density of a detector. Here we produce some
simple detectability criteria for the recently adopted, space based
gravitational wave detector LISA [4, 33], but note that similar estimates
apply to any similar milli-Hz GW instrument [e.g. TianQin 87].

As a proxy for LISA’s noise power spectral density, we use the convenient fit
from [101] and evaluate Eq. 35 for a range of perturbation parameters
$B^{\prime}$ and $\beta^{\prime}$. We repeat the process for an array of
source total masses, mass-ratios and redshifts. The $\delta$SNR results for
the baseline perturbations, for typical massive BHB sources, are visualised as
a contour plot in the top panel of Fig. 4. The $\delta$SNR contours have a
strong dependence on the parameter $\beta^{\prime}$, and generally align with
the expectations set from the monochromatic case, i.e. they scale inversely
with the source’s total SNR and linearly with the perturbation amplitude
$B^{\prime}$. Here we highlight some results for the case of a constant
amplitude for the phase perturbation, $\beta^{\prime}=0$, and the case for a
constant amplitude in radius or frequency of the perturbing force,
$\beta^{\prime}=-4/3$. For the former, the perturbations are detectable for a
wide range of sources whenever $B^{\prime}\sim 10^{-3}$ to $10^{-2}$. For the
latter, detectability requires $B^{\prime}\sim 0.1$. The following formulae
provide a simple fit to the results of Fig. 4 at $z=1$:

$\displaystyle\frac{\delta{\text{SNR}}}{B^{\prime}}\approx\sqrt{\frac{2\mu}{M}}\begin{cases}\left(10^{1.2\beta^{\prime}+2.5}+5.9\right)&{\text{for
}}M=10^{4}\,\rm{M}_{\odot}\\\
\left(10^{1.8\beta^{\prime}+3.5}+3.0\right)&{\text{for
}}M=10^{5}\,\rm{M}_{\odot}\\\
\left(10^{1.8\beta^{\prime}+3.8}+3.0\right)&{\text{for
}}M=10^{6}\,\rm{M}_{\odot}\\\
\left(10^{1.2\beta^{\prime}+2.9}+1.0\right)&{\text{for
}}M=10^{7}\,\rm{M}_{\odot}\end{cases}$ (60)

where we separated typical LISA binary mergers into several mass decades.

We confirm these estimates with a series of numerical tests, also partially
shown in the bottom panel of Fig. 4. We take a representative source with
$M=2.2\times 10^{6}$ M⊙ and define three detectability threshold lines at
$\delta$SNR=1, 8 and 64 according to Eq. 35. We then perform several MCMC
parameter recovery experiments with post-Newtonian waveforms, in which all
parameters detailed in Eq. 45 are sampled freely. We plot the standard
deviations of the marginalised posterior distributions for $\beta^{\prime}$
and $B^{\prime}$ as errorbars in the bottom panel of Fig. 4 for a series of
injected values of $B^{\prime}$ and $\beta^{\prime}$ that lie along the
corresponding $\delta$SNR lines. We find that all injected perturbations can
be recovered, provided that they lie above a $\delta$SNR of a few. In general
we find the following qualitative correspondence between the $\delta$SNR and
the MCMC based parameter recovery:

$\delta$SNR MCMC Results $\displaystyle\sim 1$
$\displaystyle\implies\text{"Null hypothesis" excluded at 1$\sigma$.}$
$\displaystyle\sim 8$ $\displaystyle\implies\text{1$\sigma$ uncertainty < 0.5
for $\log_{10}(B^{\prime})$ and $\beta^{\prime}$}$ $\displaystyle\sim 64$
$\displaystyle\implies\text{1$\sigma$ uncertainty <0.05 for
$\log_{10}(B^{\prime})$ and $\beta^{\prime}$}.$

Here the ”null hypothesis” represents the complete absence of a perturbation
within the waveform. We have tested the results showcased in Fig. 4 for a
range of total masses between $M=10^{3}$ M⊙ and $M=10^{8}$ M⊙, with mass-
ratios from $q=0.1$ to $q=1$ and for a few redshifts between $z=0.5$ and
$z=3$. Since the aforementioned qualitative trends seem to hold without fail,
we conclude that our $\delta$SNR based detection criteria are robust, and can
provide a good estimate for detectability even when accounting for parameter
degeneracies and post-Newtonian modelling for the unperturbed vacuum
waveforms. We emphasise that this is in stark contrast to secular de-phasing,
as further discussed below.

Figure 4: Top panel: The detectability contours for a periodic dephasing with
amplitude $B^{\prime}$ and frequency scaling $\beta^{\prime}$, for a wide
range of LISA sources with total mass $M$, mass-ratio $q=1$ at redhsift $z=1$.
The dashed lines highlight the cases of a constant perturbation amplitude
($\beta^{\prime}=0$) or a constant force ($\beta^{\prime}=-4/3$) throughout
the chirp. Bottom panel: The recovered marginalised posteriors of
$\beta^{\prime}$ and $B^{\prime}$ from several MCMC parameter estimation
tests, with injected perturbations that lie along the $\delta$SNR = 8 (blue)
and $\delta$SNR = 64 (green) threshold lines. The grey dots represent the
density of samples of an MCMC run testing the ”null hypothesis”, i.e. the
absence of a perturbation. The latter is excluded at $1\sigma$ for
perturbations with $\delta$SNR above 1, which is represented by the red dashed
line.

###  IV.3 Detectability for future ground based detectors.

Here we perform the same analysis as above, while adopting the prospective
sensitivity curve of the Einstein Telescope (ET) as a proxy for the next
generation of ground based GW detectors [84]. The implications of the latter
on the astrophysics of stellar-mass black holes will be further discussed in
Sec. V.

The results of our analysis are showcased in Fig. 5, and similar
considerations apply here as in the previous section. The highest
detectability is achieved for sources with total masses of $\sim$few$\times
10^{2}$ M⊙, which maximise the SNR of the detection. This is comparable to the
highest observed masses in the current LIGO catalogue of sources, GW190521
[2]. Interestingly, detecting EE in latter type of sources may also be
possible with LISA, during the very early inspiral stages [120, 63, 136, 118].
We highlight that with the given ET sensitivity, the detectability of periodic
perturbations in chirping sources is plausible for a much larger range of
$\beta^{\prime}$ than for LISA sources, reaching down to
$\beta^{\prime}\sim-3$ for high SNR sources, due to the different form of the
noise power spectral density. Similarly to the previous section, we provide
fits to the resulting $\delta$SNR, which may be used as a simple estimate:

$\displaystyle\frac{\delta{\text{SNR}}}{B^{\prime}}\approx\sqrt{\frac{2\mu}{M}}\begin{cases}\left(10^{1.4\beta^{\prime}+1.4}+\beta^{\prime}+5\right)&{\text{for
}}M=3\,\rm{M}_{\odot}\\\
\left(10^{0.8\beta^{\prime}+2.2}+1.0\right)&{\text{for
}}M=30\,\rm{M}_{\odot}\\\
\left(10^{0.6\beta^{\prime}+2.1}+1.7\right)&{\text{for
}}M=300\,\rm{M}_{\odot}\end{cases}$ (61)

Figure 5: Same as Fig. 4, using the noise power spectral density of ET [84] as
a proxy for future ground based detectors.

###  IV.4 Test case: Degeneracies in secular vs periodic dephasing

Figure 6: Full showcase of the marginalised posteriors of a high SNR post-
Newtonian waveform. The injected periodic dephasing has a Newtonian
$\delta$SNR of approximately 40, and can be recovered very precisely. Note
that the frequency of the perturbation can be recovered with a precision of
0.1$\%$. Crucially, no correlations between the periodic perturbation and the
vacuum binary parameters appear. The MCMC is performed with 32 parallel
walkers for 50’000 steps, requiring approximately 45 minutes of compute on a
single core. See Eqs. 45 for the meaning of the nine parameters.

In Fig 6, we showcase the full 9 parameter posteriors (in black) from a
realisation of an MCMC parameter estimation test for a representative high-SNR
LISA source, where we injected a strong periodic de-phasing with
$B^{\prime}=0.1$, $n=1$ and $\beta^{\prime}=-0.5$, resulting in a $\delta$SNR
of approximately 40. All of the two dimensional marginalised parameter
posteriors are well behaved, clearly localised and centered on the true value.
The only correlations visible in the binary parameters are limited to the
familiar chirp mass and redshift degeneracy, as well as the initial phase and
effective spin parameter degeneracy. The only degeneracy that appears in the
recovered perturbation is between the amplitude $B^{\prime}$ and the frequency
scaling $\beta^{\prime}$. Investigating the latter degeneracy in more detail,
we find that by defining two new parameters:

$\displaystyle\left[\beta^{\prime},\log_{10}(B^{\prime})\right]\to\left[\kappa_{1},\kappa_{2}\right]$
(62)

where:

$\displaystyle\beta^{\prime}=-\kappa_{1}\times\log_{10}(B^{\prime})+\kappa_{2}$
(63)

one can completely characterise the correlation in $B^{\prime}$ and
$\beta^{\prime}$. Here we do not complete a full analysis of this, but note
that by performing the parameter estimation on $\kappa_{1}$ and $\kappa_{2}$,
one can determine the latter within approximately $1\sigma$ of the true
injected value already for sources with $\delta$SNR $\sim 1$. Therefore, while
the parameters $\kappa_{1}$ and $\kappa_{2}$ are less informative than
$B^{\prime}$ and $\beta$, they can be used to distinguish the presence of a
perturbation from its absence, and are therefore connected to the red dashed
line shown in Fig. 4 and described in the text.

The test case showcased in Fig. 6 is representative of all $\delta$SNR larger
than a few, where the widths of the Gaussian posteriors roughly scale as
$1/\delta$SNR. This suggests that the additional parameters required to
describe periodic de-phasing are completely orthogonal to the usual vacuum
binary parameters. Crucially, they do not cause any visible degradation in the
recovery of the vacuum waveform.

We compare the aforementioned result with an attempt to recover the same
vacuum binary source while adding a secular de-phasing term, with a
corresponding amplitude $A^{\prime}=0.1$ and frequency scaling
$\alpha^{\prime}=-0.5$. The latter perturbation actually corresponds to
significantly higher $\delta$SNR with respect to the comparative case for a
periodic perturbation, i.e. approximately 60 instead of 40 (according to Eq.
35). In Fig. 7, we showcase the difference in the recovered posteriors between
the two cases.

We highlight some interesting differences: For the injected value of
$B^{\prime}=A^{\prime}=0.1$, the recovered value of the periodic perturbation
is extremely well constrained, yielding $B^{\prime}_{\rm rec}=0.094\pm 0.04$
at $1\sigma$. In contrast, the posterior for the amplitude of the secular
perturbation is highly degenerate and essentially completely unconstrained,
allowing for values that range all the way from $\sim 0$ to $\sim 1$. Indeed,
the recovered value turns out to be $A^{\prime}=0.069$, with $1\sigma$
uncertainties of $-0.062$ and $+0.69$, respectively (though the specifics here
strongly depend on the chosen prior, which in our case is flat in
$\log_{10}(A^{\prime})$ between $10^{-2}$ and $2\pi$). This degradation in the
reconstruction of the binary’s phase is also reflected in the posterior for
the spin parameter. Attempting to recover a secular de-phasing causes the
latter to broaden significantly, and induces a bias in the recovered value of
approximately $2\sigma$ with respect to the injected value. We note that this
is still true for values of the accumulated dephasing that exceed $2\pi$. We
note furthermore, that attempting to recover a secular dephasing degrades the
convergence of the MCMC.

While only based on post-Newtonian waveforms, this test case strongly
indicates that the periodic perturbations of the type presented in this work
are more informative and less prone to degeneracies with respect to secular
dephasing, when it comes to a sophisticated parameter estimation pipeline (see
Sec. VI for an interpretation of this result).

We are aware that a complete validation of this claim would require a more
systematic study of both secular and periodic dephasing, for many sources and
across an array of different waveform models, which is beyond the scope of
this work. Nevertheless, the results above do suggest that waveforms of the
type presented in Eqs. 27 and 28 represent a strong candidate to search for EE
in GW signals without necessarily degrading the estimation of the vacuum
parameters of the source.

Figure 7: Full showcase of the marginalised posteriors when attempting to
recover a periodic (black) and a secular (red) dephasing. The injected
periodic dephasing is characterised by $B^{\prime}=0.1$ and $\beta$ has a
$\delta$SNR of approximately 40. the injected secular dephasing has the same
amplitude and scaling, resulting in a $\delta$SNR of approximately 60.
Nevertheless, the latter cannot be reconstructed confidently as it suffers
from a degeneracy with the binary’s initial phase and spin parameter. The two
MCMCs are performed with 32 parallel walkers for 50’000 steps, requiring
approximately 90 minutes of compute on a single core. See Eqs. 45 for the
meaning of the nine parameters.

##  V Time-varying forces in astrophysical scenarios

In this section we briefly explore a range of scenarios in which binary
sources of GWs will be affected by time-varying forces. While we comment on
the detectability of the resulting perturbations for relevant near-future GW
detectors, the goal is strictly to highlight examples in which the physical
processes at play naturally lead to strong fluctuations. We start by
mentioning generic multipolar modifications to the gravitational potential,
then delve into compact object triple systems as well as compact objects
embedded in accretion discs. The discussion presented in this section is brief
and mostly qualitative, though based on extrapolations from N-body simulations
and hydro-dynamical simulations of gas-embedded binaries rather than purely
analytical arguments.

###  V.1 Generic multipole perturbations

We apply the formalism developed in this work to a very simple, heuristic
description of multi-polar modifications to the gravitational potential. We
imagine the vacuum gravitational field to be perturbed by a series of
additional components, arising either from the presence of background fluid,
field or from modifications to GR. As the binary orbits in the modified
potential, it will experience time varying forces with the following physical
scaling:

$\displaystyle B\sim J_{l}\frac{GM}{r^{2}}\cos\left(2lf_{\rm k}t\right),$ (64)

where we neglect to model any projection or inclination effect, and the
coefficients $J_{l}$ are analogous to the multi-pole moments commonly used in
planetary physics. Here $l=1$ corresponds to a dipole perturbation to the
gravitational potential, $l=2$ to a quadrupole perturbation, and so on. Then,
the magnitude of the resulting periodic phase perturbations simply reads:

$\displaystyle B^{\prime}_{n}\approx\frac{3}{4l^{2}}J_{l},$ (65)

and the $1/r^{2}$ scaling implies:

$\displaystyle\beta^{\prime}=0.$ (66)

Extrapolating from the monochromatic results discussed in section IV.1 we
conclude that by employing waveforms of the type showcased by Eqs. 27 and 28,
one should expect to be able to constrain completely general multipolar
modifications to the gravitational potential (and therefore the spacetime
metric), provided they exceed a certain typical magnitude:

$\displaystyle J_{l}>J_{l}^{\rm det}\approx\frac{l^{2}}{\text{SNR}}.$ (67)

As an example, the presence of a flattened distribution of energy density of
any kind within the binary’s orbit should be detectable uniquely by its
induced quadrupole, provided that the total enclosed mass is of order
$\sim$few/SNR. Similar considerations apply for any other multi-polar
distribution.

###  V.2 Dynamical and stellar triple systems

Triple systems of compact objects are a rich laboratory of gravitational
dynamics, and are being studied in great detail in the context of the
dynamical formation channel for GW sources [98, 17, 116, 82, 104, 105, 105, 8,
94, 7, 122, 111, 112, 134, 103, 81, 67, 113, 141, 139, 66, 21, 139]. It is
known that a close-by tertiary massive body may induce several detectable
imprints on GW emission from the inner binary. In this context, the most
studied EEs are e.g. Doppler modulations that arise as a consequence of Römer
delays or line-of-sight accelerations, as well as the secular influence of
tidal fields throughout the inspiral [153, 20, 109, 69, 102, 27, 100, 147,
128, 50, 137, 126, 150, 76, 133, 115, 115]. Here we neglect those effects,
focusing instead on the instantaneous forces experienced by the inner binary
components, as they orbit in the time-varying gravitational potential induced
by the tertiary.

In the top panel of Fig. 8 we show the time series of the instantaneous torque
experienced by an inner binary due to its tertiary companion. In this specific
case, the inner binary has masses of 25 and 20 $M_{\odot}$, while the outer
companion has a mass of 30 $M_{\odot}$. The inner binary has a semi-major axis
of $2.5\times 10^{-4}$ AU, and emits GWs at 0.1 Hz. The outer binary has a
semi-major axis of $1.5\times 10^{-3}$, yielding a period ratio of
approximately 11. The simulations are performed with the regularized, direct
$N$-body code tsunami [140, 138], which includes post-Newtonian corrections to
the acceleration up to the 3rd order.

As we can see, the torque time series presents several oscillating components
at different harmonics of the inner binary’s orbital frequency. The amplitude
of the oscillations typically exceeds the orbit-averaged value by a factor
$\sim 5$. The two major oscillatory components can by attributed to inner
binary’s and the outer binary’s orbital motions, respectively. A more detailed
view of the various components is presented in the bottom panel of Fig. 8, in
which we express the torque time series as a periodogram in the form of Eq.
II.1. We observe that the torque spectrum is composed by several distinctive
peaks, the most prominent of which occur at $\sim 1.8f_{\rm K}$ and $\sim
0.086f_{\rm K}$ respectively. The exact value and relation of these two peaks
encodes information about the relative orbital parameters of the binary and
the tertiary. In particular, resolving such torque peaks would yield
information regarding the mass, distance and eccentricity of the latter, and
allow to break the perturber mass-distance degeneracy that exists for secular
de-phasing in triples. We note that this specific realisation of a triple
system can be easily scaled to arbitrary masses and separation by
approximating the potential of the tertiary as a tidal field:

$\displaystyle
B\sim\frac{Gm_{3}}{(R+p_{0})^{2}}-\frac{Gm_{3}}{(R-p_{0})^{2}}\approx\frac{4Gm_{3}p_{0}}{R^{3}},$
(68)

where $m_{3}$ is the tertiary’s mass and $R$ is the distance of the binary to
the tertiary. Like other secular tidal effects, time varying oscillations will
therefore be stronger when the perturber is massive and close-by. Such time-
varying forces always accompany any other effect caused by the presence of a
tertiary body.

Figure 8: Top panel: Time series of the simulated torque on a binary of
compact objects with masses $m_{1}$ and $m_{2}$, orbiting in the gravitational
potential of a tertiary companion with mass $m_{3}$ (see Sec. V.2). The outer
orbit is given an eccentricity of $e_{\rm{out}}=0.1$, which is reflected in
the modulation of the torques. The magnitude is normalised by its own orbit-
averaged value, showcasing that strong fluctuations are present. Bottom panel:
The spectrum of force Fourier-components (in the form of Eq. II.1) that
correspond to the simulated torque, normalised by the binary’s own
gravitational acceleration. We highlight three dominant force peaks at
different frequencies. The dashed lines denote the approximate SNR required to
detect individual force peaks, according to the monochromatic criterion (see
Sec. IV and Eq. 48), serving as a demonstration that such peaks may be
detectable for certain configurations. Note that for this realisation the
$\delta$SNR of the total force spectrum is dominated by the highest peaks.

###  V.3 Massive black holes in circumbinary discs

Figure 9: Top panel: Time series of the simulated hydro-dynamical torque on an
equal-mass SMBH binary, orbiting in a viscous accretion disc (see Sec. V.3).
The magnitude is normalised by its own orbit-averaged value, showcasing that
strong fluctuations are present. Bottom panel: The spectrum of force Fourier-
components (in the form of Eq. II.1) that arise from the simulated torque,
normalised by the binary’s own gravitational acceleration. We highlight three
dominant force peaks at different frequencies. The dashed lines denote the
approximate SNR required to detect individual force peaks, according to the
monochromatic criterion (see Sec. IV and Eq. 48). Here the SNR lines are
scaled with the disc’s scale-mass to binary-mass-ratio, $q_{\rm d}$, defined
in Eq. 69. This example serves as a demonstration that such peaks may be
detectable for certain binary-disc configurations. Note that for this
realisation of a torque spectrum, the total $\delta$SNR of the resulting phase
perturbation is significantly larger than what is suggested by the largest
individual peaks.

The loudest and highest-quality GW sources of future space-based detectors
will come from comparable-mass, massive black holes coalescing in the m-Hz
band. Such binaries are expected to form following galaxy mergers, which are
known to efficiently funnel material into the resulting galactic nucleus [14,
125]. Therefore, many SMBH binaries are expected to undergo gas-embedded
accretion phases which may play an important role in shepherding them from
initially large separations down to GW-dominated regimes [15, 89, 65]. Modern
studies have also suggested that circum-binary material can remain coupled to
the binary significantly in to a GW-dominated regime [55, 45, 75, 9, 31] and
that secular gas-induced changes to the binary’s orbital elements could
feasibly be detected by LISA [62, 61, 131]. However, the study of the time-
variability in such forces has largely been omitted.

The dynamics of gas-embedded comparable-mass systems are highly non-trivial
and have primarily been studied numerically [for a recent review see 77]. A
primary focus of such investigations has been the periodicity of the accretion
rate onto the binary (which is not always dominated by the binary’s orbital
frequency) and its applications to electromagnetic searches for compact
massive binaries [e.g. 48, 57]. In [106] a spectral analysis of the torques in
such systems was performed, but nearly all recent analyses of the gas forcing
and resultant change in orbital elements have focused on the time-averaged,
secular evolution of the orbital elements.

In the top panel of Fig. 9 we illustrate a raw time-series of the torques
experienced by an equal-mass binary fixed on a circular orbit over 40 periods.
We see that the resultant torques in a quasi-steady configuration are in fact
highly variable, with features beyond those solely at the binary orbital
frequency. Here the magnitude of the torque shows fluctuations by a factor of
order $\sim 30$ with respect to its orbit-averaged value. The presented data
is taken from a 2D, isothermal-hydrodynamics simulation performed with
Sailfish [c.f. 143]. The setup is exactly that presented in [52] for a thin,
co-planar disc of characteristic scale-height $h/r\approx\mathcal{M}^{-1}=0.1$
and global kinematic viscosity $\nu=\nu_{\rm d}a^{2}\Omega_{b}$ with $\nu_{\rm
d}=10^{-3}$ (roughly equivalent to a turbulo-viscous parameter $\alpha=0.1$ at
$r=a$; [121]).

The bottom panel of Fig. 9 shows a periodogram of the components of the
torque. In comparison to three-body systems, hydro-dynamical torques are
significantly more complex and spread over a wide range of harmonics.
Therefore, we expect the $\delta$SNR of the resulting phase perturbations to
be significantly larger than what individual torque component peaks may
suggest (following Eq. 51). Nevertheless, we can identify several peaks at
approximately $\sim 0.38f_{\rm K}$, $\sim 1.6f_{\rm{K}}$ and $\sim 3.9f_{\rm
K}$. Therefore, a high-SNR observation of GWs from such a system would in
principle allow one to reconstruct features in the flow of gas around in-
spiralling black holes. We note that this should be the case for any type of
gas-embedded binary, including stellar-mass binaries in accretion discs [e.g.
5, 79, 37, 107, 144, 43].

The results displayed in the bottom panel of Fig. 9 scale with a dimensionless
parameter $q_{\rm d}$:

$\displaystyle q_{\rm d}=\frac{1}{3\pi\nu_{\rm d}}\frac{\dot{M}}{M}(2\pi
f_{\rm K})^{-1},$ (69)

where $\dot{M}$ is the accretion rate onto the binary. The parameter $q_{\rm
d}$ can be roughly understood as the ratio of the disc’s scale mass and the
binary total mass, which depends on the binary separation [compare with e.g.
131]. For realistic accretion discs this parameter will be small when sources
enter the m-Hz band, of order $\sim 10^{-5}$ to few $\sim 10^{-1}$, depending
on the total binary mass and the properties of the disc.

We also note that while the above analysis and the majority of numerical work
on this problem has focused on discs of characteristic scale-height $h/r=0.1$,
steady-state modeling suggests that discs around such massive central objects
ought to be much thinner with $h/r\sim\mathcal{O}(10^{-2}-10^{-3})$ [59, 123,
130]. A series of recent studies have documented how the binary evolution and
accretion morphology are sensitive to the disc scale height [132, 44, 95, 58].
Of primary note, both the net gravitational torque on the binary and the
magnitude of torque fluctuations have been observed to grow meaningfully with
decreasing $h/r$ [132, 39]. Therefore, the results in Fig. 9 might be regarded
as conservative if thin, radiatively efficient discs around massive binaries
resemble those around their single-BH AGN-counterparts777But see also [31] who
find that binary heating of the inner disc can sustain thicker discs akin to
the $h/r=0.1$ case presented here..

###  V.4 EMRIs and IMRIs in accretion discs

Extreme and intermediate mass-ratio inspirals (E/IMRI) consist in a stellar-
mass BH orbiting a central SMBH. Such systems may form in the accretion discs
of active galactic nuclei at potentially high rates, via the process of BH-
accretion and disc interaction [127, 54, 93, 36, 107] and/or in-situ formation
of embedded stars [78, 40]. For this particular type of GW source,
characterising the environmental properties is critical, as a
misinterpretation of EEs can interfere with precision measurements of BH spin
and of the central metric’s multipolar structure [12, 124].

Simulations exploring E/IMRI coalescence in accretion discs find that the gas
exerts torques with fluctuating magnitude on the embedded BH, which roughly
averages out to match the expectations of linear perturbation theory and
planetary migration over many orbits [129, 64, 38]. In Fig. 10, we show a
time-series of the instantaneous torque experienced by an IMRI embedded in a
steady-state, iso-thermal and viscous accretion disc. The data are taken from
un-smoothed results of the simulations detailed in [38], i.e. 2D hydro-
dynamical Newtonian simulation performed with the code DISCO [51], in which an
IMRI slowly in-spirals on a prescribed, quasi-circular orbit. The parameters
for the simulation are a binary mass-ratio $q=10^{-3}$, a disc alpha viscosity
parameter $\alpha=0.03$, and a uniform disc scale height $h/r\approx 0.03$.
Here we additionally choose a conservative surface density scale for the disc
of $0.1$ g/cm2 [compare with e.g., 38], which already suffices to produce very
strong effects. Note that both the secular torques and the fluctuations scale
linearly with the disc’s scale density.

During the in-spiral, the IMRI experiences torque fluctuations on various
timescales comparable to the orbital timescale, which typically exceed their
orbit-averaged values by a factor of several hundred. Such strong fluctuations
are interpreted to be a consequence of gravitational forces arising from the
asymmetry in the mini-disc around the secondary BH. In the bottom panel of
Fig. 10, we showcase a periodogram of the acceleration suffered by the
embedded IMRI, normalised by the gravitational acceleration due to the central
object. We see that the force spectrum is rich with fluctuations at multiple
harmonics, though presents some distinct peaks at $\sim 0.13f_{\rm K}$, $\sim
1.1f_{\rm K}$ and $\sim 2.1f_{\rm K}$. Detecting such components in an IMRI
signal would provide information regarding the gas dynamics around the
embedded BH, as well as the global response of the disc to a small perturber.
This is particularly useful, as GW distortions are probes of the nearby gas
mass distribution, while all of our current information from EM observations
only probes the surface of such optically thick discs. Similarly to the
circumbinary disc case, the total torque spectrum on an embedded IMRI will
produce phase perturbations with a significantly larger $\delta$SNR than the
highest peaks in the periodogram would suggest (according to Eq. 51). We note
that at later stages of the in-spiral, the magnitude of fluctuations grows
(albeit remains weak compared to the dominance of GW induced evolution), while
the periodicity remains correlated with the binary orbital period. Note also
that the torque fluctuations for this system are enhanced by the imposed GW
in-spiral, an effect which is not explored in our demonstrated example of MBHB
interaction (Fig. 9).

Aside from the nonlinear torque evolution discovered in laminar accretion
flows, we comment here that accretion discs are also expected to be magnetized
and turbulent [10]. This may cause stochastic migration for lower mass BHs
(EMRIs), which would result in additional periodic torque fluctuations. For an
example, we refer the reader to the simulation results of [148], which explore
migration in turbulent accretion disc models. In short, it is expected that
gas-embedded IMRIs and EMRIs will be subject to a wide range of physical
processes that produce torque fluctuations in addition to the secular energy
and angular momentum fluxes (see also [156] for a thorough discussion).

Figure 10: Top panel: Time series of the simulated hydro-dynamical torque on
an intermediate mass-ratio inspiral ($q=10^{-3}$) in a viscous accretion disc
(see Sec. V.4). The magnitude is normalised by its own orbit-averaged value,
showcasing how fluctuations dominate the short time-scale dynamics of the
system. Bottom panel: The spectrum of force Fourier-components (in the form of
Eq. II.1) that arise from the simulated torque. We highlight three dominant
force peaks at different frequencies, related to features of the gas flow
around the perturber. The dashed lines denote the approximate SNR required to
detect individual force peaks, according to the monochromatic criterion (see
Sec. IV and Eq. 48), serving as a demonstration that such peaks may be
detectable for high SNR IMRIs and EMRIs. Note that for this realisation the
total $\delta$SNR of the total force spectrum is be significantly larger than
what is suggested by the individual peaks.

##  VI Conclusions

In this work, we have characterised aspects of the motion of perturbed
circular binaries (Sec. II) and the resulting GW emission, focusing
specifically on time-varying forces that do not cause any secular evolution in
the orbital elements. We constructed a waveform model, which we later use to
produce detectability criteria for the resulting periodic perturbations in the
GW phase (Sec. III and IV). Finally, we discussed a range of astrophysical
scenarios in which such time-varying forces occur naturally, and often at
significantly larger amplitudes that secular effects (Sec. V).

A major unanswered question remains: Are sufficiently large periodic phase
perturbations realised in nature, such that they may indeed be extracted from
realistic GW signals with near-future instruments? In Sec. V, we have shown
how many astrophysical scenarios will naturally produce strong force
fluctuations that greatly exceed the magnitude of orbit-averaged forces,
basing our arguments on the results of N-body and hydro-dynamical simulations.
Our examples show that the potential exists to detect such perturbations in
astrophysical sources that are expected for near future detectors. We
nevertheless stress that answering the aforementioned question requires a
thorough estimate of detection rates. Beyond the intrinsic modelling
uncertainties for the perturbations in question, the latter would have to
entail the following elements:

  * •

Population synthesis models for each of the GW source types and GW detectors
mentioned in this work.

  * •

A complete data analysis pipeline that takes all GW parameters into account,
including spins, various orientation angles and potentially the ulterior GW
perturbations discussed in Sec. II.

  * •

Accounting for the full force spectra rather than individual harmonics. This
is in particular relevant for the complex forces expected to act on gas-
embedded sources.

  * •

A thorough analysis of the odds of periodic perturbations being related to
astrophysical EEs rather than some unknown sources of noise.

While such an analysis is beyond the scope of this paper, we do stress that
the detection of even a single binary source that is clearly perturbed by
periodic EEs would be momentous, as it would represent unequivocal and
informative evidence of the interaction of a binary in the GW dominated phase
with its surroundings. Additionally, the results of this work highlight the
potential of considering periodic perturbations in GWs, in comparison with
e.g. secular ones. To reiterate the argument: In Sec. IV.4 we have discussed
how a hypothetically constant secular de-phasing, i.e. a phase offset, suffers
from being entirely degenerate with the source’s initial phase (for values
below $2\pi$). In other words, it is only the rate of variation of de-phasing
that is informative, rather than the magnitude of de-phasing itself.
Crucially, the periodic phase perturbations of the type analysed in this work
essentially maximise the rate of variation of the perturbation in the GW phase
by construction. They are therefore maximally detectable as an EE, given a
typical amplitude. Furthermore, they are less prone to result in ulterior
degeneracies, or in the degradation of the recovery of vacuum parameters (see
Fig. 7 and Sec. IV.4). Indeed, the question of when the accumulation of
secular de-phasing over long timescales out-competes the advantages of
periodic de-phasing also requires more investigation.

In the meantime, we suggest that the latter be included in future studies of
EEs by means of Eqs. 27 and 28 (or any other extension), in addition to
secular de-phasing prescriptions.

###### Acknowledgements.

L.Z. acknowledges support from ERC Starting Grant No. 121817–BlackHoleMergs.
A.A.T. acknowledges support from the European Union’s Horizon 2020 and Horizon
Europe research and innovation programs under the Marie Skłodowska-Curie grant
agreements No. 847523 and 101103134. D.J.D. and C.T. received support from the
Danish Independent Research Fund through Sapere Aude Starting grant No.
121587. Z.H. acknowledges support from NSF grant AST-2006176 and NASA grants
80NSSC22K0822 and 80NSSC24K0440. AD acknowledges support from NSF grant
AST-2319441.

## References

  * Abadie et al. [2010] Abadie, J., Abbott, B. P., Abbott, R., et al. 2010, Classical and Quantum Gravity, 27, 173001, doi: 10.1088/0264-9381/27/17/173001
  * Abbott et al. [2020] Abbott, R., Abbott, T. D., Abraham, S., et al. 2020, Phys. Rev. Lett. , 125, 101102, doi: 10.1103/PhysRevLett.125.101102
  * Amaro-Seoane et al. [2017] Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786.  https://arxiv.org/abs/1702.00786
  * Amaro-Seoane et al. [2022] Amaro-Seoane, P., Andrews, J., Arca Sedda, M., et al. 2022, arXiv e-prints, arXiv:2203.06016.  https://arxiv.org/abs/2203.06016
  * Antoni et al. [2019] Antoni, A., MacLeod, M., & Ramirez-Ruiz, E. 2019, Astrophys. J. , 884, 22, doi: 10.3847/1538-4357/ab3466
  * Antonini & Rasio [2016] Antonini, F., & Rasio, F. A. 2016, Astrophys. J. , 831, 187, doi: 10.3847/0004-637X/831/2/187
  * Antonini et al. [2017] Antonini, F., Toonen, S., & Hamers, A. S. 2017, Astrophys. J. , 841, 77, doi: 10.3847/1538-4357/aa6f5e
  * Askar et al. [2017] Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36, doi: 10.1093/mnrasl/slw177
  * Avara et al. [2023] Avara, M. J., Krolik, J. H., Campanelli, M., et al. 2023, arXiv e-prints, arXiv:2305.18538, doi: 10.48550/arXiv.2305.18538
  * Balbus & Hawley [1991] Balbus, S. A., & Hawley, J. F. 1991, Astrophys. J. , 376, 214, doi: 10.1086/170270
  * Barack & Cutler [2007] Barack, L., & Cutler, C. 2007, Phys. Rev. D, 75, 042003, doi: 10.1103/PhysRevD.75.042003
  * Barausse et al. [2014] Barausse, E., Cardoso, V., & Pani, P. 2014, Phys. Rev. D, 89, 104059, doi: 10.1103/PhysRevD.89.104059
  * Barausse & Rezzolla [2008] Barausse, E., & Rezzolla, L. 2008, Phys. Rev. D, 77, 104027, doi: 10.1103/PhysRevD.77.104027
  * Barnes & Hernquist [1996] Barnes, J. E., & Hernquist, L. 1996, Astrophys. J. , 471, 115, doi: 10.1086/177957
  * Begelman et al. [1980] Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Nature (London), 287, 307, doi: 10.1038/287307a0
  * Belczynski et al. [2002] Belczynski, K., Kalogera, V., & Bulik, T. 2002, Astrophys. J. , 572, 407, doi: 10.1086/340304
  * Blaes et al. [2002] Blaes, O., Lee, M. H., & Socrates, A. 2002, Astrophys. J. , 578, 775, doi: 10.1086/342655
  * Blanchet [2014a] Blanchet, L. 2014a, Living Reviews in Relativity, 17, 2, doi: 10.12942/lrr-2014-2
  * Blanchet [2014b] —. 2014b, Living Rev. Relat., 17, 2, doi: 10.12942/lrr-2014-2
  * Bonetti et al. [2017] Bonetti, M., Barausse, E., Faye, G., Haardt, F., & Sesana, A. 2017, Classical and Quantum Gravity, 34, 215004, doi: 10.1088/1361-6382/aa8da5
  * Camilloni et al. [2023] Camilloni, F., Grignani, G., Harmark, T., Orselli, M., & Pica, D. 2023, arXiv e-prints, arXiv:2310.06894, doi: 10.48550/arXiv.2310.06894
  * Cardoso & Duque [2020] Cardoso, V., & Duque, F. 2020, Phys. Rev. D, 101, 064028, doi: 10.1103/PhysRevD.101.064028
  * Cardoso & Maselli [2020a] Cardoso, V., & Maselli, A. 2020a, A&A, 644, A147, doi: 10.1051/0004-6361/202037654
  * Cardoso & Maselli [2020b] —. 2020b, A&A, 644, A147, doi: 10.1051/0004-6361/202037654
  * Chakrabarti [1993] Chakrabarti, S. K. 1993, Astrophys. J. , 411, 610, doi: 10.1086/172863
  * Chamberlain et al. [2019a] Chamberlain, K., Moore, C. J., Gerosa, D., & Yunes, N. 2019a, Phys. Rev. D, 99, 024025, doi: 10.1103/PhysRevD.99.024025
  * Chamberlain et al. [2019b] —. 2019b, Phys. Rev. D, 99, 024025, doi: 10.1103/PhysRevD.99.024025
  * Chandramouli & Yunes [2022] Chandramouli, R. S., & Yunes, N. 2022, Phys. Rev. D, 105, 064009, doi: 10.1103/PhysRevD.105.064009
  * Chen [2021] Chen, X. 2021, in Handbook of Gravitational Wave Astronomy, 39, doi: 10.1007/978-981-15-4702-7_39-1
  * Chirikov [1979] Chirikov, B. V. 1979, Phys. Rep., 52, 263, doi: 10.1016/0370-1573(79)90023-1
  * Cocchiararo et al. [2024] Cocchiararo, F., Franchini, A., Lupi, A., & Sesana, A. 2024, arXiv e-prints, arXiv:2402.05175, doi: 10.48550/arXiv.2402.05175
  * Cole et al. [2022] Cole, P. S., Coogan, A., Kavanagh, B. J., & Bertone, G. 2022, arXiv e-prints, arXiv:2207.07576.  https://arxiv.org/abs/2207.07576
  * Colpi et al. [2024] Colpi, M., Danzmann, K., Hewitson, M., et al. 2024, arXiv e-prints, arXiv:2402.07571, doi: 10.48550/arXiv.2402.07571
  * Cutler & Flanagan [1994] Cutler, C., & Flanagan, É. E. 1994, Phys. Rev. D, 49, 2658, doi: 10.1103/PhysRevD.49.2658
  * Damour et al. [2008] Damour, T., Nagar, A., Hannam, M., Husa, S., & Brügmann, B. 2008, Phys. Rev. D, 78, 044039, doi: 10.1103/PhysRevD.78.044039
  * DeLaurentiis et al. [2022] DeLaurentiis, S., Epstein-Martin, M., & Haiman, Z. 2022, arXiv e-prints, arXiv:2212.02650, doi: 10.48550/arXiv.2212.02650
  * Dempsey et al. [2022] Dempsey, A. M., Li, H., Mishra, B., & Li, S. 2022, Astrophys. J. , 940, 155, doi: 10.3847/1538-4357/ac9d92
  * Derdzinski et al. [2021a] Derdzinski, A., D’Orazio, D., Duffell, P., Haiman, Z., & MacFadyen, A. 2021a, MNRAS, 501, 3540, doi: 10.1093/mnras/staa3976
  * Derdzinski et al. [2021b] —. 2021b, MNRAS, 501, 3540, doi: 10.1093/mnras/staa3976
  * Derdzinski & Mayer [2023] Derdzinski, A., & Mayer, L. 2023, In situ extreme mass ratio inspirals via subparsec formation and migration of stars in thin, gravitationally unstable AGN discs, doi: 10.1093/mnras/stad749
  * Derdzinski et al. [2019] Derdzinski, A. M., D’Orazio, D., Duffell, P., Haiman, Z., & MacFadyen, A. 2019, MNRAS, 486, 2754, doi: 10.1093/mnras/stz1026
  * Destounis et al. [2022] Destounis, K., Kulathingal, A., Kokkotas, K. D., & Papadopoulos, G. O. 2022, arXiv e-prints, arXiv:2210.09357, doi: 10.48550/arXiv.2210.09357
  * Dittmann et al. [2024] Dittmann, A. J., Dempsey, A. M., & Li, H. 2024, Astrophys. J. , 964, 61, doi: 10.3847/1538-4357/ad23ce
  * Dittmann & Ryan [2022] Dittmann, A. J., & Ryan, G. 2022, MNRAS, 513, 6158, doi: 10.1093/mnras/stac935
  * Dittmann et al. [2023] Dittmann, A. J., Ryan, G., & Miller, M. C. 2023, ApJ, 949, L30, doi: 10.3847/2041-8213/acd183
  * Dominik et al. [2013] Dominik, M., Belczynski, K., Fryer, C., et al. 2013, Astrophys. J. , 779, 72, doi: 10.1088/0004-637X/779/1/72
  * D’Orazio & Charisi [2023] D’Orazio, D. J., & Charisi, M. 2023, arXiv e-prints, arXiv:2310.16896, doi: 10.48550/arXiv.2310.16896
  * D’Orazio et al. [2024] D’Orazio, D. J., Duffell, P. C., & Tiede, C. 2024, arXiv e-prints, arXiv:2403.05629, doi: 10.48550/arXiv.2403.05629
  * D’Orazio & Loeb [2020a] D’Orazio, D. J., & Loeb, A. 2020a, Phys. Rev. D, 101, 083031, doi: 10.1103/PhysRevD.101.083031
  * D’Orazio & Loeb [2020b] —. 2020b, Phys. Rev. D, 101, 083031, doi: 10.1103/PhysRevD.101.083031
  * Duffell [2016] Duffell, P. C. 2016, ApJS, 226, 2, doi: 10.3847/0067-0049/226/1/2
  * Duffell et al. [2024] Duffell, P. C., Dittmann, A. J., D’Orazio, D. J., et al. 2024, arXiv e-prints, arXiv:2402.13039, doi: 10.48550/arXiv.2402.13039
  * Ewing et al. [2021] Ewing, B., Sachdev, S., Borhanian, S., & Sathyaprakash, B. S. 2021, Phys. Rev. D, 103, 023025, doi: 10.1103/PhysRevD.103.023025
  * Fabj et al. [2020] Fabj, G., Nasim, S. S., Caban, F., et al. 2020, MNRAS, 499, 2608, doi: 10.1093/mnras/staa3004
  * Farris et al. [2015] Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z. 2015, MNRAS, 447, L80, doi: 10.1093/mnrasl/slu184
  * Foreman-Mackey et al. [2013] Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Franchini et al. [2024] Franchini, A., Bonetti, M., Lupi, A., & Sesana, A. 2024, arXiv e-prints, arXiv:2401.10331, doi: 10.48550/arXiv.2401.10331
  * Franchini et al. [2022] Franchini, A., Lupi, A., & Sesana, A. 2022, ApJ, 929, L13, doi: 10.3847/2041-8213/ac63a2
  * Frank et al. [2002] Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition
  * Gair et al. [2013] Gair, J. R., Vallisneri, M., Larson, S. L., & Baker, J. G. 2013, Living Rev. Relat., 16, 7, doi: 10.12942/lrr-2013-7
  * Garg et al. [2024] Garg, M., Derdzinski, A., Tiwari, S., Gair, J., & Mayer, L. 2024, arXiv e-prints, arXiv:2402.14058, doi: 10.48550/arXiv.2402.14058
  * Garg et al. [2022] Garg, M., Derdzinski, A., Zwick, L., Capelo, P. R., & Mayer, L. 2022, MNRAS, 517, 1339, doi: 10.1093/mnras/stac2711
  * Gerosa et al. [2019] Gerosa, D., Ma, S., Wong, K. W. K., et al. 2019, Phys. Rev. D, 99, 103004, doi: 10.1103/PhysRevD.99.103004
  * Goldreich & Tremaine [1980] Goldreich, P., & Tremaine, S. 1980, Astrophys. J. , 241, 425, doi: 10.1086/158356
  * Haiman et al. [2009] Haiman, Z., Kocsis, B., Menou, K., Lippai, Z., & Frei, Z. 2009, Classical and Quantum Gravity, 26, 094032, doi: 10.1088/0264-9381/26/9/094032
  * Hayashi et al. [2023] Hayashi, T., Suto, Y., & Trani, A. A. 2023, Astrophys. J. , 958, 26, doi: 10.3847/1538-4357/acf4f6
  * Hoang et al. [2019] Hoang, B.-M., Naoz, S., Kocsis, B., Farr, W. M., & McIver, J. 2019, ApJ, 875, L31, doi: 10.3847/2041-8213/ab14f7
  * Inayoshi et al. [2017a] Inayoshi, K., Hirai, R., Kinugawa, T., & Hotokezaka, K. 2017a, MNRAS, 468, 5020, doi: 10.1093/mnras/stx757
  * Inayoshi et al. [2017b] Inayoshi, K., Tamanini, N., Caprini, C., & Haiman, Z. 2017b, Phys. Rev. D, 96, 063014, doi: 10.1103/PhysRevD.96.063014
  * Katz et al. [2021] Katz, M. L., Chua, A. J. K., Speri, L., Warburton, N., & Hughes, S. A. 2021, Phys. Rev. D, 104, 064047, doi: 10.1103/PhysRevD.104.064047
  * Kavanagh et al. [2020] Kavanagh, B. J., Nichols, D. A., Bertone, G., & Gaggero, D. 2020, Phys. Rev. D, 102, 083006, doi: 10.1103/PhysRevD.102.083006
  * Kejriwal et al. [2023] Kejriwal, S., Speri, L., & Chua, A. J. K. 2023, arXiv e-prints, arXiv:2312.13028, doi: 10.48550/arXiv.2312.13028
  * Kimball et al. [2021] Kimball, C., Talbot, C., Berry, C. P. L., et al. 2021, ApJ, 915, L35, doi: 10.3847/2041-8213/ac0aef
  * Kocsis et al. [2011] Kocsis, B., Yunes, N., & Loeb, A. 2011, Phys. Rev. D, 84, 024032, doi: 10.1103/PhysRevD.84.024032
  * Krauth et al. [2023] Krauth, L. M., Davelaar, J., Haiman, Z., et al. 2023, MNRAS, 526, 5441, doi: 10.1093/mnras/stad3095
  * Laeuger et al. [2023] Laeuger, A., Seymour, B., Chen, Y., & Yu, H. 2023, arXiv e-prints, arXiv:2310.16799, doi: 10.48550/arXiv.2310.16799
  * Lai & Muñoz [2022] Lai, D., & Muñoz, D. J. 2022, arXiv e-prints, arXiv:2211.00028.  https://arxiv.org/abs/2211.00028
  * Levin [2007] Levin, Y. 2007, MNRAS, 374, 515, doi: 10.1111/j.1365-2966.2006.11155.x
  * Li & Lai [2022] Li, R., & Lai, D. 2022, MNRAS, 517, 1602, doi: 10.1093/mnras/stac2577
  * Liu et al. [2022] Liu, B., D’Orazio, D. J., Vigna-Gómez, A., & Samsing, J. 2022, Phys. Rev. D, 106, 123010, doi: 10.1103/PhysRevD.106.123010
  * Liu & Lai [2018] Liu, B., & Lai, D. 2018, Astrophys. J. , 863, 68, doi: 10.3847/1538-4357/aad09f
  * Liu et al. [2015] Liu, B., Muñoz, D. J., & Lai, D. 2015, MNRAS, 447, 747, doi: 10.1093/mnras/stu2396
  * Maggiore [2018] Maggiore, M. 2018, Gravitational Waves: Volume 2: Astrophysics and Cosmology, doi: 10.1093/oso/9780198570899.001.0001
  * Maggiore et al. [2020] Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, J. Cosmology Astropart. Phys., 2020, 050, doi: 10.1088/1475-7516/2020/03/050
  * Mapelli [2021] Mapelli, M. 2021, in Handbook of Gravitational Wave Astronomy, 16, doi: 10.1007/978-981-15-4702-7_16-1
  * Mardling [2008] Mardling, R. A. 2008, in The Cambridge N-Body Lectures, ed. S. J. Aarseth, C. A. Tout, & R. A. Mardling, Vol. 760, 59, doi: 10.1007/978-1-4020-8431-7_3
  * Mei et al. [2021] Mei, J., Bai, Y.-Z., Bao, J., et al. 2021, Progr. Theor. Exp. Phys., 2021, 05A107, doi: 10.1093/ptep/ptaa114
  * Meiron et al. [2017] Meiron, Y., Kocsis, B., & Loeb, A. 2017, Astrophys. J. , 834, 200, doi: 10.3847/1538-4357/834/2/200
  * Milosavljević & Merritt [2003] Milosavljević, M., & Merritt, D. 2003, Astrophys. J. , 596, 860, doi: 10.1086/378086
  * Moore et al. [2017] Moore, C. J., Chua, A. J. K., & Gair, J. R. 2017, Classical and Quantum Gravity, 34, 195009, doi: 10.1088/1361-6382/aa85fa
  * O’Leary et al. [2007] O’Leary, R. M., O’Shaughnessy, R., & Rasio, F. A. 2007, Phys. Rev. D, 76, 061504, doi: 10.1103/PhysRevD.76.061504
  * Owen et al. [2023] Owen, C. B., Haster, C.-J., Perkins, S., Cornish, N. J., & Yunes, N. 2023, arXiv e-prints, arXiv:2301.11941.  https://arxiv.org/abs/2301.11941
  * Pan et al. [2021] Pan, Z., Lyu, Z., & Yang, H. 2021, Phys. Rev. D, 104, 063007, doi: 10.1103/PhysRevD.104.063007
  * Park et al. [2017] Park, D., Kim, C., Lee, H. M., Bae, Y.-B., & Belczynski, K. 2017, MNRAS, 469, 4665, doi: 10.1093/mnras/stx1015
  * Penzlin et al. [2022] Penzlin, A. B. T., Kley, W., Audiffren, H., & Schäfer, C. M. 2022, A&A, 660, A101, doi: 10.1051/0004-6361/202141399
  * Poisson & Will [2014] Poisson, E., & Will, C. M. 2014, Gravity
  * Polcar et al. [2022] Polcar, L., Lukes-Gerakopoulos, G., & Witzany, V. 2022, Phys. Rev. D, 106, 044069, doi: 10.1103/PhysRevD.106.044069
  * Portegies Zwart & McMillan [2000] Portegies Zwart, S. F., & McMillan, S. L. W. 2000, Astrophys. J. , 528, L17
  * Randall & Xianyu [2019a] Randall, L., & Xianyu, Z.-Z. 2019a, arXiv e-prints, arXiv:1902.08604, doi: 10.48550/arXiv.1902.08604
  * Randall & Xianyu [2019b] —. 2019b, Astrophys. J. , 878, 75, doi: 10.3847/1538-4357/ab20c6
  * Robson et al. [2019] Robson, T., Cornish, N. J., & Liu, C. 2019, Class. Quantum Gravity, 36, 105011, doi: 10.1088/1361-6382/ab1101
  * Robson et al. [2018] Robson, T., Cornish, N. J., Tamanini, N., & Toonen, S. 2018, Phys. Rev. D, 98, 064012, doi: 10.1103/PhysRevD.98.064012
  * Rodriguez & Antonini [2018] Rodriguez, C. L., & Antonini, F. 2018, Astrophys. J. , 863, 7, doi: 10.3847/1538-4357/aacea4
  * Rodriguez et al. [2016a] Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016a, Phys. Rev. D, 93, 084029, doi: 10.1103/PhysRevD.93.084029
  * Rodriguez et al. [2016b] Rodriguez, C. L., Haster, C.-J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016b, ApJ, 824, L8, doi: 10.3847/2041-8205/824/1/L8
  * Roedig et al. [2012] Roedig, C., Sesana, A., Dotti, M., et al. 2012, A&A, 545, A127, doi: 10.1051/0004-6361/201219986
  * Rowan et al. [2022] Rowan, C., Boekholt, T., Kocsis, B., & Haiman, Z. 2022, arXiv e-prints, arXiv:2212.06133, doi: 10.48550/arXiv.2212.06133
  * Ryan [1995] Ryan, F. D. 1995, Phys. Rev. D, 52, 5707, doi: 10.1103/PhysRevD.52.5707
  * Ryu et al. [2018] Ryu, T., Perna, R., Haiman, Z., Ostriker, J. P., & Stone, N. C. 2018, MNRAS, 473, 3410, doi: 10.1093/mnras/stx2524
  * Sadowski et al. [2008] Sadowski, A., Belczynski, K., Bulik, T., et al. 2008, Astrophys. J. , 676, 1162, doi: 10.1086/528932
  * Samsing [2018] Samsing, J. 2018, Phys. Rev. D, 97, 103014, doi: 10.1103/PhysRevD.97.103014
  * Samsing & D’Orazio [2018] Samsing, J., & D’Orazio, D. J. 2018, MNRAS, doi: 10.1093/mnras/sty2334
  * Samsing et al. [2020] Samsing, J., D’Orazio, D. J., Kremer, K., Rodriguez, C. L., & Askar, A. 2020, Phys. Rev. D, 101, 123010, doi: 10.1103/PhysRevD.101.123010
  * Samsing et al. [2024a] Samsing, J., Hendriks, K., Zwick, L., D’Orazio, D. J., & Liu, B. 2024a, arXiv e-prints, arXiv:2403.05625, doi: 10.48550/arXiv.2403.05625
  * Samsing et al. [2024b] —. 2024b, arXiv e-prints, arXiv:2403.05625, doi: 10.48550/arXiv.2403.05625
  * Samsing et al. [2014] Samsing, J., MacLeod, M., & Ramirez-Ruiz, E. 2014, Astrophys. J. , 784, 71, doi: 10.1088/0004-637X/784/1/71
  * Santini et al. [2023] Santini, A., Gerosa, D., Cotesta, R., & Berti, E. 2023, Phys. Rev. D, 108, 083033, doi: 10.1103/PhysRevD.108.083033
  * Sberna et al. [2022] Sberna, L., Babak, S., Marsat, S., et al. 2022, Phys. Rev. D, 106, 064056, doi: 10.1103/PhysRevD.106.064056
  * Sesana [2013] Sesana, A. 2013, Classical and Quantum Gravity, 30, 224014, doi: 10.1088/0264-9381/30/22/224014
  * Sesana [2016] —. 2016, Phys. Rev. Lett. , 116, 231102, doi: 10.1103/PhysRevLett.116.231102
  * Shakura & Sunyaev [1973] Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Silsbee & Tremaine [2017] Silsbee, K., & Tremaine, S. 2017, Astrophys. J. , 836, 39, doi: 10.3847/1538-4357/aa5729
  * Sirko & Goodman [2003] Sirko, E., & Goodman, J. 2003, MNRAS, 341, 501, doi: 10.1046/j.1365-8711.2003.06431.x
  * Speri et al. [2022] Speri, L., Antonelli, A., Sberna, L., et al. 2022, arXiv e-prints, arXiv:2207.10086.  https://arxiv.org/abs/2207.10086
  * Springel et al. [2005] Springel, V., Di Matteo, T., & Hernquist, L. 2005, MNRAS, 361, 776, doi: 10.1111/j.1365-2966.2005.09238.x
  * Strokov et al. [2022] Strokov, V., Fragione, G., Wong, K. W. K., Helfer, T., & Berti, E. 2022, Phys. Rev. D, 105, 124048, doi: 10.1103/PhysRevD.105.124048
  * Tagawa et al. [2020] Tagawa, H., Haiman, Z., & Kocsis, B. 2020, Astrophys. J. , 898, 25, doi: 10.3847/1538-4357/ab9b8c
  * Tamanini et al. [2020] Tamanini, N., Klein, A., Bonvin, C., Barausse, E., & Caprini, C. 2020, Phys. Rev. D, 101, 063002, doi: 10.1103/PhysRevD.101.063002
  * Tanaka et al. [2002] Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, Astrophys. J. , 565, 1257, doi: 10.1086/324713
  * Thompson et al. [2005] Thompson, T. A., Quataert, E., & Murray, N. 2005, Astrophys. J. , 630, 167, doi: 10.1086/431923
  * Tiede et al. [2024] Tiede, C., D’Orazio, D. J., Zwick, L., & Duffell, P. C. 2024, Astrophys. J. , 964, 46, doi: 10.3847/1538-4357/ad2613
  * Tiede et al. [2020] Tiede, C., Zrake, J., MacFadyen, A., & Haiman, Z. 2020, Astrophys. J. , 900, 43, doi: 10.3847/1538-4357/aba432
  * Tiwari et al. [2024] Tiwari, A., Vijaykumar, A., Kapadia, S. J., Fragione, G., & Chatterjee, S. 2024, MNRAS, 527, 8586, doi: 10.1093/mnras/stad3749
  * Toonen et al. [2018] Toonen, S., Perets, H. B., & Hamers, A. S. 2018, A&A, 610, A22, doi: 10.1051/0004-6361/201731874
  * Torres-Orjuela et al. [2019] Torres-Orjuela, A., Chen, X., Cao, Z., Amaro-Seoane, P., & Peng, P. 2019, Phys. Rev. D, 100, 063012, doi: 10.1103/PhysRevD.100.063012
  * Toubiana et al. [2021a] Toubiana, A., Sberna, L., Caputo, A., et al. 2021a, Phys. Rev. Lett. , 126, 101105, doi: 10.1103/PhysRevLett.126.101105
  * Toubiana et al. [2021b] —. 2021b, Phys. Rev. Lett. , 126, 101105, doi: 10.1103/PhysRevLett.126.101105
  * Trani et al. [2024] Trani, A. A., Quaini, S., & Colpi, M. 2024, A&A, 683, A135, doi: 10.1051/0004-6361/202347920
  * Trani et al. [2022] Trani, A. A., Rastello, S., Di Carlo, U. N., et al. 2022, MNRAS, 511, 1362, doi: 10.1093/mnras/stac122
  * Trani & Spera [2023] Trani, A. A., & Spera, M. 2023, IAU Symposium, 362, 404, doi: 10.1017/S1743921322001818
  * Vigna-Gómez et al. [2021] Vigna-Gómez, A., Toonen, S., Ramirez-Ruiz, E., et al. 2021, ApJ, 907, L19, doi: 10.3847/2041-8213/abd5b7
  * Vitale et al. [2017] Vitale, S., Lynch, R., Sturani, R., & Graff, P. 2017, Classical and Quantum Gravity, 34, 03LT01, doi: 10.1088/1361-6382/aa552e
  * Westernacher-Schneider et al. [2022] Westernacher-Schneider, J. R., Zrake, J., MacFadyen, A., & Haiman, Z. 2022, Phys. Rev. D, 106, 103010, doi: 10.1103/PhysRevD.106.103010
  * Whitehead et al. [2023] Whitehead, H., Rowan, C., Boekholt, T., & Kocsis, B. 2023, arXiv e-prints, arXiv:2309.11561, doi: 10.48550/arXiv.2309.11561
  * Will [2014a] Will, C. M. 2014a, Classical and Quantum Gravity, 31, 244001, doi: 10.1088/0264-9381/31/24/244001
  * Will [2014b] —. 2014b, Phys. Rev. D, 89, 044043, doi: 10.1103/PhysRevD.89.044043
  * Wong et al. [2019] Wong, K. W. K., Baibhav, V., & Berti, E. 2019, MNRAS, 488, 5665, doi: 10.1093/mnras/stz2077
  * Wu et al. [2024] Wu, Y., Chen, Y.-X., & Lin, D. N. C. 2024, MNRAS, 528, L127, doi: 10.1093/mnrasl/slad183
  * Xuan et al. [2022] Xuan, Z., Naoz, S., & Chen, X. 2022, arXiv e-prints, arXiv:2210.03129, doi: 10.48550/arXiv.2210.03129
  * Xuan et al. [2023] —. 2023, Phys. Rev. D, 107, 043009, doi: 10.1103/PhysRevD.107.043009
  * Yamada et al. [2019] Yamada, K., Narikawa, T., & Tanaka, T. 2019, Progress of Theoretical and Experimental Physics, 2019, 103E01, doi: 10.1093/ptep/ptz103
  * Yan et al. [2023] Yan, H., Chen, X., & Torres-Orjuela, A. 2023, Phys. Rev. D, 107, 103044, doi: 10.1103/PhysRevD.107.103044
  * Yunes et al. [2011] Yunes, N., Miller, M. C., & Thornburg, J. 2011, Phys. Rev. D, 83, 044030, doi: 10.1103/PhysRevD.83.044030
  * Zevin et al. [2021] Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, Astrophys. J. , 910, 152, doi: 10.3847/1538-4357/abe40e
  * Zwick et al. [2023] Zwick, L., Capelo, P. R., & Mayer, L. 2023, MNRAS, 521, 4645, doi: 10.1093/mnras/stad707
  * Zwick et al. [2022] Zwick, L., Derdzinski, A., Garg, M., Capelo, P. R., & Mayer, L. 2022, MNRAS, 511, 6143, doi: 10.1093/mnras/stac299
  * Zwick & Mayer [2023] Zwick, L., & Mayer, L. 2023, MNRAS, 526, 2754, doi: 10.1093/mnras/stad2942