# Accuracy Requirements: Assessing the Importance of First Post-Adiabatic
Terms for Small-Mass-Ratio Binaries

Ollie Burke<EMAIL_ADDRESS>Laboratoire des 2 Infinis – Toulouse
(L2IT-IN2P3), Universitè de Toulouse, CNRS, UPS, F-31062 Toulouse Cedex 9,
France  Gabriel Andres Piovano<EMAIL_ADDRESS>School of Mathematics
and Statistics, University College Dublin, Belfield, Dublin 4, Ireland  Niels
Warburton  School of Mathematics and Statistics, University College Dublin,
Belfield, Dublin 4, Ireland  Philip Lynch  Max Planck Institute for
Gravitational Physics (Albert Einstein Institute) Am Mühlenberg 1, 14476
Potsdam, Germany  Lorenzo Speri  Max Planck Institute for Gravitational
Physics (Albert Einstein Institute) Am Mühlenberg 1, 14476 Potsdam, Germany
Chris Kavanagh  School of Mathematics and Statistics, University College
Dublin, Belfield, Dublin 4, Ireland  Barry Wardell  School of Mathematics and
Statistics, University College Dublin, Belfield, Dublin 4, Ireland  Adam Pound
School of Mathematical Sciences and STAG Research Centre, University of
Southampton, Southampton, SO17 1BJ, United Kingdom  Leanne Durkan  Center of
Gravitational Physics, University of Texas at Austin, Austin, Texas, USA,
78712  Jeremy Miller  Shamoon College of Engineering, 56 Bialik Street, Be’er
Sheva 8410802, Israel.

###### Abstract

We investigate the impact of post-adiabatic (1PA) terms on parameter
estimation for extreme and intermediate mass-ratio inspirals using state-of-
the-art waveform models. Our analysis is the first to employ Bayesian
inference to assess systematic errors for 1PA waveforms. We find that
neglecting 1PA terms introduces significant biases for the (small) mass ratio
$\epsilon\gtrsim 10^{-5}$ for quasi circular orbits in Schwarzschild
spacetime, which can be mitigated with resummed 3PN expressions at 1PA order.
Moreover, we show that the secondary spin is strongly correlated with the
other intrinsic parameters, and it can not be constrained for
$\epsilon\lesssim 10^{-5}$. Finally, we highlight the need for addressing
eccentric waveform systematics in the small-mass-ratio regime, as they yield
stronger biases than the circular limit in both intrinsic and extrinsic
parameters.

Gravitational Waves, Extreme/Intermediate Mass Ratio Inspirals, Gravitational
Self-Force, Accuracy Requirements, Waveform Systematics, Bayesian Inference

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

One of the most promising sources for the Laser Interferometer Space Antenna
(LISA), the first planned space-based gravitational-wave detector, are
extreme-mass-ratio inspirals (EMRIs). An EMRI is the slow inspiral of a
stellar-mass compact object (CO) of mass $\mu\sim 10^{0-2}M_{\odot}$ into a
massive black hole (MBH) with mass $M\sim 10^{5-7}M_{\odot}$. The successful
detection of an EMRI within the LISA data stream will be a groundbreaking
achievement, offering outstanding tests of general relativity and unique
insights into the fundamental nature of the central MBH barack2007using ;
gair2013testing ; Maselli:2020zgv ; Maselli:2021men ; Barsanti:2022ana .
Extracting the maximum information from LISA will only be possible with
sufficiently high-fidelity waveform models. The source modeling of the
gravitational-wave (GW) signal is exceptionally complicated, requiring
sophisticated mathematical techniques arising from black hole perturbation
theory and gravitational self-force (GSF) theory barack2018self ;
Pound:2021qin . Work in these areas dates back to the late 1950s, but, to this
day, the mathematical details for generic orbits (eccentric and precessing
inspirals) into rotating BHs, when accounting for GW emission, are still being
developed. Recent progress in this waveform-modeling program includes rapid
generation of waveforms Speri:2023jte ; katz2021fast ; chua2021rapid ;
VanDeMeent:2018cgn ; McCart:2021upc ; Lynch:2021ogr ; Lynch:2023gpu , the
inclusion of transient resonances Isoyama:2021jjd ; Gupta:2022fbe ;
speri2021assessing ; Lynch2022  and of the small companion’s spin
Akcay:2019bvk ; Piovano:2021iwv ; Skoupy:2021asz ; Mathews:2021rod ;
Drummond:2022xej ; Drummond:2022efc ; Skoupy:2023lih , generation of adiabatic
waveforms for generic orbits around a Kerr BH Hughes:2021exa ; Isoyama:2021jjd
and of some specific post-adiabatic corrections on generic orbits
vandeMeent:2017bcc , and generation of waveforms including all effects at
first post-adiabatic order in the case of quasicircular, nonspinning binaries
Wardell:2021fyy .

Not only are EMRIs challenging to model, but they prove to be one of the
hardest problems in LISA data analysis MockLISADataChallengeTaskForce:2009wir
; Chua:2019wgs ; chua2022nonlocal . Due to the rich structure of the waveform
and sheer volume of parameter space, grid-based searches such as those used by
LIGO will be impossible gair2004event . Furthermore, EMRI signals are
typically long-lived and may remain within the LISA sensitivity band for
multiple years. The number of cycles scales with the inverse of the (small)
mass ratio, $\epsilon=\mu/M$, implying that one could observe hundreds of
thousands of orbits. Thus, we expect to constrain parameters to sub-percent
level precision Babak:17aa ; burke2021extreme .

The detection and eventual characterization of an EMRI signal require waveform
models that are faithful to the GW signal buried within the noise of the
instrument. A detailed knowledge of the GSF is required, which involves
solving the Einstein field equations order by order in powers of the small
perturbative variable $\epsilon\ll 1$. Using a two-timescale expansion
Hinderer:2008dm ; Pound:2021qin , it has been shown that the CO’s orbital
phase has an expansion of the form111In this expansion, we have neglected the
effect of resonances that appear at fractional powers of the mass ratio
$\mathcal{O}(\epsilon^{-1/2})$. For the class of orbits considered in this
work, resonant orbits do not exist and thus will be neglected.

$\phi=\underbrace{\epsilon^{-1}\phi_{0}}_{\text{0PA}}+\underbrace{\phi_{1}}_{\text{1PA}}+\underbrace{\epsilon\phi_{2}}_{\text{2PA}}+O(\epsilon^{2}).$
(1)

The leading order in this expansion represents the _adiabatic_ (0PA)
evolution, which can be understood as a slow evolution of geodesic orbital
parameters due to the orbit-averaged dissipative piece of the first-order term
in the self-force (arising from the $\sim\epsilon$ term in the metric
perturbation). The forcing functions that drive this evolution have been
computed at 0PA order for generic bound orbits Fujita:2009us ; Flanagan:2012kg
; Isoyama:2021jjd ; Hughes:2021exa . However, it is known that we also require
the _first post–adiabatic_ (1PA) contribution, $\phi_{1}$, to ensure
sufficiently accurate waveform models Hinderer:2008dm ; LISA:2022kgy ;
LISA:2022yao . This 1PA term involves not only the conservative and
dissipative-oscillatory corrections at first order in the mass ratio, but also
the orbit-averaged dissipative effects of the second-order self-force (arising
from the $\sim\epsilon^{2}$ term in the metric perturbation). Finally, linear
effects from the secondary spin and evolving spin and mass of the primary feed
into this 1PA term Mathews:2021rod ; Miller:2020bft . Since the _second post-
adiabatic_ (2PA) corrections induce a contribution to the orbital phase which
scales linearly with the mass ratio, they are deemed sufficiently small to be
unnecessary for parameter estimation of EMRIs. Thus, knowledge of the second-
order self-force corrections (and linear corrections from the secondary spin)
is thought to be both necessary _and sufficient_ for EMRI data analysis. In
2021, a major breakthrough in accuracy was achieved: the second-order metric
perturbation was computed, allowing the construction of the first complete 1PA
waveforms for quasicircular orbits in Schwarzschild spacetimes Wardell:2021fyy
(building on Refs. Pound:2019lzj ; Miller:2020bft ; Warburton:2021kwk ). This
waveform showed remarkable agreement with numerical relativity (NR) waveforms,
even for $\epsilon\sim 10^{-1}$, far outside the EMRI regime; we refer to Ref.
Albertini:2022rfe  for a thorough accuracy analysis.

It is crucial that EMRI waveforms are not only accurate but also fast and
computationally efficient, as typical Bayesian analysis requires $\sim
10^{5-6}$ waveform evaluations to infer the parameters that govern the
underlying true model. The vast parameter space of EMRIs (given by fourteen
parameters, excluding the small companion spin) further adds to the complexity
and computational burden of Bayesian inference. In order to develop approaches
to tackle the huge parameter space, and because self-force models were not
available, several fast-to-evaluate but approximate “kludge” models
barack2004lisa ; babak2007kludge ; chua2017augmented  were developed. Despite
not being as faithful as self-force waveforms, kludge models were essential
for scoping out early LISA science. These models have been commonly employed
in parameter estimation studies in conjunction with very efficient, but
unfortunately non-robust systematic tests. Such systematic studies include the
Lindblom criteria lindblom2008model , mismatches between waveforms, and the
requirement that the orbital dephasing between two different trajectories is
$\lesssim$ 1 radian. More sophisticated data analysis studies on EMRIs
employed Fisher Matrix-based estimates to quote precision measurements of
parameters vallisneri2008use ; Babak:17aa ; Piovano:2021iwv ; burke2021extreme
; Piovano:2022ojl ; gupta2022modeling ; gair2013testing ;
moore2017gravitational ; maselli2022detecting  and biases on parameters
arising from waveform modeling errors cutler2007lisa ; speri2021assessing .
Fisher matrices are cheap to compute, requiring a tiny number of waveform
evaluations compared to Bayesian inference. However, they are prone to severe
numerical instabilities vallisneri2008use , and they assume that the
underlying probability distribution of the parameters is approximately
Gaussian. The use of Bayesian inference is not a silver bullet either: un-
converged posteriors can result in false conclusions, which can turn out to be
quite dangerous when forecasting LISA science. However, if performed with
care, the results can be interpreted as more robust than other systematic
tests discussed in this paragraph.

Waveform models based on self-force calculations are now emerging. Typically,
these require an expensive offline step involving the calculation of the self-
force or other post-adiabatic effects, and a more rapid online step that
computes the inspiral trajectory and the associated waveform. The first
waveform models that contained partial post-adiabatic phasing information took
minutes to hours to evaluate Osburn:2015duj ; Warburton:2017sxk ;
Drummond:2023loz , but more recently near-identity averaging
VanDeMeent:2018cgn ; Lynch:2021ogr ; Lynch:2023gpu  and two-timescale
approaches Miller:2020bft  have reduced the calculation of the inspiral
trajectory to a few seconds. Rapidly generating waveforms that include the
full mode content then requires additional optimizations katz2021fast ;
chua2021rapid .

The work presented here is a first of its kind. It is the first waveform
accuracy study to assess the requirement of post-adiabatic terms for LISA-
focused studies based entirely on Bayesian inference. We employ the state-of-
the-art FastEMRIWaveforms (FEW) framework katz2021fast ; chua2021rapid , which
exploits the EMRI multiscale structure to quickly generate sub-second
waveforms with accuracy suitable for LISA data analysis. We also incorporate
state-of-the-art second-order GSF results into the FEW framework. Moreover,
our analysis includes the latest LISA response function katz2022assessing ,
yielding second-generation Time-Delay Interferometry (TDI) variables with
realistic orbits generated by the European Space Agency martens2021trajectory
. Finally, we also use the most recent instrumental noise model given by the
LISA mission requirements team LISAsr:18aa . All computations presented in
this work exploit Graphical Processing Units (GPUs) to accelerate the waveform
and LISA response evaluation time to $\lesssim$ 1 second, making Bayesian
inference possible. Thus, we employ the most accurate EMRI waveforms present
in the literature to date. Armed with these tools, our goal is to understand
the importance of the first post-adiabatic corrections when performing
parameter estimation. Finally, we extended our analysis for a class of
intermediate mass-ratio inspirals (IMRIs), in particular assuming the primary
is a massive black hole (MBH) with mass $M\sim 10^{5-7}M_{\odot}$ (like an
EMRI) while the smaller companion has mass $\mu\sim 10^{3}M_{\odot}$222We
choose these masses for an IMRI to ensure that the corresponding GW signal is
within the LISA band. In general, both components of an intermediate mass-
ratio binary may be much lighter mandel2009can ; mandel2008rates ..

The paper is organized as follows: in Sec. II we describe the approximate and
exact waveform models used throughout this work; Sec. III outlines the data
analysis schemes. The results are presented in Sec IV, more specifically Sec.
IV.1 on mismodeling post-adiabatic templates and Sec. IV.2 on constraining the
spin of the smaller companion. Finally, we discuss the impact of mismodeling
adiabatic templates for eccentric orbits in Sec. IV.3. A comprehensive summary
of the results is given in Sec.V. Finally, we present a discussion alongside
the scope for future work in Sec. VI and Sec. VII respectively.

Expert readers who do not wish to dig through the details of the paper can
instead skip straight to the main plots, namely Fig. 1, circular orbit corner
plots Figs. 2–4 and eccentric orbit corner plots Figs. 5–6. We use $G=c=1$
units throughout the paper unless otherwise stated.

##  II Waveform models

We assume that the central body is a spinless BH described by the
Schwarzschild metric with line element

$\displaystyle ds^{2}$
$\displaystyle=-\bigg{(}1-\frac{2M}{r}\bigg{)}dt^{2}+\bigg{(}1-\frac{2M}{r}\bigg{)}^{-1}dr^{2}$
$\displaystyle\quad+r^{2}\left(d\theta^{2}+\sin\theta^{2}d\phi^{2}\right)\,.$
(2)

The binary’s smaller companion is a generic CO endowed with spin; we do not
make any assumption on its internal composition.

To conform with Ref. Wardell:2021fyy , we find it convenient in our
computations to use the symmetric mass ratio
$\nu=\mu/(M+\mu)^{2}=\epsilon/(1+\epsilon)^{2}=\epsilon+\mathcal{O}(\epsilon^{2})$.
The symmetric mass ratio ranges from $0<\nu\leq 1/4$ and has been shown to
provide a better agreement for binding energy, fluxes, and waveforms when
compared with numerical relativity simulations for comparable mass binaries
Pound:2019lzj ; Warburton:2021kwk ; Wardell:2021fyy ; vandeMeent:2020xgc .
However, it does not materially affect accuracy at the mass ratios we consider
here. We implement in FEW two 1PA waveforms specialized to quasicircular
orbits. The first one is a hybrid, approximate waveform, where the 0PA fluxes
were computed exactly (up to negligible numerical error) using linear BH
perturbation theory while the post-adiabatic second-order self-force
corrections are obtained by resumming a third-Post-Newtonian-order (3PN)
expansion. We label this model as cir0PA+1PA-3PN. The second 1PA waveform is a
state-of-the-art model that includes all relevant 1PA terms: the second-order
self-force fluxes and binding energy corrections, computed in Ref.
Warburton:2021kwk  and Ref. Pound:2019lzj , respectively, and the
contributions from the secondary spin, given in Refs. Piovano:2021iwv ;
Akcay:2019bvk . This second waveform model is labeled as cir1PA. Finally, to
allow for comparison between 0PA and 1PA waveform models, we implement an
adiabatic template that we call cir0PA. This is identical to the cir1PA model
but with the post-adiabatic terms removed and setting the spin of the
secondary to zero.

To assess the potential impact of eccentricity on EMRI data analysis, we
additionally compare two other adiabatic models: the fully relativistic 0PA
waveform available in the FEW package of the BHPToolkit BHPToolkit  and an
approximate waveform that includes GW fluxes known at 9PN Munna:2020juq . We
label the former as ecc0PA and the latter as ecc0PA-9PN. The two models have
the same waveform amplitudes and geodesic orbital frequencies and differ only
in their expressions for the energy and angular momentum fluxes that drive the
inspiral. For more details on the ecc0PA model, we refer the reader to Refs.
Chua:2020stf ; Katz:2021yft , while we provide the evolution equations for the
ecc0PA-9PN model in Sec. II.2.

In the following subsections, we summarize the cir1PA model and describe the
approximations introduced in cir0PA+1PA-3PN and ecc0PA-9PN. For all waveform
models, we follow the convention adopted in Ref. Katz:2021yft  for the source
reference frame, namely the orbital angular momentum is set to be aligned to
the z-axis of a Cartesian frame centered on the MBH. In this way, the motion
is confined to the equatorial plane $\theta=\pi/2$, and the orbital plane does
not precess.

###  II.1 Quasicircular inspirals with spinning secondary to 1PA order

The equations of motion for EMRIs can be conveniently expanded in the small
mass-ratio $\epsilon$. At zeroth order, the motion of the binary components is
approximated by a free-falling point-particle in a fixed background spacetime.
Post-geodesic corrections arise from the self-force (see Refs. Pound:2021qin ;
Harte:2014wya  for a comprehensive review) and the coupling between the
curvature and spin of the smaller companion, called spin-curvature force
Tanaka:1996ht ; Dixon:1964NCim ; Dixon:1970I ; Dixon:1970II . The model we
consider is specifically based on a two-timescale formulation of this
small-$\epsilon$ expansion. See Refs. Hinderer:2008dm ; Miller:2020bft  and
Ref. Pound:2021qin  for a comprehensive review. As the name suggests, a two-
timescale expansion assumes the existence of two disparate timescales. The
short one is the orbital timescale, i.e. the orbital period of the secondary.
The long timescale is the radiation-reaction timescale. During the inspiral,
orbital quantities like the frequency, radius, and waveform amplitudes evolve
on the slow timescale, at a rate of order $\mathcal{O}(\epsilon)$. By
contrast, orbital phases evolve on the fast timescale, at a rate of order
$\mathcal{O}(1)$.

In this section we summarize the orbital evolution in this formulation through
1PA order, dividing the discussion into conservative corrections and
dissipative evolution.

####  II.1.1 Conservative corrections to orbital motion

We first summarize the conservative corrections to the (slowly varying)
constants of motion in the case of quasicircular, equatorial orbits in the
Schwarzschild spacetime. For more details on the interplay between the self-
force and spin-curvature force, see Ref. Mathews:2021rod , while for more
details on the dynamics of spinning particles in circular equatorial orbits,
see Tanaka:1996ht ; Piovano:2020zin .

In general, if dissipation is neglected, a spinning particle in Schwarzschild
spacetime admits four conserved quantities, which are the (normalized) energy
$\check{E}=E/\mu$ and the components of the (normalized) total angular
momentum $\vec{J}/(\mu M)=(J_{x},J_{y},J_{z})/(\mu M)$ Rudiger:1981a ;
Rudiger:1981b . With our choice of reference frame, $\vec{J}/(\mu
M)=(0,0,J_{z})/(\mu M)$, we define accordingly $\check{J}=J_{z}/(\mu M)$. As
initial conditions, we choose $\theta=\pi/2$ and set the secondary spin
(anti-)aligned to the $z$-axis, which implies that $\chi>0$ is aligned
($\chi<0$ is anti-aligned) to $\vec{J}$, where $\chi$ is the secondary’s
dimensionless spin parameter. Such conditions ensure that neither the orbital
plane nor the secondary spin precess Tanaka:1996ht ; Witzany:2023bmq .

Hereafter, hatted quantities refer to dimensionless variables normalized by
$M$ (as opposed to quantities with checks, such as $\check{E}$, which are
normalized with $\mu$). For example, $\widehat{\Omega}_{\phi}=M\Omega_{\phi}$
and $\hat{r}=r/M$. The binding energy $\check{E}$ is the only first integral
we need to model the inspiral at 1PA order for our quasicircular orbital
configurations. It is convenient to parameterize the orbit in terms of its
orbital frequency $\widehat{\Omega}_{\phi}$. The orbital radius of the
particle, $\hat{r}_{p}$, can then be written to linear order in $\nu$, at
fixed frequency $\widehat{\Omega}_{\phi}$, as

$\hat{r}_{p}=\hat{r}+\nu\chi\delta\hat{r}^{\chi}(\widehat{\Omega}_{\phi})+\nu\delta\hat{r}^{\text{SF}}(\widehat{\Omega}_{\phi}).$
(3)

The leading-order term is the _geodesic_ orbital radius,
$\hat{r}=\widehat{\Omega}_{\phi}^{-2/3}$, while the corrections
$\delta\hat{r}^{\chi}(\widehat{\Omega}_{\phi})$ and
$\delta\hat{r}^{\text{SF}}(\widehat{\Omega}_{\phi})$ are the linear shifts due
to the secondary spin and to the conservative first-order self-force, given by

$\delta\hat{r}^{\chi}=-\frac{1}{\sqrt{\hat{r}}}$ (4)

and by Eq. (A8) of Ref. Miller:2020bft , respectively. We can similarly expand
$\check{E}$ to first order in $\nu$ at fixed frequency:

$\displaystyle\check{E}$
$\displaystyle=\check{E}_{0}+\nu\chi\check{E}_{1}^{\chi}+\nu\check{E}_{1}^{\text{SF}}\
,$ (5)

where $\check{E}_{0}$ is the geodesic binding energy and
$\check{E}_{1}^{\chi}$ is the shift induced by the secondary spin
Jefremov:2015gza , given by

$\displaystyle\check{E}_{0}=\frac{\hat{r}-2}{\sqrt{\hat{r}}\sqrt{\hat{r}-3}}-1\
,\qquad\check{E}_{1}^{\chi}=-\frac{1}{\hat{r}^{2}\sqrt{\hat{r}-3}}\ .$ (6)

Finally, $\check{E}_{1}^{\text{SF}}$ is the correction to the binding energy
induced by the conservative piece of the first-order self-force LeTiec:2011dp
(see Refs. Pound:2019lzj ; Albertini:2022rfe  for discussion).

####  II.1.2 Evolution equations

We now consider the binary evolution through 1PA order, accounting for the
slow evolution of due to dissipation.

In our models, we neglect the evolution of the primary mass and spin and
second-order GSF horizon fluxes, the latter of which are currently unknown.
While these effects appear formally at 1PA order, they have been shown to have
a numerically small effect compared to the overall 1PA contribution to the
inspiral Martel:2003jj ; Hughes:2018qxz ; Albertini:2022rfe , so we can safely
neglect them.

The azimuthal orbital phase $\Phi_{\phi}$ and the leading-order orbital radius
$\hat{r}$ evolve according to the following coupled set of equations:

$\displaystyle\frac{\mathrm{d}\Phi_{\phi}}{\mathrm{d}\hat{t}}$
$\displaystyle=\widehat{\Omega}_{\phi}(\hat{r})\ ,$ (7)
$\displaystyle\frac{\mathrm{d}\hat{r}}{\mathrm{d}\hat{t}}$
$\displaystyle=-\nu\big{[}F_{0}(\hat{r})+\nu F_{1}(\hat{r})\big{]}\ .$ (8)

Here, $F_{0}(\hat{r})$ is the leading-order, adiabatic forcing term given by
$F_{0}(\hat{r})=\big{(}\partial\check{E}_{0}/\partial\hat{r}\big{)}^{-1}\mathcal{F}_{0}(\hat{r})$,
where $\mathcal{F}_{0}(\hat{r})$ is the leading-order energy flux. The
subleading, 1PA term $F_{1}(\hat{r})$ is what differs between our two circular
orbit models. For the complete 1PA model, cir1PA, we have

$F_{1}(\hat{r})=F_{1}^{\text{SF}}(\hat{r})+\chi F_{1}^{\chi}(\hat{r}).$ (9)

where

$\displaystyle F_{1}^{\text{SF}}(\hat{r})$
$\displaystyle=a(\hat{r})\mathcal{F}_{1}(\hat{r})+a(\hat{r})^{2}\bigg{(}\frac{\partial\check{E}_{1}^{\text{SF}}}{\partial\hat{r}}\bigg{)}\mathcal{F}_{0}(\hat{r})\
,$ (10) $\displaystyle F_{1}^{\chi}(\hat{r})$
$\displaystyle=a(\hat{r})\mathcal{F}_{1}^{\chi}(\hat{r})+a(\hat{r})^{2}\bigg{(}\frac{\partial\check{E}_{1}^{\chi}}{\partial\hat{r}}\bigg{)}\mathcal{F}_{0}(\hat{r})\
,$ (11)

and $a(\hat{r})=\big{(}\partial\check{E}_{0}/\partial\hat{r}\big{)}^{-1}$. The
post-adiabatic corrections to the fluxes
$\mathcal{F}_{1}^{\text{SF}}(\hat{r})$ and $\mathcal{F}_{1}^{\chi}(\hat{r})$
are, respectively, generated by nonlinear (quadratic) terms in the field
equations Warburton:2021kwk  and by the secondary spin Piovano:2021iwv .

In the cir0PA+1PA-3PN model, $F_{1}(\hat{r})$ is approximated as

$F_{1}(\hat{r})=F_{1}^{3\text{PN}}(\hat{r})+\chi F_{1}^{\chi}(\hat{r}),$ (12)

where $F_{1}^{3\text{PN}}$ has the same functional form as Eq. (10), but it is
constructed using analytic post-Newtonian expressions for the fluxes at 0PA
and 1PA order, including the self-force correction to the binding energy, in
which we retain relative 3PN order accuracy for each. The adiabatic and post-
adiabatic fluxes can be extracted by taking the leading and next-to-leading-
order-in-$\nu$ terms in the 4PN fluxes recently computed in Ref.
Blanchet:2023sbv . The binding energy is computed using post-Newtonian self-
force expansions Fujita:2012cm . Explicitly, these are given in App. B. It is
crucially important that we treat the 3PN flux terms and binding energy as
polynomial approximants when used in Eq. (10), without further expanding the
geodesic quantity $a(\hat{r})$. This is because $a(\hat{r})$ is divergent at
the lightring, rendering its large-radius Taylor expansion extremely
inaccurate in the strong field, which would corrupt the entire model.

The evolution equations for our adiabatic model are obtained by simply setting
the forcing function $F_{1}(\hat{r})$ to zero in Eqs. (7-8). We label this
model as cir0PA. Finally, we remark that for all mass-ratios considered in
this work, each inspiral terminates far (in radial coordinate distance) from
where the transition to plunge begins ori2000transition ; o2003transition ;
kesden2011transition ; apte2019exciting ; burke2020transition ;
compere2020transition ; compere2021self ,
$\hat{r}-\hat{r}_{\text{isco}}\sim\epsilon^{2/5}$. This is important since our
two-timescale evolution equations cease to accurately describe the trajectory
around that point Albertini:2022rfe .

###  II.2 Non-spinning eccentric 0PA inspirals

In the eccentric case, the orbital radius $\hat{r}$ oscillates on the orbital
timescale between a maximum value $\hat{r}_{\text{max}}$ and a minimum value
$\hat{r}_{\text{min}}$, meaning it is not an ideal quantity to evolve
directly. Instead, we parameterize the eccentric orbit using the slowly
evolving semi-latus rectum $p$ and eccentricity $e$, which are defined in
terms of the slowly evolving maximum and minimum values of the orbital radius:

$\displaystyle
p=\frac{2\hat{r}_{\text{max}}\hat{r}_{\text{min}}}{(\hat{r}_{\text{max}}+\hat{r}_{\text{min}})}\
,\qquad
e=\frac{\hat{r}_{\text{max}}-\hat{r}_{\text{min}}}{\hat{r}_{\text{max}}+\hat{r}_{\text{min}}}\
.$ (13)

Note that in the circular orbit limit, $e\rightarrow 0$ and
$p\rightarrow\hat{r}$.

Due to the orbital-timescale radial motion, the waveform picks up an
additional phase $\Phi_{r}$ which evolves with the Boyer-Lindquist fundamental
radial frequency $\widehat{\Omega}_{r}$. Expressions for the fundamental
frequencies $\widehat{\Omega}_{\phi}(p,e)$ and $\widehat{\Omega}_{r}(p,e)$ for
eccentric geodesic orbits in Schwarzschild spacetime can be found in Ref.
Cutler:1994pb . In practice we make use of the semi-analytic expressions in
terms of elliptic integrals in Refs. Fujita:2009bp  which have been
implemented in the KerrGeodesics package KerrGeodesicsPackage .

The eccentric equations of motion at adiabatic order can be written as

$\displaystyle\begin{split}\frac{\mathrm{d}\Phi_{\phi}}{\mathrm{d}\hat{t}}&=\widehat{\Omega}_{\phi}(p,e)\
,\end{split}$ (14)
$\displaystyle\begin{split}\frac{\mathrm{d}\Phi_{r}}{\mathrm{d}\hat{t}}&=\widehat{\Omega}_{r}(p,e)\
,\end{split}$ (15)
$\displaystyle\begin{split}\frac{\mathrm{d}p}{\mathrm{d}\hat{t}}&=-\epsilon\left(\frac{\partial
p}{\partial\check{E}_{0}}\mathcal{F}_{0}^{E}(p,e)+\frac{\partial
p}{\partial\check{J}_{0}}\mathcal{F}_{0}^{J}(p,e)\right),\end{split}$ (16)
$\displaystyle\begin{split}\frac{\mathrm{d}e}{\mathrm{d}\hat{t}}&=-\epsilon\left(\frac{\partial
e}{\partial\check{E}_{0}}\mathcal{F}_{0}^{E}(p,e)+\frac{\partial
e}{\partial\check{J}_{0}}\mathcal{F}_{0}^{J}(p,e)\right),\end{split}$ (17)

where $\mathcal{F}_{0}^{E}$ and $\mathcal{F}_{0}^{J}$ are the leading-order
total fluxes of energy and angular momentum, respectively. We switch back to
the small mass ratio $\epsilon$ in the above equations since both the vanilla
FEW and ecc0PA-9PN models are implemented in terms of $\epsilon$ instead of
the symmetric mass ratio $\nu$. By employing the geodesic relations for
$\check{E}_{0}(p,e)$ and $\check{J}_{0}(p,e)$, i.e., Cutler:1994pb

$\displaystyle\check{E}_{0}(p,e)$
$\displaystyle=\sqrt{\frac{(p-2)^{2}-4e^{2}}{p(p-e^{2}-3)}}\ ,$ (18)
$\displaystyle\check{J}_{0}(p,e)$ $\displaystyle=\frac{p}{\sqrt{p-e^{2}-3}}\
,$ (19)

one can obtain the various partial derivatives in Eqs. (14 \- 17) by
constructing a Jacobian and inverting so that

$\displaystyle\begin{pmatrix}\frac{\partial
p}{\partial\check{E}_{0}}&\frac{\partial p}{\partial\check{J}_{0}}\\\
\frac{\partial e}{\partial\check{E}_{0}}&\frac{\partial
e}{\partial\check{J}_{0}}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\frac{\partial\check{E}_{0}}{\partial
p}&\frac{\partial\check{E}_{0}}{\partial e}\\\
\frac{\partial\check{J}_{0}}{\partial p}&\frac{\partial\check{J}_{0}}{\partial
e}\end{pmatrix}^{-1}$ (20)
$\displaystyle=\frac{1}{\frac{\partial\check{E}_{0}}{\partial
p}\frac{\partial\check{J}_{0}}{\partial
e}-\frac{\partial\check{E}_{0}}{\partial
e}\frac{\partial\check{J}_{0}}{\partial
p}}\begin{pmatrix}\frac{\partial\hat{J}_{0}}{\partial
e}&-\frac{\partial\check{E}_{0}}{\partial e}\\\
-\frac{\partial\hat{J}_{0}}{\partial p}&\frac{\partial\check{E}_{0}}{\partial
p}\end{pmatrix}\ .$ (21)

Explicit expressions are available in Ref. Hughes:2021exa . Note that this
procedure will give a singular expression for the rate of change of $e$ in the
circular orbit limit, $e=0$. However, it is well established that for
adiabatic inspirals Kennefick:1995za , quasicircular orbits remain
quasicircular and so we can simply use $de/d\hat{t}=0$ in this special case.

The only difference between the two eccentric models is the expressions for
the total energy and angular momentum fluxes. For the ecc0PA model, we
numerically solve the Teukolsky equations for the energy and angular momentum
fluxes to infinity and down the horizon generated by a point particle on a
geodesic orbit with a given value of $p$ and $e$. We do this for multiple
points in the $(p,e)$ parameter space and interpolate using splines so that
the fluxes can be rapidly evaluated when solving the eccentric equations of
motion. For more details, see Sec. IV A of Ref. katz2021fast .

By contrast, the ecc0PA-9PN model uses analytic PN expansions of the fluxes of
energy and angular momentum to infinity. These are valid to leading order in
the mass ratio and relative 9PN order Munna:2020juq , and are available in the
PostNewtonianSelfForce package of the BHPToolkit PostNewtonianSelfForcePackage
. The PN series contain expansions in eccentricity to at least order $e^{10}$,
with some PN components having even higher order expansions up to $e^{30}$.

We highlight that the cir0PA model and ecc0PA model will not agree in the
limit as $e\rightarrow 0$ because the former is expanded in $\nu$ whereas the
latter is expanded in $\epsilon$. For our purposes, this is not a problem: we
will only directly compare waveforms expanded in the same small perturbative
variable, e.g., cir0PA with cir1PA ($\nu$) and ecc0PA with ecc0PA-9PN
($\epsilon$). We checked that recovery of a cir0PA model with an ecc0PA model
(at $e=0$) yields a bias on the secondary mass $\mu$. In turn, such bias
results in a biased mass ratio $\epsilon_{\text{bias}}$ given by

$\epsilon_{\text{bias}}\approx\epsilon_{\text{true}}-2\epsilon_{\text{true}}^{2}.$
(22)

Since
$\nu_{\text{true}}=\epsilon_{\text{true}}/(1+\epsilon_{\text{true}})^{2}\sim\epsilon_{\text{true}}-2\epsilon_{\text{true}}^{2}+\mathcal{O}(\epsilon^{3})$,
we recovered the symmetric mass-ratio, as expected. If the 1PA components were
known for eccentric orbits, then it would not matter whether the equations of
motion were expanded in $\epsilon$ or $\nu$ (for the range of mass ratios we
consider).

###  II.3 Waveform models

In the Teukoslky formalism, the gravitational waveform detected by an observer
at infinity in the Schwarzschild spacetime can be written as Hughes:2021exa

$\displaystyle h_{+}-ih_{\times}$ $\displaystyle=\frac{\mu}{D_{\rm
S}}\displaystyle\sum_{\ell,m,n}\mathcal{A}_{\ell mn}(p(t),e(t))Y_{\ell
m}(\vartheta,\varphi)e^{-i\Phi_{mn}(t)}\,,$ (23)

where $D_{\rm S}$ is the source’s luminosity distance from the detector,
$\mathcal{A}_{\ell mn}\equiv 2\hat{Z}^{\infty}_{\ell
mn}/\hat{\omega}^{2}_{mn}$, and $\hat{Z}^{\infty}_{\ell
mn}=M^{2}Z^{\infty}_{\ell mn}$, the latter being the inhomogenous solution of
the radial Teukoslky equation in the limit $\hat{r}\to\infty$. Here the GW
frequency $\hat{\omega}_{mn}$ is defined as
$\hat{\omega}_{mn}=m\widehat{\Omega}_{\phi}+n\widehat{\Omega}_{r}$. In
practice, the mode amplitudes $\mathcal{A}_{\ell mn}$ are interpolated across
the $p-e$ space using a neural network chua2021rapid ; katz2021fast . To
determine the waveform mode content, the modes are cumulatively summed in
order of decreasing power. This summation is truncated at a threshold value of
$(1-\kappa)$ times the total power of the modes, with $\kappa$ a tunable
parameter. Only the modes that pass this threshold are included in the
waveform computation (for further details see chua2021rapid ; katz2021fast ).

For our circular orbit runs, we have checked that the mode content of injected
waveforms and model template waveforms are identical. This is to ensure
consistency between waveforms and to make sure that biases are not a feature
of missing important modes between waveform models in the analysis. When
generating waveforms, we use the default mode selection parameter
$\kappa=10^{-5}$ in FEW, resulting in the same 12 $\ell m$-modes _for all_
circular orbit-based waveforms present in this work. The number of modes
scales considerably with increasing eccentricity (see Fig. 4 in katz2021fast
); this will be discussed in Sec. IV.3.

Since our waveform represents an evolving source, we do not have a
decomposition into discrete frequencies, but a multi-voice decomposition where
the evolving “voices” are found by solving the equations of motion for the
phases $\Phi_{\phi}(t)$ and $\Phi_{r}(t)$ and summing together as
$\Phi_{mn}(t)=m\Phi_{\phi}(t)+n\Phi_{r}(t)$.

The functions $Y_{\ell m}(\vartheta,\varphi)$ appearing in Eq. (23) are the
spin-weighted spherical harmonics for spin weight $s=-2$ Goldberg:1967 , while
the angles $(\vartheta,\varphi)$ identify the direction of the detector in the
source reference frame. Due to the LISA constellation’s motion, the sources’
sky positions are measured with respect to a suitable solar system barycentric
(SSB) frame attached to the ecliptic Barack:2003fp . We adopt the same
convention as Ref. Katz:2021yft  and label the binary’s sky position and the
direction of binary’s angular momentum as $(\theta_{S},\phi_{S})$ and
$(\theta_{K},\phi_{K})$, respectively. The viewing angle $\varphi$ is set to
$\varphi=-\pi/2$, which implies that $\cos\vartheta$ can then be written in
terms of the constant angles $(\theta_{S},\phi_{S})$ and
$(\theta_{K},\phi_{K})$ as

$\cos\vartheta=-\cos\theta_{S}\cos\theta_{K}-\sin\theta_{S}\sin\theta_{K}\cos(\phi_{S}-\phi_{K})\
.$ (24)

The initial phases of the waveform are then entirely determined by the angles
$\Phi_{\phi_{0}}$ and $\Phi_{r_{0}}$. Notice that the angles
$(\theta_{S},\phi_{S})$ and $(\theta_{K},\phi_{K})$ also appear in the LISA
response function katz2022assessing  and katz2021fast .

####  II.3.1 1PA waveforms for circular orbits

In the circular orbit case, we no longer have a radial phase and its
corresponding harmonic contribution, so the waveform model simplifies to

$\displaystyle h_{+}-ih_{\times}$ $\displaystyle=\frac{\mu}{D_{\rm
S}}\displaystyle\sum_{\ell,m}\mathcal{A}_{\ell m}(\hat{r}(t))Y_{\ell
m}(\vartheta,\varphi)e^{-im\Phi_{\phi}(t)}\ .$ (25)

Note that at 1PA order, there will also be order-mass-ratio corrections to the
waveform amplitudes $\mathcal{A}_{\ell m}=\mathcal{A}_{\ell
m}^{(0)}+\epsilon\mathcal{A}_{\ell m}^{(1)}+\mathcal{O}(\epsilon^{2})$. GW
interferometers are much more sensitive to fluctuations in frequency, rather
than amplitude, especially for asymmetric binaries with $\epsilon\ll 1$. As
such, including sub-leading corrections to the orbital phase is more important
than to the amplitudes. Since the largest mass ratio we consider is
$\epsilon\sim 10^{-3}$, we can safely neglect 1PA corrections to the
amplitudes in our parameter estimation studies. For this reason, we employed
Eq. (25) for our all circular orbit models: cir1PA, cir0PA+1PA-3PN and cir0PA
with amplitudes computed at adiabatic order.

##  III LISA data analysis

GW data analysis relies on three crucial ingredients: GW waveforms, a
description of the time-evolving instrument response to the incoming waves,
and an accurate characterization of the noise. In this section, we describe
our model for the LISA response function and the noise process, ultimately
leading to the Whittle-likelihood.

###  III.1 The data stream of LISA

The projection of the polarised GW signals onto the LISA arms depends on the
geometry of the instrument, which changes in time due to the spacecraft’s
motion. The LISA response function is then a dynamical quantity in both the
time and frequency domain, and clearly much more complicated than the “static”
response of ground-based detectors. Such a feature introduces severe
complexities to the accurate modeling and efficient computation of the LISA
response function. By projecting the incoming GWs onto the arm-lengths of the
detector, one can model the deformations across the six LISA links between the
individual craft. It is then possible to construct a first set of time-shifted
second-generation TDI variables, which massively suppress the laser noise (by
$\sim 8$ orders of magnitude). These variables can be linearly combined to
construct a further set of TDI variables $(A,E,T)$ armstrong1999time ;
Tinto:2004wu , which are uncorrelated in their noise properties. In our work,
we use the TDI variables $(A,E,T)$ Tinto:2004wu ; vallisneri2005synthetic .

The data streams can then be written as

$d^{(X)}(t)=h_{\rm e}^{(X)}(t;\bm{\theta}_{\text{tr}})+n^{(X)}(t)\,,\quad
X=\\{A,E,T\\},$ (26)

where $h^{(X)}_{\rm e}$ denotes the true deterministic signal with true
parameters $\bm{\theta}_{\text{tr}}$ for the $X$ TDI observable. We use lisa-
on-gpu, a framework to compute the response of the LISA instrument to the
incoming GWs available at katz2022assessing , which generates the three TDI
variables. We employ the most realistic orbits of the LISA crafts generated by
the European Space Agency (ESA) martens2021trajectory . These numerically
computed orbits take into account gravitational attractions from relevant
celestial bodies in our solar system and are built to minimize the “breathing”
of the LISA interferometer arms. Such orbits have _approximately_ equal and
constant arm-lengths, resulting in minor correlations between noise components
vallisneri2012non  that we neglect in favor of computational efficiency
katz2022assessing . If these correlations are mismodeled, then the injected
mismodeled noise realizations will impact the statistical errors of the
recovered parameters burke2021extreme ; talbot2021inference ; edy2021issues .
This is not a problem in our analysis because, as we will explain later on, we
consider zero-noise injections. Furthermore, the signal-to-noise ratio is
weakly impacted when mismodeling the noise process333This was demonstrated in
the case of gaps, where correlations between noise components are severe and
cannot be neglected. For more information, see Chap. 8 in burke2021extreme ..

###  III.2 Noise process and likelihood

We assume that the noise in each channel is a zero-mean, weakly stationary,
ergodic, Gaussian random process with noise covariance matrix (in the domain
of positive frequencies) wiener1930generalized ;
khintchine1934korrelationstheorie

$\langle\tilde{n}^{(X)}(f)[\tilde{n}^{(X)}(f^{\prime})]^{\star}\rangle=\frac{1}{2}\delta(f-f^{\prime})S^{(X)}_{n}(f^{\prime})\,,$
(27)

with tilded quantities representing Fourier transforms of the time-domain
data, i.e.,

$\tilde{a}(f)=\int_{-\infty}^{\infty}\text{d}t\,a(t)e^{-2\pi ift}\,.$ (28)

In Eq. (27), $\langle\cdot\rangle$ denotes an ensemble averaging process,
$\delta$ is the Dirac delta function and $S^{(X)}_{n}(f)$ is the (one-sided)
power spectral density (PSD) of the instrumental noise process for each
channel $X=\\{A,E,T\\}$. In our analysis, we use the latest SciRDv1 model
noise PSD LISAScienceRequirementsDocument  and second-generation TDI variables
with a python implementation available in KatzTools . Finally, we assume that
the PSDs for each channel $X=\\{A,E,T\\}$ are all known and fixed.

Furthermore, we neglect noise correlations between each TDI channel, i.e.,
$\langle\hat{n}^{X}(f)[\hat{n}^{Y}(f^{\prime})]^{\star}\rangle=0$ for $X\neq
Y$ armstrong1999time ; vallisneri2005synthetic ; katz2022assessing . Under our
assumptions, we can write the Whittle-likelihood for a known PSD
$S^{(X)}_{n}(f)$ as whittle:1957 ; finn1992detection

$\log
p(d|\bm{\theta})\propto-\frac{1}{2}\sum_{X=\\{A,E,T\\}}(d-h_{m}|d-h_{m})_{(X)}$
(29)

with approximate model templates $h^{(X)}_{m}$ and noise-weighted inner
product $(a|b)_{(X)}$ given by finn1992detection

$(a|b)_{X}=4\mathrm{Re}\int_{0}^{\infty}\frac{\tilde{a}^{(X)}(f)(\tilde{b}^{(X)})^{\star}(f)}{S^{(X)}_{n}(f)}\,\text{d}f.$
(30)

Given a model template $h_{m}$, we define the _effective_ SNR with respect to
the true signal $h_{e}$ as

$\rho^{\text{eff}}_{AET}(\bm{\theta})=\left[\sum_{X=\\{A,E,T\\}}\frac{(h_{e}|h_{m}(\bm{\theta}))^{2}_{(X)}}{(h_{m}(\bm{\theta})|h_{m}(\bm{\theta}))_{(X)}}\right]^{1/2},$
(31)

with exact templates $h_{e}$ evaluated at the true parameters
$\bm{\theta}_{\text{tr}}$. The maximum of Eq. (31) is given when
$\bm{\theta}=\bm{\theta}_{\text{tr}}$ and there are no mismodeling errors,
i.e., $h_{m}\equiv h_{e}$ for all $\bm{\theta}$. This maximum is the _optimal_
matched filtering SNR, given by woodward2014probability

$\rho^{\text{opt}}_{AET}=\left[\sum_{X=\\{A,E,T\\}}(h_{e}|h_{e})_{(X)}\right]^{1/2},$
(32)

where $h_{e}$ is evaluated at the true parameters $\bm{\theta}_{\text{tr}}$.
The optimal SNR over each TDI stream $X=\\{A,E,T\\}$ denotes the average power
of the signal when compared to the root-mean-square average of the noise
floor. The greater the SNR of the signal $h_{e}$ in Eq. (26), the greater the
likelihood of claiming detection. Previous works Babak:17aa ;
babak2009algorithm  set, rather arbitrarily, the SNR threshold to claim
detection as $\rho_{AET}\gtrsim 20$. The sources selected in this work are
well above this threshold, with $\rho_{AET}\sim 70$ in the EMRI regime and
$\rho_{AET}\sim 340$ in the IMRI regime.

###  III.3 Bayesian parameter estimation

Parameter estimation in GW astronomy is typically performed using Bayesian
inference. At the heart of Bayesian statistics lies Bayes’ theorem, which, up
to a normalization factor, is given by

$\log p(\bm{\theta}|d)\propto\log p(d|\bm{\theta})+\log p(\bm{\theta})\,.$
(33)

On the right-hand side, $p(d|\bm{\theta})$ is the likelihood function, a
probability distribution that describes the probability of observing the data
stream given the parameters. Under our assumptions about the noise, we can use
the Whittle-likelihood in Eq. (29). The density $p(\bm{\theta})$ is the prior
probability distribution representing a-priori beliefs on the parameter set
$\bm{\theta}$ before observing the data stream $d$. We opt for uninformative,
uniform prior distributions for $\bm{\theta}$. Finally, the quantity
$p(\bm{\theta}|d)$ is the sought-for posterior distribution reflecting our
beliefs on the parameters after our observations. The goal is then to generate
auto-correlated samples from the posterior density $p(\bm{\theta}|d)$ in order
to compute summary statistics.

Our analysis employs Markov-Chain Monte-Carlo (MCMC) techniques to sample from
the posterior distribution. In particular, we use the eryn sampler
Karnesis:2023ras ; michael_katz_2023_7705496 , which is based on the emcee
Foreman-Mackey:2013  code. Our proposal distribution is the default stretch
proposal goodman2010ensemble . MCMC algorithms for EMRI inference are non-
trivial, so we describe in more details our codes in App. A. In Sec. A.1 we
give a brief overview of the eryn and emcee samplers, discussing their
strengths and weaknesses for the simulations presented in this work.

The initial samples are chosen so that
$\bm{\theta}^{(0)}\approx\bm{\theta}_{\text{tr}}$ since our goal is to
identify potential biases close to source parameters, rather than perform a
search. More details on our choice of starting coordinates and prior choices
(see Eq. (46)) are given in App. A.2. Given a chain of samples
$\bm{\theta}^{(i)}\sim p(\bm{\theta}|d)$, we define our “best-fit parameters”
as the _maximum a posteriori_ (MAP) point estimate
$\bm{\theta}_{\text{bf}}=\text{argmax}_{\bm{\theta}}\\{p(\bm{\theta}|d)\\}$.
The _best-fit_ parameters give the best “match” between the model template and
the observed data stream $d$. In other words, when using non-informative
uniform priors, $\bm{\theta}_{\text{bf}}$ maximizes both the likelihood
function and posterior density. If the recovered parameters
$\bm{\theta}_{\text{bf}}\neq\bm{\theta}_{\text{tr}}$ then the MAP estimate
$\bm{\theta}$ gives a biased point estimate of the true source parameters
$\bm{\theta}_{\text{tr}}$.

We neglect the instrumental noise $n^{(X)}(t)$ present in Eq. (26), which
implies that the Whittle-likelihood reduces to

$\log
p(d|\bm{\theta})\propto-\frac{1}{2}\sum_{X=\\{A,E,T\\}}(h_{e}-h_{m}|h_{e}-h_{m})_{(X)}.$
(34)

In this way, we focus on the impact of biases on the parameters arising due to
waveform mismodeling, which may otherwise be obfuscated by nuisance
statistical fluctuations given by noise realizations. We also neglect the
confusion noise sourced by galactic dwarf binaries Cornish_2017  since it has
not been implemented yet in the latest PSD $S^{(X)}_{n}$ KatzTools ;
LISAScienceRequirementsDocument . If the white-dwarf background were included,
the SNR of the signals would be lower, resulting in wider posteriors on the
parameters. Nevertheless, the impact of the white-dwarf background on the SNR
is marginal, hence it would not significantly affect our key results.

###  III.4 Waveform systematics and detection

We now outline our systematic tests used to compare waveform models. We begin
first by describing our main systematic test defined in a Bayesian framework.
We will then describe alternate statistical tests that are present in the
literature.

For a probability density $p(\bm{\theta}|d)$, we define the 68% credible set,
$C_{p(\bm{\theta}|d)}$, of the samples as the probability $P(\bm{\theta}\in
C_{p(\bm{\theta}|d)})=0.68$. Let $\tilde{p}(\bm{\theta}|d)$ represent an
approximate model posterior, generated through inference using approximate
model waveforms. We then define the systematic test

$\mathcal{C}[i]=\begin{cases}1&\theta_{\text{tr}}^{i}\in\hat{C}_{\tilde{p}(\theta^{i}|d)},\\\
0&\text{otherwise},\end{cases}$ (35)

where $\hat{C}_{\tilde{p}(\bm{\theta}^{i}|d)}$ is the estimated _marginalized_
approximate posterior credible interval for parameter $\theta^{i}$. Equation
(35) generalizes the Cutler-Vallisneri CV criterion cutler2007lisa , since it
accounts for non-Gaussian features in the posterior that are not captured by a
Fisher Matrix-based approach. If $\mathcal{C}[i]=(0)1$ for all recovered
parameters $\theta^{i}_{\text{bf}}$, then a waveform model is (un)suitable for
parameter estimation. Put simply, if the true parameter is not contained
within the 68$\%$ credible interval generated using approximate waveforms,
then the waveform model is not suitable for statistical inference. An example
is given in Fig. 7 in App. A. We stress that Eq. (35) is SNR dependent, since
the size of the interval will increase (decrease) with a decrease (increase)
of SNR. Since these analyses are computationally expensive, we adopted
astrophysically motivated SNRs. One can deduce, as a very rough approximation,
how the criterion (35) changes for brighter or dimmer sources by reducing the
size of the credible interval as the SNR increases.

We now describe further quantities that are used to draw comparisons between
waveform models. The overlap function between two waveforms $h^{(X)}_{1}$ and
$h^{(X)}_{2}$ is defined as

$\mathcal{O}(h^{(X)}_{1},h^{(X)}_{2})=\frac{(h_{1}|h_{2})_{(X)}}{\sqrt{(h_{1}|h_{1})_{(X)}(h_{2}|h_{2})_{(X)}}}.$
(36)

Then, over all channels $X=\\{A,E,T\\}$, the total mismatch between two
waveform models $h_{1}$ and $h_{2}$ is

$\displaystyle\mathcal{M}_{AET}(h_{1},h_{2})=1-\sqrt{\frac{1}{3}\sum_{X=\\{A,E,T\\}}\mathcal{O}^{2}(h_{1},h_{2})}\,,$

where $\mathcal{M}=0\ (\mathcal{M}=1)$ indicates that $h_{1}$ and $h_{2}$ are
identical (orthogonal). In a similar way, we define

$\displaystyle\mathcal{M}^{\text{(inj)}}$
$\displaystyle=\mathcal{M}_{AET}(h_{e}(\bm{\theta}_{\text{tr}}),h_{m}(\bm{\theta}_{\text{tr}}))\
,$ (37) $\displaystyle\mathcal{M}^{(\text{bf})}$
$\displaystyle=\mathcal{M}_{AET}(h_{e}(\bm{\theta}_{\text{tr}}),h_{m}(\bm{\theta}_{\text{bf}}))\
.$ (38)

We remark here that (38) is the usual _fitting factor_ , computed after
stochastically identifying parameters $\bm{\theta}_{\text{bf}}$ that minimize
the mismatch function.

It is usual to perform systematic studies between waveform models by analyzing
the difference of orbital phase, called dephasing, between two trajectories of
the CO. We define two types of dephasing in $\Phi_{\phi}$: one between two
trajectories at the injected parameters $\bm{\theta}_{\text{tr}}$ (Eq. (39));
and one between the two trajectories at inferred parameters
$\bm{\theta}_{\text{bf}}$ (Eq. (40)):

$\displaystyle\Delta\Phi^{\text{(inj)}}$
$\displaystyle=\text{Max}\\{\Phi^{\text{exact}}_{\phi}\\}_{\bm{\theta}=\bm{\theta}_{\text{tr}}}-\text{Max}\\{\Phi^{\text{model}}_{\phi}\\}_{\bm{\theta}=\bm{\theta}_{\text{tr}}}\,,$
(39) $\displaystyle\Delta\Phi^{\text{(bf)}}$
$\displaystyle=\text{Max}\\{\Phi^{\text{exact}}_{\phi}\\}_{\bm{\theta}=\bm{\theta}_{\text{tr}}}-\text{Max}\\{\Phi^{\text{model}}_{\phi}\\}_{\bm{\theta}=\bm{\theta}_{\text{bf}}}\
.$ (40)

Eq. (39) is the usual quantity used to estimate the accuracy requirements of
EMRI waveforms. We comment that the two equations above are evaluated over the
same time of observation.

Finally, we refer to the accumulated SNR normalised by the optimal SNR as the
quantity

$\displaystyle\rho^{\text{(inj)}}/\rho^{\text{(opt)}}$
$\displaystyle=\rho_{AET}^{\text{eff}}(\bm{\theta}_{\text{tr}})/\rho_{AET}^{\text{opt}}\,,$
(41) $\displaystyle\rho^{\text{(bf)}}/\rho^{\text{(opt)}}$
$\displaystyle=\rho_{AET}^{\text{eff}}(\bm{\theta}_{\text{bf}})/\rho_{AET}^{\text{(opt)}}\,,$
(42)

for $\rho^{\text{eff}}_{AET}(\bm{\theta})$ and $\rho^{\text{opt}}$ defined in
Eq. (31) and Eq. (32), respectively. The equations above represent the
fraction of SNR (normalised between 0 and 1) accumulated throughout the
inspiral. The quantity in Eq. (42) is useful to determine whether a reference
true waveform model is detectable by an approximate model. The quantities in
Eqs. (37 \- 42) will be useful to compare against our Bayesian inference
results.

##  IV Results

We now present our results on the Bayesian parameter inference with complete
circular 1PA and mismodeled eccentric 0PA waveforms. In Sec. IV.1, we will
investigate the impact of mismodeling EMRI templates by neglecting some (or
all) of the post-adiabatic corrections. In Sec. IV.2, we will then focus our
attention on constraining the secondary spin parameter. Finally, in Sec. IV.3
we will describe the impact of mismodeling 0PA eccentric templates.

The inferred parameters are the following: the red-shifted masses $M$ and
$\mu$, the initial radial coordinate and initial phase $r_{0}/M$ and
$\Phi_{\phi_{0}}$, respectively (both defined at time $t=0$), the luminosity
distance to the source $D_{\rm S}$, the source sky position $(\theta_{\rm
S},\phi_{\rm S})$ and the orientation of the orbital angular momentum
$(\theta_{\rm K},\phi_{\rm K})$. Waveforms are generated in the source frame
and transformed to the SSB frame via the response function. For the sake of
clarity, if a model includes the spin of the secondary as a parameter it will
be denoted “w/ spin” and otherwise “w/o spin”. In Sec. IV.1, we will only
infer the secondary spin $\chi$ using the cir1PA w/ spin model, whereas in
Sec. IV.2 we will infer the secondary spin with all approximate models. For
sections IV.1 and IV.2, configurations for each of the mass ratios and the
relative SNR can be found in Table 1.

Config. | $\epsilon$ | $M\,[M_{\odot}]$ | $r_{0}/M$ | $D_{\rm S}\,[\text{Gpc}]$ | $T_{\text{obs}}\,[\text{yrs}]$ | $\rho_{AET}$
---|---|---|---|---|---|---
(1) | $10^{-5}$ | $10^{6}$ | $10.6025$ | 1.0 | 2.0 | 70
(2) | $10^{-4}$ | $10^{6}$ | $15.7905$ | 2.0 | 1.5 | 65
(3) | $10^{-3}$ | $5\cdot 10^{6}$ | $16.8123$ | 1.0 | 1.0 | 340
Table 1: Here we tabulate injection EMRI/IMRI parameters for all waveforms in
sections IV.1 and IV.2. The IMRI configuration is given by the last row in the
table. For each case, we use identical extrinsic parameters: $\theta_{\rm
S}=\pi/3,\ \phi_{\rm S}=\pi/4,\ \theta_{\rm K}=2,\phi_{\rm
K}=5,\Phi_{\phi_{0}}=1.5$. The secondary spin is fixed to the fiducial value
$\chi=0.5$.

Details of our sampling algorithm can be found in App. A, including starting
points and prior bounds. In this analysis, we will use the emcee sampler due
to the simplicity of the likelihood structure in the case of circular orbits.

###  IV.1 Systematic biases — missing 1PA terms for quasicircular orbits

We first consider the impact of using mismodeled waveform templates on LISA
science when attempting to extract full 1PA waveforms within the data stream.

For each mass ratio $\epsilon=\\{10^{-5},10^{-4},10^{-3}\\}$ with respective
parameters given by Table 1, we perform four parameter estimation simulations.
We inject a true reference signal cir1PA with (w) spin and recover with the
following models each without (w/o) secondary spin

$h_{m}=\begin{cases}\text{cir1PA w/o spin},\\\ \text{cir0PA+1PA-3PN w/o
spin},\\\ \text{cir0PA w/o spin}.\end{cases}$ (43)

The results for each studied mass ratio
$\epsilon=\\{10^{-5},10^{-4},10^{-3}\\}$ are displayed (from top to bottom
respectively) in Fig. 1. In each of the three panels, the top rows (blue
posteriors) are exact marginalized posteriors, generated when injecting and
recovering with the exact model cir1PA w/ spin where the spin on the secondary
is sampled over. We do not present the posteriors on the extrinsic parameters
as they display near-to-zero bias with respect to the true parameters. The
non-Gaussian features and shifts to the true posterior are a feature of the
secondary spin, which will be discussed in Sec. IV.2. We begin by discussing
the case with the smallest mass ratio, $\epsilon=10^{-5}$.

Referring to the top panel of Fig. 1, we see that both the approximate models
cir1PA w/o spin (green) and cir0PA + 1PA-3PN w/o spin (red) are suitable for
parameter estimation of full 1PA waveforms. Each posterior shows statistically
insignificant biases at $\rho_{AET}\sim 70$ with Eq. (35) resulting in
$\mathcal{C}=1$ for all parameters. Remarkably, the parameters of the exact
cir1PA w/ spin model can be correctly inferred, with statistically
insignificant biases, by the cir0PA+1PA-3PN model. We remind the reader that
the latter contains 0PA information with a 1PA term approximated by a resummed
3PN expansion. The last row of the top panel in Fig. 1 shows a cir1PA w/ spin
model recovered with our least faithful model, cir0PA w/o spin. The intrinsic
parameters show statistically significant biases with $\mathcal{C}=0$, but the
recovered parameters are very similar to the true ones. For example, given the
true primary mass, $M=10^{6}M_{\odot}$, our best-fit parameter is only $\sim
10M_{\odot}$ away in magnitude. This quantitatively confirms that adiabatic
models would be fine for detection purposes, but not for statistical
inference.

Figure 1: Marginalized posteriors with shaded 68% credible intervals generated by injecting a true reference model cir1PA and recovering using different models with different mass ratios. The top panel is with $\epsilon=10^{-5}$, middle panel $\epsilon=10^{-4}$ and bottom panel $\epsilon=10^{-3}$. (Blue:) recovery with the true reference model cir1PA and sampling over the secondary spin parameter. (Green:) recovery using the cir1PA without spin. (Red:) recovery with cir0PA+1PA-3PN without spin. (Purple:) recovery with cir0PA (purely adiabatic) waveform. The black vertical dashed line indicates the true parameters. These simulations used parameters given by Table 1 for Configs. (1), (2) and (3) are given by the top row, middle row, and bottom row respectively. Orbital dephasings, mismatches, accumulated SNRs, and maximum log-likelihood values can be found in Table 2. $\epsilon$ | Model Waveform | $\Delta\Phi^{(\text{inj})}$ | $\Delta\Phi^{(\text{bf})}$ | $\mathcal{M}^{(\text{inj})}$ | $\mathcal{M}^{(\text{bf})}$ | $\rho^{(\text{inj})}/\rho^{(\text{opt})}$ | $\rho^{(\text{bf})}/\rho^{(\text{opt})}$ | $\log\mathcal{L}^{(\text{inj})}$ | $\log\mathcal{L}^{(\text{bf})}$
---|---|---|---|---|---|---|---|---|---
$10^{-5}$ |  |  |  |  |  |  |  |  |
Cir1PA w/o spin | $0.779$ | $0.0165$ | $0.143$ | $4.497\times 10^{-5}$ | $83.4\%$ | $99.9\%$ | $-846$ | -0.250
Cir0PA 1PA-3PN w/o spin | $0.786$ | $0.00179$ | $0.163$ | $4.293\times 10^{-6}$ | $81.5\%$ | $99.8\%$ | $-943$ | -0.0324
Cir0PA w/o spin | $3.002$ | $0.00532$ | $0.889$ | $2.412\times 10^{-6}$ | $6.4\%$ | $99.8\%$ | $-4800$ | -0.0234
$10^{-4}$ |  |  |  |  |  |  |  |  |
Cir1PA w/o spin | 3.994 | 0.00702 | 0.511 | $8.601\times 10^{-6}$ | $30.3\%$ | $99.9\%$ | -5019 | -0.336
Cir0PA 1PA-3PN w/o spin | 4.310 | 0.0179 | 0.486 | $1.26\times 10^{-4}$ | $34.2\%$ | $99.9\%$ | -4799 | -0.441
Cir0PA w/o spin | 13.093 | 0.0354 | 0.653 | $2.573\times 10^{-5}$ | $19.0\%$ | $99.9\%$ | -5506 | -0.122
$10^{-3}$ |  |  |  |  |  |  |  |  |
Cir1PA w/o spin | 4.518 | 0.00559 | 0.922 | $3.643\times 10^{-6}$ | $3.3\%$ | $99.9\%$ | -112938 | -0.226
Cir0PA 1PA-3PN w/o spin | 4.882 | 0.0218 | 0.949 | $3.443\times 10^{-5}$ | $3.4\%$ | $99.9\%$ | -112827 | -2.132
Cir0PA w/o spin | 14.958 | 0.153 | 0.938 | $6.854\times 10^{-3}$ | $4.9\%$ | $99.1\%$ | -122173 | -524.798
Table 2: Here we present a summary of computed statistics for various mass
ratios $\epsilon=\\{10^{-5},10^{-4},10^{-3}\\}$ (first column) when comparing
an injected cir1PA waveform and approximate model templates (second column).
We compute the orbital dephasings (Eqs. 39-40) (third and fourth columns);
mismatch (Eqs. 37-38) (fifth and sixth columns); accumulated SNRs (Eqs. 41-42)
(seventh and eighth columns); and, finally the log-likelihood function, Eq.
(34), at the injected/recovered parameters. The top, middle, and bottom panels
of this table correspond to the top, middle, and bottom panels of Fig. 1,
respectively.

We now discuss our simulations for $\epsilon\sim 10^{-4}$, given by the middle
panel in Fig. 1. The marginalized Gaussians in the top row of the middle panel
all exhibit heavy tails. This is due to correlations between the intrinsic
parameters and secondary spin, which will be discussed in Sec. IV.2. Since the
approximate models do not contain spin, there are no such correlations and the
resultant posteriors resemble their familiar Gaussian shapes for high SNRs.
With all approximate model templates, we see significant biases
($\mathcal{C}=0$) across the intrinsic parameters when neglecting post-
adiabatic terms at $\rho_{AET}\sim 65$. The most dramatic bias arises when
using the cir0PA w/o spin model to recover the exact cir1PA w/ spin reference
model. We conclude here that our approximate models summarized in Eq. (43),
without secondary spin, are unsuitable for parameter estimation.

Finally, we considered an IMRI with $\epsilon=10^{-3}$, given by the bottom
panel in Fig. 1. There are clear differences between the top row
(inject/recover with exact model) and the bottom three rows where we recover
the reference model with approximations (43) absent of spin. In the top row of
the bottom panel, the marginalized posteriors are not Gaussian and are non-
trivially skewed. This is a result of strong correlations between the
intrinsic parameters and the secondary spin. The Gaussian posteriors on the
second and third rows, representing model templates cir1PA w/o spin and cir0PA
+ 1PA-3PN w/o spin, show significant deviations from the true parameters with
respect to their own statistical uncertainty ($1\sigma$ width). From the
approximate model distributions, the uncertainties on the recovered parameters
are exceptionally small, reflected by the tight marginalized distributions.
This is a consequence of neglecting the important correlations with the
secondary spin. The recovered parameters are undoubtedly biased: the true
intrinsic parameters do not lie within their 68% credible interval.
Interestingly, the recovered parameters are contained within the 68% credible
interval of the marginalized true posterior distribution, implying that they
are consistent with the true parameter distribution. This is alarming: not
only are the incorrect parameters recovered, but our confidence that they are
the “correct” ones is largely inflated due to the tightness of the posteriors.
Finally, we see that the cir0PA w/o spin model features much stronger biases
and constraints than both the cir0PA + 1PA-3PN w/o spin and cir1PA w/o spin
models. In light of Eq. (35), we conclude that all models are unsuitable for
parameter inference of the cir1PA w/ spin model at $\rho_{AET}\sim 340$.

In Table 2, we give a summary of details regarding the individual MCMC
simulations for each small-mass-ratio binary configuration presented in Table
1. The details of the specific computations can be found in the caption of the
table. One of the main features of this table is the small mismatch, and
accumulated SNR normalized by the optimal SNR. In the worst case,
$\mathcal{M}\sim 10^{-3}$ and $\rho^{\text{(bf)}}/\rho^{\text{(opt)}}\sim
99.1\%$ for $\epsilon=10^{-3}$, between the injected cir1PA waveform and
adiabatic cir0PA w/o spin model evaluated at the recovered parameters. This
approximate template nearly matches the optimal matched filtering SNR, the SNR
that would be attained if the exact model was used during inference. This is
further evidence that, for quasicircular binaries, adiabatic models could be
used for detection purposes. Finally, we remark from Table 2 that all unbiased
results _satisfy_ the condition $\Delta\Phi^{\text{(inj)}}\lesssim 1$ radian.

###  IV.2 Constraining the secondary spin

We now focus our attention on constraining the spin $\chi$ of the CO. Similar
to Sec. IV.1, we study each mass ratio
$\epsilon=\\{10^{-5},10^{-4},10^{-3}\\}$ with parameters given by Table 1, and
perform three parameter estimation simulations. The injection is a cir1PA w/
spin model, and approximate waveforms are similar to (43) but with spin
included:

$h_{m}=\begin{cases}\text{cir0PA+1PA-3PN w/ spin},\\\ \text{cir0PA w/
spin}.\end{cases}$ (44)

Our corner plots for each of the $\epsilon=\\{10^{-5},10^{-4},10^{-3}\\}$ are
displayed in Figs. 2, 3 and 4, respectively. We will first discuss the
$\epsilon=10^{-5}$ case.

From Fig. 2, we see that the parameter $\chi$ _cannot_ be constrained for the
$\epsilon=10^{-5}$ case at $\rho_{AET}\sim 70$. The marginalized posterior
distribution for $\chi$ is almost flat. This implies that our posterior
information is not dominated by the likelihood (a function of the data), but
instead dominated by the prior (a function of the parameters, irrespective of
the data). We have tested various values of $\chi=\\{-1,0.5,0,0.5,1\\}$, and
in no situation can the secondary spin be constrained. The exact model cir1PA
w/ spin and approximate cir0PA + 1PA-3PN w/ spin model are indistinguishable.
When recovering the exact cir1PA w/ spin with the exact model itself and the
approximate cir0PA + 1PA-3PN w/ spin model, statistically insignificant biases
to the intrinsic parameters are observed. For the case with mass ratio
$\epsilon=10^{-5}$, the spin on the secondary can be compensated by the minor
tweaking of the intrinsic parameters. Finally, we report statistically
significant biases when employing the cir0PA w/ spin model to extract
parameters from a cir1PA waveform. These biases are consistent with the top
panel and fourth row of Fig. 1. To conclude, neither the exact model nor
approximate model templates are able to detect the presence of the spin on the
smaller companion for a mass ratio $\epsilon=10^{-5}$ and $\rho_{AET}\sim 70$.

Our results for mass ratio $\epsilon=10^{-4}$ with $\rho_{AET}\sim 65$ are
given in Fig. 3. The secondary spin has a more noticeable impact, and it can
be constrained when using model templates given by cir1PA w/ spin and cir0PA +
1PA-3PN w/ spin. From the 2D marginalized posteriors, it is evident that the
distributions on the intrinsic parameters and secondary spin are not Gaussian
due to the presence of heavy tails, highlighting strong correlations between
these parameters. Recall the red posterior in the middle panel of Fig. 1,
where biases on parameters were observed if we recovered a cir1PA w/ spin
template using a cir0PA + 1PA-3PN w/o spin model. We see from the red
posterior in Fig. 3 that including the spin on the cir0PA + 1PA-3PN model
eliminates the biases on parameters, making it completely indistinguishable
from the true cir1PA w/ spin model. This indicates that the cir0PA + 1PA-3PN
would be suitable for parameter estimation, _only if_ the secondary spin
parameter was included in the approximate model Piovano:2020ooe . Finally, we
see that the cir0PA w/ spin model fails to constrain the secondary spin, with
biases consistent with the second panel and fourth row of Fig. 1. Thus,
neglecting 1PA components of the GSF will have a detrimental effect on
recovering the spin of the smaller companion.

To conclude this section, we now discuss the impact of the secondary spin on
IMRIs with a mass ratio $\epsilon=10^{-3}$ with $\rho_{AET}\sim 340$. Our
results are displayed in Fig. 4 for each of the various model templates. In
contrast to the previous results for smaller mass-ratios, here we find the
only waveform suitable for parameter estimation is the exact model cir1PA w/
spin. The cir0PA + 1PA-3PN w/spin model yields biases on the intrinsic
parameters and the secondary spin, although it accounts for correlations
between the parameters. By contrast, the cir0PA w/ spin model exhibits
significantly stronger biases because it does not correctly represent such
correlations due to the lack of 1PA information. This can be seen in the 2D
marginalized posteriors. The “hard cut-offs” observed in the two posteriors
for cir1PA w/ spin and cir1PA + 1PA-3PN w/ spin are due to the spin on the
secondary reaching the prior bounds. These cut-offs are not physical but
merely a sampling artifact. We remark here that had we chosen a uniform prior
with more support, say $\chi\in[-2,2]$, then we could have made a wrong
conclusion on the nature of the spinning companion using the cir0PA + 1PA-3PN
w/ spin model. This analysis suggests that the bias on the secondary spin is
unacceptable for all approximate models and one must be careful, in this case,
of using PN results to approximate the 1PA components of the self-force. For
IMRIs, it is essential that we have full access to 1PA waveforms when
performing parameter estimation.

We conclude this section by briefly discussing the impact of first post-
adiabatic effects on small-mass-ratio binaries with $\epsilon<10^{-5}$.
Although the main sources considered in this work are small-mass-ratio
binaries with $\epsilon\in[10^{-5},10^{-3}]$, we have also studied a strong-
field EMRI with $\epsilon=10^{-6}$ at an SNR $\rho_{AET}\sim 23$. Our analysis
indicates that the secondary parameter cannot be constrained and that
parameter estimation studies can be conducted with 0PA waveforms. That is, the
1PA components of the self-force are negligible for such small mass-ratios. We
do not present our posterior results here, but the result can be extrapolated
from 0PA results of Fig. 1. The primary mass $M$ shows the strongest level of
bias, with value approximately $\sim\epsilon M$. For a mass ratio
$\epsilon=10^{-6}$, the bias on the primary mass $M$ is
$\mathcal{O}(1M_{\odot})$. This bias is well contained within the statistical
error given by the approximate posterior, making 0PA waveforms suitable for
parameter estimation. This implies that the true parameter is contained within
the 68% credible interval of the 0PA approximate posterior, satisfying Eq.
(35) with $\mathcal{C}=1$. We conclude that adiabatic waveforms are suitable
for _both_ search and characterizing 1PA waveforms at $\epsilon=10^{-6}$, at
least for quasicircular systems with nonspinning primaries.

In the analysis above we have tested three specific EMRI/IMRI configurations
at mass-ratios $\epsilon\in\\{10^{-6},10^{-5},10^{-4},10^{-3}\\}$ with SNR
$\rho_{AET}\in\\{22,70,65,340\\}$ respectively. We remind the reader that the
conclusions in the four paragraphs above are governed by two quantities: (1)
the specific configuration of parameters describing the system, the resultant
SNR and (2) the _geometry_ of the inspiral and resultant waveform. Our
Schwarzschild inspirals are strong-field orbits with trajectories evolved to
within $\hat{r}\sim 6.27$. However, the orbits could evolve much closer to the
horizon for a spinning primary. The increase in the number of orbits,
compounded with closer-horizon geometry, would increase significantly the
precision in parameters _and_ the SNR. This is strongly supported by work in
Babak:17aa , see Fig. 6 and Fig. 11. See also Ref. Burke:20aa  for a detailed
discussion on enhancements on parameter precision for circular orbits into
rotating MBHs. For more complex orbits, it is then expected that one could
achieve SNRs exceeding those of what we present here. As a result, the
accuracy requirements on EMRIs would become more stringent.

On the other hand, some might view our choices of luminosity distances as
overly optimistic. We assume a flat-$\Lambda$CDM cosmological model with
matter density $\Omega_{\rm m}=0.274$, dark energy density
$\Omega_{\Lambda}=0.726$ and Hubble constant $H_{0}=70.5\text{km
s}^{-1}\,\text{Mpc}$. Our choices for the luminosity distances $D_{\rm
S}\in\\{1,2\\}$Gpc correspond to small redshifts
$z_{\text{S}}=\\{0.203,0.371\\}$, less than $z_{\rm S}=1$. By comparison,
$z_{\rm S}=1$ corresponds to a luminosity distance of $D_{\rm S}(z_{\rm
S}=1)=6.716\,$Gpc. Fig. 9 in Ref. Babak:17aa  represents the redshift
distributions of detected EMRI events assuming 12 distinct astrophysical
models. All such distributions are peaked between $z_{\rm S}\in[1,2]$. Thus,
our choices represent _golden_ EMRIs: strong-field in their orbital
characteristics and placed at low redshifts $z_{S}\ll 1$. Assuming an
astrophysically relevant luminosity distance $D_{\rm S}=6.67$, the SNR of our
sources would decrease significantly since $\rho\sim 1/D_{L}$. In fact, our
$\epsilon=10^{-5}$ case in Table 1 would only reach an SNR $\rho_{AET}\sim
11$, which is lower than the detection threshold for EMRIs.

This highlights a fundamental point. Since the statistical error on the
parameters scales with the SNR, one can choose an SNR $\sim 20$ for
$\epsilon=10^{-5}$ such that 0PA waveforms are suitable for parameter
estimation of 1PA waveforms. For a more relevant, but still quite small,
luminosity distance $D_{S}\sim 3\,$Gpc, corresponding to $z_{\rm S}=0.52$,
$\rho_{AET}\sim 20$ for $\epsilon=10^{-5}$. The true parameter could then be
within the 68% credible interval of the purple posterior in Fig. 1. For such a
system we could then conclude that adiabatic waveforms are suitable for _both_
detection and parameter extraction of 1PA waveforms. Similarly, one could
place the source at exceptionally low redshift $z_{\rm S}=2.5\times 10^{-4}$,
giving a luminosity distance $D_{\rm S}=0.01\,$ Gpc and increasing the SNR by
a factor of 100 in the $\epsilon=10^{-5}$ case. The spinning secondary may be
observable in such a situation, but the probability that such an EMRI will be
observed is essentially zero according to Babak:17aa .

In conclusion, we have chosen parameter sets such that all of our orbits are
astrophysically sound, allowing us to draw reasonable conclusions. Within the
range of realistic sources, we have focused on ones that are loud enough to be
of most interest for high-precision science. The high-accuracy models built by
the self-force community naturally aim for these golden EMRIs (and IMRIs),
which have high SNRs and are in the strong-field regime.

Figure 2: (Mass ratio $\epsilon=10^{-5}$): Here we inject an EMRI waveform
with _spinning_ CO on a circular orbit with parameters $M=10^{6}M_{\odot}$,
$\mu=10M_{\odot}$, $r_{0}/M=10.6025$, $D_{S}=1\,$Gpc and extrinsic parameters
given in the caption of Table 1. The magnitude of the spinning secondary is
$\chi=0.5$, the SNR is $\rho_{AET}\sim 70$ and the time of observation is
$T_{\text{obs}}=2$ years. The blue, red, and purple parameter posteriors are
generated when recovering with a cir1PA w/ spin model, cir0PA + 1PA-3PN w/
spin model, and cir0PA w/ spin model, respectively. The black vertical lines
indicate the true parameters. The take-home message is that the spin of the
secondary _cannot_ be constrained for $\rho_{AET}=70$ and the models
considered here.

Figure 3: (Mass ratio $\epsilon=10^{-4}$): Here we inject an EMRI waveform
with _spinning_ CO on a circular orbit with parameters $M=10^{6}M_{\odot}$,
$\mu=100M_{\odot}$, $r_{0}/M=15.7905$, $D_{S}=2\,$Gpc and extrinsic parameters
given in the caption of Table 1. The magnitude of the secondary spin is
$\chi=0.5$, the SNR is $\rho_{AET}\sim 65$ and the time of observation is
$T_{\text{obs}}=1.5$ years. The blue, red, and purple parameter posteriors are
generated when recovering with a cir1PA w/ spin model, cir0PA + 1PA-3PN w/
spin model, and cir0PA w/ spin models respectively. The black vertical lines
indicate the true parameters. The take-home message here is that the spin of
the secondary _can_ be constrained using either the cir1PA w/ spin or cir1PA +
1PA-3PN w/ spin model. The cir0PA w/ spin model yields significant biases and
_cannot_ constrain the spin of the secondary.

Figure 4: (Mass ratio $\epsilon=10^{-3}$): Here we inject an EMRI waveform
with _spinning_ CO on a circular orbit with parameters $M=5\cdot
10^{6}M_{\odot}$, $\mu=5000M_{\odot}$, $r_{0}/M=16.81230$, $D_{S}=1\,$Gpc and
extrinsic parameters given in the caption of Table 1. The magnitude of the
secondary spin is $\chi=0.5$, the SNR is $\rho_{AET}\sim 340$ and the time of
observation is $T_{\text{obs}}=1$ year. The blue, red, and purple parameter
posteriors are generated when recovering with a cir1PA w/ spin model, cir0PA +
1PA-3PN w/ spin model, and cir0PA w/ spin model, respectively. The black
vertical lines indicate the true parameters. The take-home message here is
that all 1PA terms (including the spin on the secondary) must be included to
perform parameter estimation. The PN approximated waveform cir0PA + 1PA-3PN w/
spin shows significant biases and yields biased results on the secondary spin.
The cir0PA w/ spin model is unable to account for correlations and constrain
the spin on the secondary body.

###  IV.3 Systematic biases — mismodeling evolution of eccentric orbits

Our previous analyses focus on a very restricted class of binary
configurations, and it is not obvious how our results will extend to more
generic systems. In this section we begin to explore that question by
assessing the potential impact of mismodeling EMRI waveforms for eccentric
orbits. As 1PA models do not yet exist for eccentric orbits, we inject
waveforms with our adiabatic ecc0PA model (see Sec. II.2) and attempt to
recover them with our approximate ecc0PA-9PN model. We consider injected
waveforms with various values of eccentricity $e=\\{0,0.01,0.1,0.2\\}$.

The trajectories are evolved using a low-eccentricity 9PN expansion,
exhibiting slow convergence as the eccentricity increases. We do not consider
eccentricities $e>0.2$ for two reasons. The first is the potential for the PN
expansion to break down, yielding artificially inflated biases that could be
less severe in reality. The second reason is the enormous difficulty for the
sampler to converge to the point of highest likelihood. For $e=0.3$, we
noticed that the number of secondary maxima in the likelihood surface grew
significantly in comparison to $e\in\\{0.1,0.2\\}$. This means that the
individual chains of the sampler got stuck for $\sim\mathcal{O}(10^{3})$
iterations even with $\sim 10$ temperatures (see App. A for more details).
Whether this is an artifact of eccentricity itself, or the PN expansion
breaking down, is unclear. For the aforementioned reasons, we consider weak-
field orbits with low eccentricities to avoid exaggerating the impact of
eccentricity when mismodeling templates.

For all cases presented in this section, the injected waveforms’ parameters
are set to $M=10^{6}M_{\odot}$, $\mu=10M_{\odot}$, and initial semi-latus
rectum $p_{0}=9.86$. The trajectories are evolved for $T_{\text{obs}}=1$ year
until a final semi-latus rectum of $p\sim 8.1$ is reached. We do not evolve
the trajectories further than this point due to the breakdown of the PN
expansions in the strong-field regime. We choose the same extrinsic parameters
as presented in Table 1, but with luminosity distance $D_{L}=0.7\,$Gpc and
initial radial phase $\Phi_{r_{0}}=3$ for all cases. For each eccentricity
$e=\\{0,0.01,0.1,0.2\\}$, Eq. (39) gives a dephasing on the order
$\Delta\Phi^{\text{(inj)}}\approx\\{14,15,19,35\\}$ radians, respectively,
between the two models. Finally, the number of modes in the waveform for both
models, $\\#\text{modes}$, for each chosen eccentricity $e$ is
$(e,\\#\text{modes})=\\{(0,12),(0.01,12),(0.1,60),(0.2,94)\\}$ respectively.
Finally, both circular and eccentric waveform models yield similar SNRs on the
order of $\rho_{AET}\sim 70$. The details of how we used eryn for our
eccentric parameter estimation simulations are presented in App. A.

We begin with the circular orbit case with the result shown in Fig. 5. The
details on the individual runs can be found in the caption. We treat the
parameters $e=0$ and $\Phi_{r_{0}}=3$ as known, and we do not sample over
them. The intrinsic parameters exhibit statistically significant biases,
whereas the extrinsic parameters are unbiased. Notice that the “directions” of
the biases are similar to those presented when recovering (circular) post-
adiabatic waveforms with adiabatic templates in Fig. 1. Figure 5 will be our
reference figure when making direct comparisons with eccentric orbits.

We perform a series of similar simulations of recovering an exact adiabatic
model (ecc0PA), with an approximate (ecc0PA-9PN) model template with small to
moderate eccentricities $e=\\{0.01,0.1,0.2\\}$. The case with $e=0.01$ is
qualitatively similar to Fig. 5, so we will not present it here. The results
for $e=0.2$ are shown in Fig. 6. All the intrinsic parameters show severe
levels of biases, stronger compared to the circular orbits. Furthermore,
unlike the circular case, the angular parameters
$\\{\theta_{S},\phi_{S},\theta_{K},\phi_{K}\\}$ and initial phases
$\\{\Phi_{\phi_{0}},\Phi_{r_{0}}\\}$ show statistically significant levels of
bias. It should be noted that the magnitude of the biases on all parameters
increases significantly between the $e=0$ and $e=0.2$ cases presented in Fig.
5 and Fig. 6, respectively.

The biases on the extrinsic parameters for eccentric orbits stem from two
effects. The first is the correlations between the parameters, which are more
pronounced compared to circular orbits. Unlike the circular orbit case, the
intrinsic and extrinsic parameter spaces are not orthogonal: minor tweaks to
the intrinsic parameters can be compensated by minor tweaks in the extrinsic
parameters. The second reason is that the orbital evolution is more complex,
resulting in a LISA-responsed waveform with a richer structure in comparison
to its circular counterpart. A biased result on the intrinsic parameters
(mainly eccentricity) will induce modulations to the approximate model
template. Minor shifts to the angular parameters will induce further
modulations due to the presence of the LISA response function. The tweaking of
the angular parameters in response to the biases in the intrinsic parameters
appears to minimize the mismatch between the two signals. This can be explored
by comparing the ecc0PA waveform at the true parameters and ecc0PA-9PN
waveform at the recovered parameters, but fixing the angular parameters
$\\{\theta_{S},\phi_{S},\theta_{K},\phi_{K}\\}$ to their true values. We
obtain mismatches $\mathcal{M}\sim 0.283$ and accumulated SNR $\rho^{\text{(bf
w/ inj angles)}}/\rho^{(\text{opt})}\sim 68\%$, indicating that the two
waveforms quickly go out of phase. Comparing the true ecc0PA signal with
ecc0PA-9PN evaluated at the recovered parameters, including the recovered
angular parameters, yields $\mathcal{M}\sim 10^{-2}$ and accumulated SNR
$\rho^{\text{(bf)}}/\rho^{\text{(opt)}}\sim 98\%$, indicating that the
approximate model template remains in phase with the true reference signal for
a much longer duration.

For circular orbits, there are weak correlations between the intrinsic and
extrinsic parameters, demonstrated by the marginalized 2D distributions in
Fig. 5. Minor shifts to the intrinsic parameters will not affect the harmonic
structure, and thus there will be no biases across the extrinsic parameters.

Figure 5:  (Circular Orbits:) Here we consider as injection an ecc0PA model
(vanilla adiabatic FEW) with parameters $M=10^{6}M_{\odot}$,
$\mu=10M_{\odot}$, $r_{0}/M=9.6$, $e=0$. The trajectory evolves for one year
and terminates at a radial coordinate of $r_{0}/M\approx 8.1$, giving an SNR
$\rho_{AET}\sim 70$. We treat parameters related to eccentricity
$\\{e=0,\Phi_{r_{0}}=3\\}$ as known. The blue posterior is generated through
recovery using the injected model and the red posterior is built using the
approximate ecc0PA-9PN model. The black vertical lines indicate the location
of the true parameters.

Figure 6: (Eccentric orbits:) The same parameters and set up as Fig. 5 but
with $e=0.2$ and $\Phi_{r_{0}}=3$. Compared to Fig. 5, the biases on the
intrinsic parameters are more severe across the parameter space. We also
observe severe biases on the extrinsic parameters, a feature not seen for the
circular orbit case presented in Fig. 5.

##  V Summary

In this section, we summarize the results presented in Sec. IV. For the
configuration of parameters in this work (see Table 1), adiabatic templates
_are not suitable_ for EMRI parameter estimation, where the exact model
contains full post-adiabatic information. Adiabatic templates though, are
suitable for detection purposes in all cases. We also highlight the
performance of suitably re-summed third-order PN expansions when approximating
the 1PA components of the self-force. For mass ratios $\epsilon=10^{-5}$ and
$\epsilon=10^{-4}$, such approximate PN-based models are suitable for EMRI
data analysis (at least for the simple binary configurations we consider),
assuming that the spin on the smaller companion is included. We also found
that such PN-based approximate waveforms break down within the IMRI regime
$\epsilon=10^{-3}$, where full post-adiabatic information on the model
templates is required. Neglecting the spin on the secondary and/or the 1PA
components of the self-force results in both significant biases and
fictitiously tighter constraints on parameters. This is due to neglecting the
significant correlations between the secondary spin and intrinsic parameters.
We conclude that the spin of the secondary should always be retained and only
post-adiabatic templates should be used to characterize post-adiabatic models.

The feasibility of constraining the secondary spin was then discussed. For
$\epsilon=10^{-5}$, the secondary spin cannot be constrained for our choice of
SNR $\rho_{AET}\sim 70$, whereas, for $\epsilon=10^{-4}$ and
$\epsilon=10^{-3}$, it could be constrained at $\rho_{AET}\sim 65$ and
$\rho_{AET}\sim 340$, respectively. In the IMRI regime, the secondary spin and
intrinsic parameters exhibit strong correlations, which significantly degrade
the precision measurement of the intrinsic parameters.

Finally, we investigated the impact of mismodeling eccentric binaries. By
injecting an adiabatic eccentric waveform and recovering with a ninth-order
PN-based waveform, we observed severe biases across _both_ the intrinsic and
extrinsic parameters, whereas only the intrinsic parameters are biased for
purely circular orbits. It is impossible to say, for now, how this will
translate when using 0PA waveforms for inference on 1PA waveforms. However, we
expect that the severity of the biases with respect to circular orbits will
increase across the full parameter space. We conclude here that care must be
taken for both the modeling and parameter inference of eccentric 1PA
waveforms.

##  VI Discussion

This paper presented, for the first time, a detailed Bayesian study using MCMC
with state-of-the-art _first post-adiabatic_ EMRI and IMRI waveforms for
circular orbits in the Schwarzschild spacetime. We showed that neglecting 1PA
terms will induce statistically significant biases in the parameters. However,
we can still detect and characterize the first post-adiabatic waveform in the
data stream using adiabatic waveform models with biased parameters. We have
confirmed that adiabatic waveforms are only suited for detection purposes (at
least for $10^{-5}\lesssim\epsilon\lesssim 10^{-3}$). The systematic errors
are subjectively quite small for adiabatic waveforms and small mass-
ratios444For example, the top panel of Fig. 1 shows that we can recover the
primary mass $M=10^{6}M_{\odot}$ with a bias on the order of $\sim
10M_{\odot}$., and might not be relevant for some astrophysical applications
like population studies 10.1093/mnras/stad1397 ; Babak:2017tow . However,
applications within fundamental physics (like testing alternative/modified
theories of gravity Canizares:2012is ; Pani:2011xj ; Yunes:2011aa ; Blazquez-
Salcedo:2016enn , investigating the nature of massive black holes
Maggio:2021uge ; Barack:2006pq ; Barack:2006pq ; Datta:2019epe , the presence
of additional fields Fell:2023mtf ; Maselli:2021men ; Barsanti:2022vvl , and
so on) crucially rely on both precise and accurate measurements of the binary
parameters. Even relatively small statistically significant systematic errors
could spoil the enormous scientific potential of EMRIs. Thus, first post-
adiabatic waveforms are _essential_ in order to reap the full scientific
rewards of EMRI data analysis.

Bayesian methods are the gold standard technique when studying waveform
systematics because they do not introduce any approximations from a
statistical point of view. Fisher matrices, the Lindblom criterion,
overlaps/mismatches, and orbital dephasings must be used with caution. Such
systematic tests are suitable for exploration and gaining insight into the
accuracy of the waveform models, but conclusions must be taken lightly with
respect to their accurate Bayesian counterpart. Fisher matrices are hard to
accurately compute for EMRIs (see the works of Refs. Burke:20aa ;
Piovano:2020zin ; speri2021assessing ; Maselli:2021men ), potentially leading
to false conclusions on biases and precision statements of parameters.
Mismatches give no information about potential biases on parameters, and fully
optimized fitting factors such as Eq. (38) can only be calculated through
stochastic sampling algorithms. Similarly, the Lindblom criterion is overly
conservative purrer2020gravitational , and if taken at face value could force
unreasonably stringent accuracy requirements on waveform templates. Finally,
comparisons of the orbital dephasings between two models are performed at the
trajectory level, and so neglect the waveform structure, SNR of the source or
correlations between parameters. These systematic studies tell an important
story, but proper Bayesian inference completes the picture. For example, we
have shown that the posterior densities of EMRIs and IMRIs can yield non-
Gaussian features when the secondary spin is included, becoming ever more
dramatic as $\epsilon\gg 10^{-5}$. No other systematic test can reveal such
interesting features.

This importance of Bayesian methods can be seen in Sec. IV.1, which presents
the parameter distributions from mis-modelling in Fig. 1 and the associated
summary statistics in Table 2. We highlight from this analysis that
statistically significant biases _are not observed_ when the orbital dephasing
$\Delta\Phi^{\text{(inj)}}\lesssim 1$ radians (see Eq. (39)). From the top
panel and bottom row of Fig. 1, for $\epsilon=10^{-5}$, the orbital dephasing
between the 0PA and 1PA waveforms is $\sim 3$ radians. We see that strong
biases are observed, yet a concrete detection is made with accumulated SNR at
the best-fit parameters $\rho^{\text{(bf)}}/\rho^{\text{(opt)}}\sim 99.8\%$.
This is further exemplified in the bottom panel, where there is a $\sim 15$
radian difference between the 0PA and 1PA waveforms and
$\rho^{\text{(bf)}}/\rho^{\text{(opt)}}\sim 99.1\%$, resulting in a clear
detection. This is comforting to see, as it implies that the requirements
often used by the modelling community (e.g., that orbital phase errors should
be $\Delta\phi\lesssim 1/{\rm SNR}$ rad Hughes:2016xwf ) may be too stringent.
However, the orbital trajectories completely neglect the SNR of the system, a
_crucial_ ingredient for both detection and parameter estimation. Relying on
only orbital dephasing should be taken with caution. Similarly, large
mismatches between the approximate and exact waveform model (for identical
parameter values) might suggest that the waveform is not detectable with the
approximate model. But this is not the case: Table 2 contradicts it, giving
small mismatches $\mathcal{M}\sim 10^{-3}$ at the recovered parameters for the
largest of mass ratios $\epsilon=10^{-3}$ even though there are large
mismatches $\mathcal{M}\sim 0.938$ at the injected parameters. In conclusion,
our results reinforce that Bayesian inference is key to making significant
claims about approximations of waveforms.

Another outcome of our work is that we have shown the importance of including
the secondary spin parameter in our waveform models. Although it appears
impossible to constrain at mass ratios $\epsilon\leq 10^{-5}$, a measurement
can be made at higher mass ratios $\epsilon\geq 10^{-4}$. Correlations between
the intrinsic parameters and the secondary spin are significant, leading to
degraded measurements for IMRIs only if the spin parameter is included. We
have demonstrated that neglecting the secondary spin causes not only a bias
but fictitiously tight constraints on the other intrinsic parameters. This
becomes more prominent as the mass ratio increases. Furthermore, we have shown
that knowledge of the 1PA components of the GSF is essential when trying to
measure the secondary spin. It may not possible to make definite conclusions
on the potential constraints on the secondary spin for genetic orbits without
including other post-adiabatic terms.

Finally, we assessed the importance of eccentricity when mismodeling eccentric
templates. At the moment, only some self-force contributions at 1PA are known
for eccentric orbits, therefore it is not possible yet to make tests similar
to Sec. IV.1 and Sec. IV.2. Instead, we injected an adiabatic waveform and
attempted to recover with an approximate adiabatic waveform with trajectories
evolved through 9PN fluxes expanded in eccentricity. For low to moderate
eccentricities, $e\in\\{0,0.01,0.1,0.2\\}$, it is clear that the severity of
the biases worsens in comparison to the circular orbit case as $e$ increases
(cf. Figs. 5 and 6). For moderate eccentricities, biases are observed across
the entire parameter space, notably in the angular parameters and initial
phases. Such biases in the sky position would be unacceptable for studies
within cosmology, where the construction of galaxy catalogs allows one to
infer cosmological parameters, such as the Hubble constant
laghi2021gravitational . Our analysis show that the inclusion of eccentricity
will complicate the picture. More work is required on the self-force and data
analysis front to understand the impact of these orbits on EMRI parameter
estimation.

##  VII Future Work

The work presented here has only scratched the surface of EMRI accuracy
requirements. Clearly, it is essential to repeat this analysis once second-
order self-force results become available for more generic orbits. Indeed,
generic orbits may break potential degeneracies with the secondary spin
parameter, leading to improved constraints at smaller mass ratios
$\epsilon\sim 10^{-5}$. We have also shown remarkable success with the use of
resummed PN expansions when approximating the 1PA components in the context of
parameter estimation. For mass ratios $\epsilon\gtrsim 10^{-4}$, one could
perform preliminary studies on the secondary spin parameter, assuming suitable
1PA-$n$PN results were available for more general orbits. Moreover, we
observed that, due to correlations, the secondary spin parameter deteriorates
the precision with which the intrinsic parameters can be recovered. It would
then be interesting to understand whether there exist degeneracies between the
secondary spin and the scalar charge due to extra scalar fields, which may
spoil the precision measurements on the latter Barsanti:2022vvl ;
Maselli:2021men ; Barsanti:2022ana .

The field of EMRI systematics can now answer some crucial practical questions
for modeling EMRI waveforms. For instance, calculating the first-order self-
force for generic Kerr inspirals is very expensive with just a single point in
the parameter space taking $\mathcal{O}(10^{4})$ CPU hours vandeMeent:2017bcc
. One then has to repeat this calculation many times across the 4-dimensional
generic Kerr parameter space to produce interpolants for the 0PA equations of
motion, and then repeat the process at 1PA. It is important to estimate the
required accuracy at each point, the minimum number of points, and their
optimal placement, for both the 0PA and 1PA contributions in order to avoid
biasing parameters significantly. The choice of a suitable interpolation
method is important as well. It is still up for debate whether to use
Chebychev-interpolation methods (which have favorable convergence properties)
or less expensive splines.

We conclude by discussing one last topic that is still relatively untouched:
the EMRI search problem. The search problem has been “solved” in extremely
simplified circumstances by various groups within the LISA community
gair2008constrained ; cornish2011detection ; gair2008improved ; babak2010mock
. The underlying noise properties were well understood and tight priors were
placed on the parameters to recover. Furthermore, many of these groups
exploited the analytical features of the injected model (the self-inconsistent
“Analytical Kludge” waveforms from Ref. Barack:2003fp ), and the known
structure of the likelihood local maxima to reach the global maximum
indicating the true parameters. It is unknown whether fully relativistic
waveforms will simplify or complicate the search problem. EMRI search is a
difficult open problem, the authors hope that work can restart on this vital
research topic thanks to the advent of FEW and access to accurate waveforms.

###### Acknowledgements.

O. Burke thanks both the University College Dublin Relativity Group and the
Albert Einstein Institute for hosting him during the preparation of this work.
He also gives thanks to Alvin Chua, Christopher Chapman-Bird, Leor Barack,
Maarten Van de Meent, Sylvain Marsat, Michael Katz and Jonathan Gair for
various insightful discussions. This project has received financial support
from the CNRS through the MITI interdisciplinary programs. GAP acknowledges
support from an Irish Research Council Fellowship under grant number
GOIPD/2022/496. He also thanks Alvin Chua and Enrico Barausse for insightful
discussions. NW acknowledges support from a Royal Society - Science Foundation
Ireland University Research Fellowship. This publication has emanated from
research conducted with the financial support of Science Foundation Ireland
under Grant number 22/RS-URF-R/3825. AP acknowledges the support of a Royal
Society University Research Fellowship and a UKRI Frontier Research Grant
under the Horizon Europe Guarantee scheme [grant number EP/Y008251/1]. He also
thanks Alvin Chua and Leor Barack for helpful discussions. CK acknowledges
support from Science Foundation Ireland under Grant number 21/PATH-S/9610.
This work has used the following python packages:emcee, eryn, lisa-on-gpu,
scipy, numpy, matplotlib, corner, cupy, astropy, chainconsumer and lisatools
(Foreman-Mackey:2013, ; Karnesis:2023ras, ; michael_katz_2023_7705496, ;
harris2020array, ; Hunter:2007, ; corner, ; Hinton2016, ; astropy:2022, ;
cupy_learningsys2017, ; KatzTools, ). This work makes use of the following
packages of the Black Hole Perturbation toolkit BHPToolkit : FastEMRIWaveforms
michael_l_katz_2023_8190418 , KerrGeodesics KerrGeodesicsPackage  and
PostNewtonianSelfForce PostNewtonianSelfForcePackage .

## Author Contributions

OB: Conceptualization, Data curation, Formal analysis: MCMC simulations,
Investigation, Methodology, Project administration, Resources, Software: All
data analysis techniques, Supervision, Validation, Visualization: All plots,
Writing – original draft.

GP: Conceptualization, Data curation: secondary spin fluxes, Methodology,
Validation: Fisher based estimates, Investigation, Writing – original draft.

NW: Conceptualization, Data curation: second-order self-force results,
Methodology, Software: implementation of the waveform models, Writing - Review
& Editing.

PL: Conceptualization, Methodology, Software: implementation of the waveform
models, Writing - Review & Editing.

LS: Conceptualization, Methodology, Software: implementation of the waveform
models, Visualization, Writing - Review & Editing.

CK: Conceptualization, Methodology, Software: PN expansions, Writing - Review
& Editing.

BW: Conceptualization, Data curation: second-order self-force results,
Methodology, Writing - Review & Editing.

AP: Conceptualization, Data curation: second-order self-force results,
Methodology, Writing - Review & Editing.

LD: Data curation: second-order self force results

JM: Data curation: second-order self force results

## References

  * (1) L. Barack and C. Cutler, “Using lisa extreme-mass-ratio inspiral sources to test off-kerr deviations in the geometry of massive black holes,” _Physical Review D_ , vol. 75, no. 4, p. 042003, 2007.
  * (2) J. R. Gair, M. Vallisneri, S. L. Larson, and J. G. Baker, “Testing general relativity with low-frequency, space-based gravitational-wave detectors,” _Living Reviews in Relativity_ , vol. 16, pp. 1–109, 2013.
  * (3) A. Maselli, N. Franchini, L. Gualtieri, and T. P. Sotiriou, “Detecting scalar fields with Extreme Mass Ratio Inspirals,” _Phys. Rev. Lett._ , vol. 125, no. 14, p. 141101, 2020.
  * (4) A. Maselli, N. Franchini, L. Gualtieri, T. P. Sotiriou, S. Barsanti, and P. Pani, “Detecting fundamental fields with LISA observations of gravitational waves from extreme mass-ratio inspirals,” _Nature Astron._ , vol. 6, no. 4, pp. 464–470, 2022.
  * (5) S. Barsanti, N. Franchini, L. Gualtieri, A. Maselli, and T. P. Sotiriou, “Extreme mass-ratio inspirals as probes of scalar fields: eccentric equatorial orbits around Kerr black holes,” 3 2022.
  * (6) L. Barack and A. Pound, “Self-force and radiation reaction in general relativity,” _Reports on Progress in Physics_ , vol. 82, no. 1, p. 016904, 2018.
  * (7) A. Pound and B. Wardell, “Black hole perturbation theory and gravitational self-force,” _Handbook of Gravitational Wave Astronomy_ , pp. 1–119, 2022.
  * (8) L. Speri, M. L. Katz, A. J. Chua, S. A. Hughes, N. Warburton, J. E. Thompson, C. E. Chapman-Bird, and J. R. Gair, “Fast and fourier: Extreme mass ratio inspiral waveforms in the frequency domain,” _arXiv preprint arXiv:2307.12585_ , 2023.
  * (9) M. L. Katz, A. J. Chua, L. Speri, N. Warburton, and S. A. Hughes, “Fast extreme-mass-ratio-inspiral waveforms: New tools for millihertz gravitational-wave data analysis,” _Physical Review D_ , vol. 104, no. 6, p. 064047, 2021.
  * (10) A. J. Chua, M. L. Katz, N. Warburton, and S. A. Hughes, “Rapid generation of fully relativistic extreme-mass-ratio-inspiral waveform templates for lisa data analysis,” _Physical Review Letters_ , vol. 126, no. 5, p. 051102, 2021.
  * (11) M. Van De Meent and N. Warburton, “Fast Self-forced Inspirals,” _Class. Quant. Grav._ , vol. 35, no. 14, p. 144003, 2018.
  * (12) J. McCart, T. Osburn, and J. Y. J. Burton, “Highly eccentric extreme-mass-ratio-inspiral waveforms via fast self-forced inspirals,” _Phys. Rev. D_ , vol. 104, no. 8, p. 084050, 2021.
  * (13) Lynch, Philip and van de Meent, Maarten and Warburton, Niels, “Eccentric self-forced inspirals into a rotating black hole,” _Class. Quant. Grav._ , vol. 39, no. 14, p. 145004, 2022.
  * (14) P. Lynch, M. van de Meent, and N. Warburton, “Self-forced inspirals with spin-orbit precession,” _arXiv preprint arXiv:2305.10533_ , 2023.
  * (15) S. Isoyama, R. Fujita, A. J. K. Chua, H. Nakano, A. Pound, and N. Sago, “Adiabatic Waveforms from Extreme-Mass-Ratio Inspirals: An Analytical Approach,” _Phys. Rev. Lett._ , vol. 128, no. 23, p. 231101, 2022.
  * (16) P. Gupta, L. Speri, B. Bonga, A. J. K. Chua, and T. Tanaka, “Modeling transient resonances in extreme-mass-ratio inspirals,” _Phys. Rev. D_ , vol. 106, no. 10, p. 104001, 2022.
  * (17) L. Speri and J. R. Gair, “Assessing the impact of transient orbital resonances,” _Physical Review D_ , vol. 103, no. 12, p. 124032, 2021.
  * (18) P. Lynch, “Efficient trajectory calculations for extreme mass-ratio inspirals using near-identity (averaging) transformations,” Ph.D. dissertation, University College Dublin, 2022, available electronically at http://hdl.handle.net/10197/13347.
  * (19) S. Akcay, S. R. Dolan, C. Kavanagh, J. Moxon, N. Warburton, and B. Wardell, “Dissipation in extreme-mass ratio binaries with a spinning secondary,” _Phys. Rev. D_ , vol. 102, no. 6, p. 064013, 2020.
  * (20) G. A. Piovano, R. Brito, A. Maselli, and P. Pani, “Assessing the detectability of the secondary spin in extreme mass-ratio inspirals with fully relativistic numerical waveforms,” _Phys. Rev. D_ , vol. 104, no. 12, p. 124019, 2021.
  * (21) V. Skoupý and G. Lukes-Gerakopoulos, “Spinning test body orbiting around a Kerr black hole: Eccentric equatorial orbits and their asymptotic gravitational-wave fluxes,” _Phys. Rev. D_ , vol. 103, no. 10, p. 104045, 2021.
  * (22) J. Mathews, A. Pound, and B. Wardell, “Self-force calculations with a spinning secondary,” _Phys. Rev. D_ , vol. 105, no. 8, p. 084031, 2022.
  * (23) L. V. Drummond and S. A. Hughes, “Precisely computing bound orbits of spinning bodies around black holes. I. General framework and results for nearly equatorial orbits,” _Phys. Rev. D_ , vol. 105, no. 12, p. 124040, 2022.
  * (24) Drummond, Lisa V. and Hughes, Scott A., “Precisely computing bound orbits of spinning bodies around black holes. II. Generic orbits,” _Phys. Rev. D_ , vol. 105, no. 12, p. 124041, 2022.
  * (25) V. Skoupỳ, G. Lukes-Gerakopoulos, L. V. Drummond, and S. A. Hughes, “Asymptotic gravitational-wave fluxes from a spinning test body on generic orbits around a kerr black hole,” _arXiv preprint arXiv:2303.16798_ , 2023.
  * (26) S. A. Hughes, N. Warburton, G. Khanna, A. J. K. Chua, and M. L. Katz, “Adiabatic waveforms for extreme mass-ratio inspirals via multivoice decomposition in time and frequency,” _Phys. Rev. D_ , vol. 103, no. 10, p. 104014, 2021, [Erratum: Phys.Rev.D 107, 089901 (2023)].
  * (27) M. van de Meent, “Gravitational self-force on generic bound geodesics in Kerr spacetime,” _Phys. Rev. D_ , vol. 97, no. 10, p. 104033, 2018.
  * (28) B. Wardell, A. Pound, N. Warburton, J. Miller, L. Durkan, and A. Le Tiec, “Gravitational Waveforms for Compact Binaries from Second-Order Self-Force Theory,” _Phys. Rev. Lett._ , vol. 130, no. 24, p. 241402, 2023.
  * (29) S. Babak _et al._ , “The Mock LISA Data Challenges: From Challenge 3 to Challenge 4,” _Class. Quant. Grav._ , vol. 27, p. 084009, 2010.
  * (30) A. J. K. Chua, N. Korsakova, C. J. Moore, J. R. Gair, and S. Babak, “Gaussian processes for the interpolation and marginalization of waveform error in extreme-mass-ratio-inspiral parameter estimation,” _Phys. Rev. D_ , vol. 101, no. 4, p. 044027, 2020.
  * (31) A. J. Chua and C. J. Cutler, “Nonlocal parameter degeneracy in the intrinsic space of gravitational-wave signals from extreme-mass-ratio inspirals,” _Physical Review D_ , vol. 106, no. 12, p. 124046, 2022.
  * (32) J. R. Gair, L. Barack, T. Creighton, C. Cutler, S. L. Larson, E. S. Phinney, and M. Vallisneri, “Event rate estimates for lisa extreme mass ratio capture sources,” _Classical and Quantum Gravity_ , vol. 21, no. 20, p. S1595, 2004.
  * (33) S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, “Science with the space-based interferometer LISA. V. Extreme mass-ratio inspirals,” _prd_ , vol. 95, no. 10, p. 103012, May 2017.
  * (34) A. I. Burke and O. Burke, “Extreme precision and extreme complexity: source modelling and data analysis development for the laser interferometer space antenna,” Ph.D. dissertation, The University of Edinburgh, 2022, available electronically at https://era.ed.ac.uk/handle/1842/38372.
  * (35) T. Hinderer and E. E. Flanagan, “Two timescale analysis of extreme mass ratio inspirals in Kerr. I. Orbital Motion,” _Phys. Rev. D_ , vol. 78, p. 064028, 2008.
  * (36) R. Fujita, W. Hikida, and H. Tagoshi, “An Efficient Numerical Method for Computing Gravitational Waves Induced by a Particle Moving on Eccentric Inclined Orbits around a Kerr Black Hole,” _Prog. Theor. Phys._ , vol. 121, pp. 843–874, 2009.
  * (37) E. E. Flanagan, S. A. Hughes, and U. Ruangsri, “Resonantly enhanced and diminished strong-field gravitational-wave fluxes,” _Phys. Rev. D_ , vol. 89, no. 8, p. 084028, 2014.
  * (38) K. G. Arun _et al._ , “New horizons for fundamental physics with LISA,” _Living Rev. Rel._ , vol. 25, no. 1, p. 4, 2022.
  * (39) P. A. Seoane _et al._ , “Astrophysics with the Laser Interferometer Space Antenna,” _Living Rev. Rel._ , vol. 26, no. 1, p. 2, 2023.
  * (40) J. Miller and A. Pound, “Two-timescale evolution of extreme-mass-ratio inspirals: waveform generation scheme for quasicircular orbits in Schwarzschild spacetime,” _Phys. Rev. D_ , vol. 103, no. 6, p. 064048, 2021.
  * (41) A. Pound, B. Wardell, N. Warburton, and J. Miller, “Second-Order Self-Force Calculation of Gravitational Binding Energy in Compact Binaries,” _Phys. Rev. Lett._ , vol. 124, no. 2, p. 021101, 2020.
  * (42) N. Warburton, A. Pound, B. Wardell, J. Miller, and L. Durkan, “Gravitational-Wave Energy Flux for Compact Binaries through Second Order in the Mass Ratio,” _Phys. Rev. Lett._ , vol. 127, no. 15, p. 151102, 2021.
  * (43) A. Albertini, A. Nagar, A. Pound, N. Warburton, B. Wardell, L. Durkan, and J. Miller, “Comparing second-order gravitational self-force, numerical relativity, and effective one body waveforms from inspiralling, quasicircular, and nonspinning black hole binaries,” _Phys. Rev. D_ , vol. 106, no. 8, p. 084061, 2022.
  * (44) L. Barack and C. Cutler, “Lisa capture sources: Approximate waveforms, signal-to-noise ratios, and parameter estimation accuracy,” _Physical Review D_ , vol. 69, no. 8, p. 082005, 2004.
  * (45) S. Babak, H. Fang, J. R. Gair, K. Glampedakis, and S. A. Hughes, ““kludge” gravitational waveforms for a test-body orbiting a kerr black hole,” _Physical Review D_ , vol. 75, no. 2, p. 024005, 2007.
  * (46) A. J. Chua, C. J. Moore, and J. R. Gair, “Augmented kludge waveforms for detecting extreme-mass-ratio inspirals,” _Physical Review D_ , vol. 96, no. 4, p. 044005, 2017.
  * (47) L. Lindblom, B. J. Owen, and D. A. Brown, “Model waveform accuracy standards for gravitational wave data analysis,” _Physical Review D_ , vol. 78, no. 12, p. 124020, 2008.
  * (48) M. Vallisneri, “Use and abuse of the fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects,” _Physical Review D_ , vol. 77, no. 4, p. 042001, 2008.
  * (49) G. A. Piovano, A. Maselli, and P. Pani, “Constraining the tidal deformability of supermassive objects with extreme mass ratio inspirals and semianalytical frequency-domain waveforms,” _Phys. Rev. D_ , vol. 107, no. 2, p. 024021, 2023.
  * (50) P. Gupta, L. Speri, B. Bonga, A. J. Chua, and T. Tanaka, “Modeling transient resonances in extreme-mass-ratio inspirals,” _Physical Review D_ , vol. 106, no. 10, p. 104001, 2022.
  * (51) C. J. Moore, A. J. Chua, and J. R. Gair, “Gravitational waves from extreme mass ratio inspirals around bumpy black holes,” _Classical and Quantum Gravity_ , vol. 34, no. 19, p. 195009, 2017.
  * (52) A. Maselli, N. Franchini, L. Gualtieri, T. P. Sotiriou, S. Barsanti, and P. Pani, “Detecting fundamental fields with lisa observations of gravitational waves from extreme mass-ratio inspirals,” _Nature Astronomy_ , vol. 6, no. 4, pp. 464–470, 2022.
  * (53) C. Cutler and M. Vallisneri, “Lisa detections of massive black hole inspirals: Parameter extraction errors due to inaccurate template waveforms,” _Physical Review D_ , vol. 76, no. 10, p. 104018, 2007.
  * (54) T. Osburn, N. Warburton, and C. R. Evans, “Highly eccentric inspirals into a black hole,” _Phys. Rev. D_ , vol. 93, no. 6, p. 064024, 2016.
  * (55) N. Warburton, T. Osburn, and C. R. Evans, “Evolution of small-mass-ratio binaries with a spinning secondary,” _Phys. Rev. D_ , vol. 96, no. 8, p. 084057, 2017.
  * (56) L. V. Drummond, A. G. Hanselman, D. R. Becker, and S. A. Hughes, “Extreme mass-ratio inspiral of a spinning body into a kerr black hole i: Evolution along generic trajectories,” _arXiv preprint arXiv:2305.08919_ , 2023.
  * (57) M. L. Katz, J.-B. Bayle, A. J. Chua, and M. Vallisneri, “Assessing the data-analysis impact of lisa orbit approximations using a gpu-accelerated response model,” _Physical Review D_ , vol. 106, no. 10, p. 103001, 2022.
  * (58) W. Martens and E. Joffre, “Trajectory design for the esa lisa mission,” _The Journal of the Astronautical Sciences_ , vol. 68, no. 2, pp. 402–443, 2021.
  * (59) LISA Science Study Team, “Lisa science requirements document, esa-l3-est-sci-rs-001. technical report 1.0,” ESA, Tech. Rep., May 2018. [Online]. Available: https://www.cosmos.esa.int/web/lisa/lisa-documents/
  * (60) I. Mandel and J. R. Gair, “Can we detect intermediate mass ratio inspirals?” _Classical and Quantum Gravity_ , vol. 26, no. 9, p. 094036, 2009.
  * (61) I. Mandel, D. A. Brown, J. R. Gair, and M. C. Miller, “Rates and characteristics of intermediate mass ratio inspirals detectable by advanced ligo,” _The Astrophysical Journal_ , vol. 681, no. 2, p. 1431, 2008.
  * (62) M. van de Meent and H. P. Pfeiffer, “Intermediate mass-ratio black hole binaries: Applicability of small mass-ratio perturbation theory,” _Phys. Rev. Lett._ , vol. 125, no. 18, p. 181101, 2020.
  * (63) “Black hole perturbation toolkit,” (bhptoolkit.org), 2018.
  * (64) C. Munna, C. R. Evans, S. Hopper, and E. Forseth, “Determination of new coefficients in the angular momentum and energy fluxes at infinity to 9PN order for eccentric Schwarzschild extreme-mass-ratio inspirals using mode-by-mode fitting,” _Phys. Rev. D_ , vol. 102, no. 2, p. 024047, 2020.
  * (65) A. J. K. Chua, M. L. Katz, N. Warburton, and S. A. Hughes, “Rapid generation of fully relativistic extreme-mass-ratio-inspiral waveform templates for LISA data analysis,” _Phys. Rev. Lett._ , vol. 126, no. 5, p. 051102, 2021.
  * (66) M. L. Katz, A. J. K. Chua, L. Speri, N. Warburton, and S. A. Hughes, “Fast extreme-mass-ratio-inspiral waveforms: New tools for millihertz gravitational-wave data analysis,” _Phys. Rev. D_ , vol. 104, no. 6, p. 064047, 2021.
  * (67) A. I. Harte, “Motion in classical field theories and the foundations of the self-force problem,” _Fund. Theor. Phys._ , vol. 179, pp. 327–398, 2015.
  * (68) T. Tanaka, Y. Mino, M. Sasaki, and M. Shibata, “Gravitational waves from a spinning particle in circular orbits around a rotating black hole,” _Phys. Rev. D_ , vol. 54, pp. 3762–3777, 1996.
  * (69) W. Dixon, “A covariant multipole formalism for extended test bodies in general relativity,” _Il Nuovo Cimento_ , vol. 34, no. 2, pp. 317–339, Oct 1964.
  * (70) W. G. Dixon, “Dynamics of extended bodies in general relativity. I. Momentum and angular momentum,” _Proc. Roy. Soc. Lond._ , vol. A314, pp. 499–527, 1970.
  * (71) ——, “Dynamics of extended bodies in general relativity. II. Moments of the charge-current vector,” _Proc. Roy. Soc. Lond._ , vol. A319, pp. 509–547, 1970.
  * (72) G. A. Piovano, A. Maselli, and P. Pani, “Extreme mass ratio inspirals with spinning secondary: a detailed study of equatorial circular motion,” _Phys. Rev. D_ , vol. 102, no. 2, p. 024041, 2020.
  * (73) R. Rüdiger, “Conserved quantities of spinning test particles in general relativity. i,” _Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences_ , vol. 375, no. 1761, pp. 185–193, 1981. [Online]. Available: http://www.jstor.org/stable/2990231
  * (74) ——, “Conserved quantities of spinning test particles in general relativity. ii,” _Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences_ , vol. 385, no. 1788, pp. 229–239, 1983. [Online]. Available: http://www.jstor.org/stable/2397480
  * (75) V. Witzany and G. A. Piovano, “Analytic solutions for the motion of spinning particles near spherically symmetric black holes and exotic compact objects,” _arXiv preprint arXiv:2308.00021_ , 2023.
  * (76) P. I. Jefremov, O. Yu. Tsupko, and G. S. Bisnovatyi-Kogan, “Innermost stable circular orbits of spinning test particles in Schwarzschild and Kerr space-times,” _Phys. Rev._ , vol. D91, no. 12, p. 124030, 2015.
  * (77) A. Le Tiec, E. Barausse, and A. Buonanno, “Gravitational Self-Force Correction to the Binding Energy of Compact Binary Systems,” _Phys. Rev. Lett._ , vol. 108, p. 131103, 2012.
  * (78) K. Martel, “Gravitational wave forms from a point particle orbiting a Schwarzschild black hole,” _Phys. Rev. D_ , vol. 69, p. 044025, 2004.
  * (79) S. A. Hughes, “Bound orbits of a slowly evolving black hole,” _Phys. Rev. D_ , vol. 100, no. 6, p. 064001, 2019.
  * (80) L. Blanchet, G. Faye, Q. Henry, F. Larrouturou, and D. Trestini, “Gravitational-wave flux and quadrupole modes from quasicircular nonspinning compact binaries to the fourth post-Newtonian order,” _Phys. Rev. D_ , vol. 108, no. 6, p. 064041, 2023.
  * (81) R. Fujita, “Gravitational Waves from a Particle in Circular Orbits around a Schwarzschild Black Hole to the 22nd Post-Newtonian Order,” _Prog. Theor. Phys._ , vol. 128, pp. 971–992, 2012.
  * (82) A. Ori and K. S. Thorne, “Transition from inspiral to plunge for a compact body in a circular equatorial orbit around a massive, spinning black hole,” _Physical Review D_ , vol. 62, no. 12, p. 124022, 2000.
  * (83) R. O’Shaughnessy, “Transition from inspiral to plunge for eccentric equatorial kerr orbits,” _Physical Review D_ , vol. 67, no. 4, p. 044004, 2003.
  * (84) M. Kesden, “Transition from adiabatic inspiral to plunge into a spinning black hole,” _Physical Review D_ , vol. 83, no. 10, p. 104011, 2011.
  * (85) A. Apte and S. A. Hughes, “Exciting black hole modes via misaligned coalescences. i. inspiral, transition, and plunge trajectories using a generalized ori-thorne procedure,” _Physical Review D_ , vol. 100, no. 8, p. 084031, 2019.
  * (86) O. Burke, J. Gair, and J. Simón, “Transition from inspiral to plunge: A complete near-extremal trajectory and associated waveform,” _Physical Review D_ , vol. 101, no. 6, p. 064026, 2020.
  * (87) G. Compère, K. Fransen, and C. Jonas, “Transition from inspiral to plunge into a highly spinning black hole,” _Classical and Quantum Gravity_ , vol. 37, no. 9, p. 095013, 2020.
  * (88) G. Compère and L. Küchler, “Self-consistent adiabatic inspiral and transition motion,” _Physical review letters_ , vol. 126, no. 24, p. 241106, 2021.
  * (89) C. Cutler, D. Kennefick, and E. Poisson, “Gravitational radiation reaction for bound motion around a Schwarzschild black hole,” _Phys. Rev. D_ , vol. 50, pp. 3816–3835, 1994.
  * (90) R. Fujita and W. Hikida, “Analytical solutions of bound timelike geodesic orbits in Kerr spacetime,” _Class. Quant. Grav._ , vol. 26, p. 135002, 2009.
  * (91) N. Warburton, B. Wardell, O. Long, S. Upton, P. Lynch, Z. Nasipak, and L. C. Stein, “Kerrgeodesics,” Jul. 2023. [Online]. Available: https://doi.org/10.5281/zenodo.8108265
  * (92) D. Kennefick and A. Ori, “Radiation reaction induced evolution of circular orbits of particles around Kerr black holes,” _Phys. Rev. D_ , vol. 53, pp. 4319–4326, 1996.
  * (93) N. Warburton, B. Wardell, C. Munna, and C. Kavanagh, “Postnewtonianselfforce,” Jul. 2023. [Online]. Available: https://doi.org/10.5281/zenodo.8112975
  * (94) J. N. Goldberg, A. J. Macfarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan, “Spin-s Spherical Harmonics and ð,” _Journal of Mathematical Physics_ , vol. 8, no. 11, pp. 2155–2161, Nov. 1967.
  * (95) L. Barack and C. Cutler, “LISA capture sources: Approximate waveforms, signal-to-noise ratios, and parameter estimation accuracy,” _Phys. Rev. D_ , vol. 69, p. 082005, 2004.
  * (96) J. Armstrong, F. Estabrook, and M. Tinto, “Time-delay interferometry for space-based gravitational wave searches,” _The Astrophysical Journal_ , vol. 527, no. 2, p. 814, 1999.
  * (97) M. Tinto and S. V. Dhurandhar, “TIME DELAY,” _Living Rev. Rel._ , vol. 8, p. 4, 2005.
  * (98) M. Vallisneri, “Synthetic lisa: Simulating time delay interferometry in a model lisa,” _Physical Review D_ , vol. 71, no. 2, p. 022001, 2005.
  * (99) M. Vallisneri and C. R. Galley, “Non-sky-averaged sensitivity curves for space-based gravitational-wave observatories,” _Classical and Quantum Gravity_ , vol. 29, no. 12, p. 124015, 2012.
  * (100) C. Talbot, E. Thrane, S. Biscoveanu, and R. Smith, “Inference with finite time series: Observing the gravitational universe through windows,” _Physical Review Research_ , vol. 3, no. 4, p. 043049, 2021.
  * (101) O. Edy, A. Lundgren, and L. K. Nuttall, “Issues of mismodeling gravitational-wave data for parameter estimation,” _Physical Review D_ , vol. 103, no. 12, p. 124061, 2021.
  * (102) N. Wiener _et al._ , “Generalized harmonic analysis,” _Acta mathematica_ , vol. 55, pp. 117–258, 1930.
  * (103) A. Khintchine, “Korrelationstheorie der stationären stochastischen prozesse,” _Mathematische Annalen_ , vol. 109, no. 1, pp. 604–615, 1934.
  * (104) “LISA Science Study Team, “LISA Science Requirements Document,”,” 2018.
  * (105) Michael Katz, “Lisaanalysistools,” AEI, Tech. Rep., May 2020. [Online]. Available: https://github.com/mikekatz04/LISAanalysistools
  * (106) P. Whittle, “Curve and periodogram smoothing,” _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , vol. 19, pp. 38–63, 1957.
  * (107) L. S. Finn, “Detection, measurement, and gravitational radiation,” _Physical Review D_ , vol. 46, no. 12, p. 5236, 1992.
  * (108) P. M. Woodward, _Probability and information theory, with applications to radar: international series of monographs on electronics and instrumentation_. Elsevier, 2014, vol. 3.
  * (109) S. Babak, J. R. Gair, and E. K. Porter, “An algorithm for the detection of extreme mass ratio inspirals in lisa data,” _Classical and quantum gravity_ , vol. 26, no. 13, p. 135004, 2009.
  * (110) N. Karnesis, M. L. Katz, N. Korsakova, J. R. Gair, and N. Stergioulas, “Eryn: A multi-purpose sampler for bayesian inference,” _arXiv preprint arXiv:2303.02164_ , 2023.
  * (111) M. Katz, N. Karnesis, and N. Korsakova, “mikekatz04/eryn: first full release,” Mar. 2023, available at https://doi.org/10.5281/zenodo.7705496.
  * (112) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, “emcee: The MCMC Hammer,” _Astronomical Society of the Pacific_ , vol. 125, no. 925, p. 306, Mar. 2013.
  * (113) J. Goodman and J. Weare, “Ensemble samplers with affine invariance,” _Communications in applied mathematics and computational science_ , vol. 5, no. 1, pp. 65–80, 2010.
  * (114) N. Cornish and T. Robson, “Galactic binary science with the new lisa design,” _Journal of Physics: Conference Series_ , vol. 840, no. 1, p. 012024, may 2017. [Online]. Available: https://dx.doi.org/10.1088/1742-6596/840/1/012024
  * (115) G. A. Piovano, A. Maselli, and P. Pani, “Model independent tests of the Kerr bound with extreme mass ratio inspirals,” _Phys. Lett. B_ , vol. 811, p. 135860, 2020.
  * (116) O. Burke, J. R. Gair, J. Simón, and M. C. Edwards, “Constraining the spin parameter of near-extremal black holes using LISA,” _prd_ , vol. 102, no. 12, p. 124054, Dec. 2020.
  * (117) C. E. A. Chapman-Bird, C. P. L. Berry, and G. Woan, “Rapid determination of LISA sensitivity to extreme mass ratio inspirals with machine learning,” _Monthly Notices of the Royal Astronomical Society_ , vol. 522, no. 4, pp. 6043–6054, 05 2023. [Online]. Available: https://doi.org/10.1093/mnras/stad1397
  * (118) S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, “Science with the space-based interferometer LISA. V: Extreme mass-ratio inspirals,” _Phys. Rev. D_ , vol. 95, no. 10, p. 103012, 2017.
  * (119) P. Canizares, J. R. Gair, and C. F. Sopuerta, “Testing Chern-Simons Modified Gravity with Gravitational-Wave Detections of Extreme-Mass-Ratio Binaries,” _Phys. Rev. D_ , vol. 86, p. 044010, 2012.
  * (120) P. Pani, V. Cardoso, and L. Gualtieri, “Gravitational waves from extreme mass-ratio inspirals in Dynamical Chern-Simons gravity,” _Phys. Rev. D_ , vol. 83, p. 104048, 2011.
  * (121) N. Yunes, P. Pani, and V. Cardoso, “Gravitational Waves from Quasicircular Extreme Mass-Ratio Inspirals as Probes of Scalar-Tensor Theories,” _Phys. Rev. D_ , vol. 85, p. 102003, 2012.
  * (122) J. L. Blázquez-Salcedo, C. F. B. Macedo, V. Cardoso, V. Ferrari, L. Gualtieri, F. S. Khoo, J. Kunz, and P. Pani, “Perturbed black holes in Einstein-dilaton-Gauss-Bonnet gravity: Stability, ringdown, and gravitational-wave emission,” _Phys. Rev. D_ , vol. 94, no. 10, p. 104024, 2016.
  * (123) E. Maggio, M. van de Meent, and P. Pani, “Extreme mass-ratio inspirals around a spinning horizonless compact object,” _Phys. Rev. D_ , vol. 104, no. 10, p. 104026, 2021.
  * (124) L. Barack and C. Cutler, “Using LISA EMRI sources to test off-Kerr deviations in the geometry of massive black holes,” _Phys. Rev. D_ , vol. 75, p. 042003, 2007.
  * (125) S. Datta, R. Brito, S. Bose, P. Pani, and S. A. Hughes, “Tidal heating as a discriminator for horizons in extreme mass ratio inspirals,” _Phys. Rev. D_ , vol. 101, no. 4, p. 044004, 2020.
  * (126) S. D. B. Fell, L. Heisenberg, and D. Veske, “Detecting Fundamental Vector Fields with LISA,” 4 2023.
  * (127) S. Barsanti, A. Maselli, T. P. Sotiriou, and L. Gualtieri, “Detecting Massive Scalar Fields with Extreme Mass-Ratio Inspirals,” _Phys. Rev. Lett._ , vol. 131, no. 5, p. 051401, 2023.
  * (128) M. Pürrer and C.-J. Haster, “Gravitational waveform accuracy requirements for future ground-based detectors,” _Physical Review Research_ , vol. 2, no. 2, p. 023151, 2020.
  * (129) S. A. Hughes, “Adiabatic and post-adiabatic approaches to extreme mass ratio inspiral,” in _14th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories_ , vol. 2, 2017, pp. 1953–1959.
  * (130) D. Laghi, N. Tamanini, W. Del Pozzo, A. Sesana, J. Gair, S. Babak, and D. Izquierdo-Villalba, “Gravitational-wave cosmology with extreme mass-ratio inspirals,” _Monthly Notices of the Royal Astronomical Society_ , vol. 508, no. 3, pp. 4512–4531, 2021.
  * (131) J. R. Gair, E. Porter, S. Babak, and L. Barack, “A constrained metropolis–hastings search for emris in the mock lisa data challenge 1b,” _Classical and Quantum Gravity_ , vol. 25, no. 18, p. 184030, 2008.
  * (132) N. J. Cornish, “Detection strategies for extreme mass ratio inspirals,” _Classical and Quantum Gravity_ , vol. 28, no. 9, p. 094016, 2011.
  * (133) J. R. Gair, I. Mandel, and L. Wen, “Improved time–frequency analysis of extreme-mass-ratio inspiral signals in mock lisa data,” _Classical and Quantum Gravity_ , vol. 25, no. 18, p. 184031, 2008.
  * (134) S. Babak, J. G. Baker, M. J. Benacquista, N. J. Cornish, S. L. Larson, I. Mandel, S. T. McWilliams, A. Petiteau, E. K. Porter, E. L. Robinson _et al._ , “The mock lisa data challenges: from challenge 3 to challenge 4,” _Classical and Quantum Gravity_ , vol. 27, no. 8, p. 084009, 2010.
  * (135) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” _Nature_ , vol. 585, no. 7825, pp. 357–362, Sep. 2020. [Online]. Available: https://doi.org/10.1038/s41586-020-2649-2
  * (136) J. D. Hunter, “Matplotlib: A 2d graphics environment,” _Computing in Science & Engineering_, vol. 9, no. 3, pp. 90–95, 2007.
  * (137) D. Foreman-Mackey, “corner.py: Scatterplot matrices in python,” _The Journal of Open Source Software_ , vol. 1, no. 2, p. 24, jun 2016. [Online]. Available: https://doi.org/10.21105/joss.00024
  * (138) S. R. Hinton, “ChainConsumer,” _The Journal of Open Source Software_ , vol. 1, p. 00045, Aug. 2016.