of 16
Parameter control for eccentric, precessing binary black
hole simulations with SpEC
Taylor Knapp ,
1,2
,*
Katerina Chatziioannou ,
1,2
,
Harald Pfeiffer ,
3
,
Mark A. Scheel ,
1
and Lawrence E. Kidder
4
,
1
Department of Physics,
California Institute of Technology
, Pasadena, California 91125, USA
2
LIGO Laboratory,
California Institute of Technology
, Pasadena, California 91125, USA
3
Max Planck Institute for Gravitational Physics (Albert Einstein Institute)
,
Am Mühlenberg 1, 14476 Potsdam-Golm, Germany
4
Cornell Center for Astrophysics and Planetary Science,
Cornell University
, Ithaca, New York 14853, USA
(Received 16 October 2024; accepted 5 December 2024; published 2 January 2025)
Numerical relativity simulations of merging black holes provide the most accurate description of the
binary dynamics and the emitted gravitational wave signal. However, practical considerations such as
imperfect initial data and initial parameters mean that achieving target parameters, such as the orbital
eccentricity or the black hole spin directions, at the beginning of the usable part of the simulation is
challenging. In this paper, we devise a method to produce simulations with specific target parameters,
namely the Keplerian orbital parameters
eccentricity, semimajor axis, mean anomaly
and the black hole
spin vectors using SpEC. The method is an extension of the current process for achieving vanishing
eccentricity and it is based on a parameter control loop that iteratively numerically evolves the system, fits
the orbit with analytical post-Newtonian equations, and calculates updated input parameters. Through SpEC
numerical simulations, we demonstrate
10
3
and
O
ð
degree
Þ
convergence for the orbital eccentricity and
the spin directions respectively in
7
iterations. These tests extend to binaries with mass ratios
q
3
,
eccentricities
e
0
.
65
, and spin magnitudes
j
χ
j
0
.
75
. Our method for controlling the orbital and spin
parameters of numerical simulations can be used to produce targeted simulations in sparsely covered regions
of the parameter space or study the dynamics of relativistic binaries.
DOI:
10.1103/PhysRevD.111.024003
I. INTRODUCTION
With no known analytic solutions, the two-body problem
in general relativity can only be solved exactly with full
numerical relativity (NR) simulations
[1]
. Such simulations
have numerous practical applications in calibrating wave-
form models
[2,3]
or serving as the basis of surrogate
models
[4]
that are used in gravitational wave data analysis,
e.g.,
[5
7]
. Moreover, they provide solutions to the full
spacetime and elucidate the properties of black holes (BHs)
and general relativistic dynamics. Modern NR codes such as
the spectral Einstein code (SpEC) have produced thousands
of simulations of coalescing BHs
[8]
.
NR solves an initial value problem where initial con-
ditions are evolved forward in time. For binary BH (BBH)
systems, one can freely specify the ratio of the BH masses,
the initial BH spin angular momenta (magnitude and
direction), and the initial coordinate positions and
velocities of the BHs. For puncture data
[9]
these are
indeed directly the free parameters. For quasiequilibrium
conformal thin sandwich data
[10
12]
which is the type
of initial-data considered here, as it provides initial lapse
and shift for evolutions and allows for nearly extremal
black holes
[13]
the positions and velocities of the BHs
are instead encoded by the binary
s orbital velocity, the BH
separation, and the BH relative radial velocity. Evolving
the initial data forward in time yields the full spacetime
dynamics, BH properties, and emitted gravitational wave
signal as a function of time.
However, in practice the initial data for a BBH simulation
do not correspond to a snapshot of a binary that has been
inspiraling since the infinite past. This is because the initial
data do not include the correct gravitational radiation that
was emitted in the past during the infinitely long inspiral,
nor do they include the appropriate tidal distortions of the
BHs. As a result, when an NR simulation begins to evolve,
initial transients occur before the BBH relaxes into a
quasiequilibrium state. The gravitational radiation emitted
during this process is known as
junk radiation.
Junk
radiation, if it is resolved by the simulation, is a physical
solution of the Einstein equations, but is not the solution of
interest. There have been multiple attempts to reduce the
*
Contact author: tknapp@caltech.edu
Contact author: kchatziioannou@caltech.edu
Contact author: harald.pfeiffer@aei.mpg.de
§
Contact author: scheel@tapir.caltech.edu
Contact author: kidder@astro.cornell.edu
PHYSICAL REVIEW D
111,
024003 (2025)
2470-0010
=
2025
=
111(2)
=
024003(16)
024003-1
© 2025 American Physical Society
amount of junk radiation in NR simulations
[14
22]
,but
currently the standard practice is to simply discard the first
few orbits until the junk radiation has decayed away
[8,23]
.
A further complication arises in selecting initial param-
eters. In Newtonian gravitation, it is straightforward to
choose the initial positions and velocities of two point
particles so that the resulting orbit has a desired semimajor
axis
a
, eccentricity
e
, and mean anomaly
l
at some
reference time. But in general relativity with no known
closed-form or ordinary differential equation solution,
achieving a BBH orbit with desired parameters is not
straightforward.
Among orbital parameters, the eccentricity is particu-
larly astrophysically relevant. Since gravitational radiation
reaction efficiently circularizes BBH orbits
[24]
, most
simulations have targeted vanishing eccentricities. Moti-
vated by this, the SpEC workflow includes an initial
eccentricity removal
stage that tunes the input param-
eters, i.e., the angular velocity, separation, and relative
radial velocity, to reduce eccentricity below a predeter-
mined threshold
[25,26]
. An iterative scheme makes initial
guesses for the input parameters, performs an NR simu-
lation for a few orbits, measures the eccentricity, and then
chooses updated input parameters for the next iteration.
However, recent hints of nonzero eccentricity in select
observed signals
[27
30]
and an emphasis of eccentricity
as a tool to study the BBH formation history
[31]
have
reinvigorated interest in eccentric BBHs
[32
34]
.
Another astrophysically relevant property is the BH
spin, which determines the length and morphology of the
signal
[35]
and whose value also carries information about
the BBH formation history. The NR spins are set as part of
the initial data construction, but the spin directions precess
if the spins are misaligned with the orbital angular
momentum
[36]
. If the first few orbits of the simulation
are discarded because of junk radiation, then the useable
part of the simulation will begin with slightly different spin
directions, whose values are difficult to predict because
they depend on how long the junk radiation lasts. Thus, it
is difficult to precisely control spin directions at the
beginning of the usable part, i.e., after junk radiation, of
a precessing NR simulation.
In this paper we extend the eccentricity removal pro-
cedure of Ref.
[25]
to nonzero eccentricity; our method
also specifies the spin directions at some chosen time
different than
t
¼
0
. After introducing the problem setup,
in Sec.
II
we describe the method. Given the target
parameters of the simulation we would like to evolve,
i.e., target Keplerian parameters and spins, we iteratively
evolve the system, extract the binary parameters at the
usable part of the data, and adjust the input parameters.
Each evolution proceeds for 3
5 orbits, at the end of which
we fit the numerical data with a post-Newtonian model and
extract the Keplerian parameters. We then use standard
root-finding techniques to compute new input parameters
that when evolved will result in a simulation that is closer
to the target parameters.
We validate this method in Secs.
III
and
IV
for systems
with aligned and precessing spins respectively. Selecting a
threshold of achieving the target eccentricity to within
7
×
10
4
motivated by the zero-eccentricity limit, we show
that
7
iterations are sufficient for binaries with mass ratio
q
3
, eccentricity
e
0
.
65
, and spin magnitude
χ
0
.
75
.
Simultaneously, the spin directions are fixed with
O
ð
degree
Þ
accuracy relative to their target directions. We
conclude in Sec.
V
.
II. INITIAL CONDITIONS FOR NUMERICAL
RELATIVITY
Consider a binary of two BHs with masses
m
A
,with
A
f
1
;
2
g
, total mass
M
m
1
þ
m
2
,massratio
q
m
1
=m
2
1
, and symmetric mass ratio
η
¼
q=
ð
1
þ
q
Þ
2
.
The initial dimensionless spins of the BHs are denoted
χ
A
.
To simulate such a binary, SpEC requires initial data
parameters
θ
ID
¼f
Ω
0
;v
r;
0
;D
0
g
where
Ω
0
is the orbital
angular frequency,
v
r;
0
is the coordinate relative radial
velocity, and
D
0
is the coordinate BH separation at
simulation time
t
¼
0
. Other SpEC papers adopt the
convention
̇
a
0
¼
v
r;
0
=D
0
instead of
v
r;
0
,e.g.,
[25,37]
.
We use
v
r;
0
here and reserve the letter
a
for the semimajor
axis of elliptical orbits.
As the evolution proceeds, the coordinate trajectories are
described by the BH positions
x
A
ð
t
Þ
. If the orbit were
Newtonian, the trajectory could be equivalently described
by a Keplerian parametrization in terms of orbital elements
θ
orb
¼f
a; e;
l
g
, where
a
is the semimajor axis,
e
is the
eccentricity, and
l
is the mean anomaly at some reference
time. The actual trajectories are however not Newtonian,
but rather post-Newtonian (PN) with an increasing number
of relevant corrections as the BHs approach merger
[38]
.As
we consider a small portion of the whole trajectory far from
merger, namely the first few orbital periods, and we still
approximately parametrize the orbit by
θ
orb
. We describe
this parametrization and correspondence between
θ
orb
and
θ
ID
in Sec.
II D
.
The goal is to devise a method of choosing appropriate
input values of the initial data parameters
θ
ID
;T
and the
initial spins
χ
A;T
ð
t
¼
0
Þ
so that we obtain simulations with
target orbital parameters
θ
orb
;T
and spins
χ
A;T
ð
t
¼
t
T
Þ
at
some time
t
T
. Throughout the paper, we use the subscript
T
to denote the
target
parameters, including both the orbital
parameters
θ
orb
;T
and the corresponding initial data param-
eters
θ
ID
;T
that when evolved will give an orbit with
θ
orb
;T
(and similarly for the spins). Figure
1
presents the problem
setup in schematic form. Input parameters
θ
ID
and initial
spins
χ
A
ð
t
¼
0
Þ
are set at the initial time
t
¼
0
and the
TAYLOR KNAPP
et al.
PHYS. REV. D
111,
024003 (2025)
024003-2
system is numerically evolved (blue). The usable portion of
the numerical data starts after the end of the junk radiation
phase (gray band), denoted at
t
T
(black dotted vertical line).
At that time, the spins
̃
χ
A
ð
t
¼
t
T
Þ
are read off the numerical
data, while
θ
orb
is obtained as described in Sec.
II D
.
When specifying the target spins
χ
A;T
ð
t
¼
t
T
Þ
at the
target time
t
T
, a suitable reference frame needs to be
selected. By convention, SpEC simulations begin at
t
¼
0
with a frame whose
z
-axis is aligned with the orbital
angular momentum and the BHs are on the
x
-axis, with the
larger BH on the positive
x
-axis. During the simulation, the
BHs move in these
inertial coordinates
, and the simulation
outputs the BH
inertial frame components
of the spins. At
time
t
T
at which we want to set the BH spins
the BHs
will in general not be on the inertial coordinate x-axis, and
will in general not move in the
xy
-plane of the inertial
coordinates. In order to specify target spins independent of
the precise inertial frame position of the black holes, we
introduce the
BH frame
, which is coorbiting in the sense
that its
x
-axis is always pointing from BH 1 to BH 2, and its
z
-axis is orthogonal to the instantaneous orbital plane. (At
t
¼
0
, the BH frame coincides with the inertial frame). For
clarity, we will denote spins in the BH frame with a vector,
e.g., the target spin is
χ
A;T
ð
t
¼
t
T
Þ
, and spins in the inertial
frame (which is the frame the NR data express them in)
with a tilde.
Two challenges arise when attempting to map from some
chosen initial parameters
θ
ID
and
χ
A
ð
t
¼
0
Þ
to orbital
parameters
θ
orb
and
χ
A
at time
t
¼
t
T
:
(1) Due to junk radiation, we have to set
t
T
100
M.
During that time, spin-precession changes the BH
spin directions in the inertial simulation frame. This
change happens on the (slow) precession timescale.
In the co-orbiting BH frame where we define the
target spins, the spins evolve on the much faster
orbital timescale. To minimize impact of this fast
timescale, we select
t
T
to correspond to an integer
number of orbits, after which the BH frame has
roughly returned to its initial position. Here, we
choose
t
T
to be one orbit. This choice is explained
further in Sec.
II E
.
(2) The mapping between
θ
ID
and
θ
orb
is not straightfor-
ward under full general relativistic dynamics, so that
initial data parameters
θ
ID
;T
that yield the desired
target orbital parameters
θ
orb
;T
are initially not
known. Therefore, we construct an iterative update
approach to obtain
θ
orb
;
T
through small updates to
θ
ID
. The iterative scheme is formalized and explained
in depth in the rest of this section.
For reference in the remainder of this text, we supply
Table
I
with definitions of parameter notations.
A. The inverse problem and root-finding
The map between the input parameters
θ
ID
at
t
¼
0
and
the orbital parameters
θ
orb
at a later time
t
T
is written as an
operator,
S
which includes
(1) The
forward
mapping from
θ
ID
at
t
¼
0
to the
numerical data at later time
t
T
; this is the SpEC
evolution.
FIG. 1. Schematic representation of the challenges in determin-
ing the initial conditions of numerical relativity simulations. We
use an example simulation and plot the derivative of the orbital
angular frequency as a function of time,
̇
Ω
ð
t
Þ
, in blue. The gray
region from
t
¼
0
to
t
¼
t
T
(dotted vertical line) is the window of
junk radiation. We overlay a post-Newtonian fitted solution
(dashed) using Eq.
(9)
. Our root-finding method is represented
by
G
and
S
as introduced in Sec.
II A
. The junk radiation
contributes a significant nonlinear effect in the first
100
Mof
the simulation.
TABLE I. A reference table for the initial parameters, orbital
parameters, and spin utilized in this work. An arrow or tilde over
each spin parameter refers to the BH or inertial frame, respect-
fully, that it is defined in.
Variable
Description
χ
A;T
ð
t
¼
t
T
Þ
Target spin vector for BH A at
t
¼
t
T
in the coorbting
BH frame
χ
A;T
ð
t
¼
0
Þ
Initial spin vector for BH A that when evolved will
result in
χ
A;T
ð
t
¼
t
T
Þ
χ
ð
i
Þ
A
ð
t
¼
0
Þ
Initial spin vector for BH A for the
i
th trial simulation
̃
χ
ð
i
Þ
A
ð
t
¼
t
T
Þ
Spin vector for BH A at
t
¼
t
T
in the inertial frame for
the
i
th trial simulation, read off the NR data
θ
orb
;T
Target orbital parameters at
t
¼
t
T
θ
ID
;T
Target initial data parameters that when evolved will
result in
θ
orb
;T
θ
ð
i
Þ
ID
Initial data parameters for the
i
th trial simulation
θ
ð
i
Þ
orb
Orbital parameters for the
i
th trial simulation
PARAMETER CONTROL FOR ECCENTRIC, PRECESSING
...
PHYS. REV. D
111,
024003 (2025)
024003-3
(2) The extraction of
θ
orb
that parametrizes the numeri-
cal data.
Collectively, these steps are denoted as,
S
ð
θ
ID
Þ
in Fig.
1
.
The inverse mapping is denoted by
G
in Fig.
1
. While
S
ð
θ
ID
Þ
can be fully calculated with SpEC,
G
is unknown,
which we will elaborate upon later.
For some target orbital parameters,
θ
orb
;T
, we seek
θ
ID
;
T
such that
S
ð
θ
ID
;
T
Þ¼
θ
orb
;T
:
ð
1
Þ
Suppose we begin with an initial data guess
θ
ð
i
Þ
ID
that yields
orbital parameters
θ
ð
i
Þ
orb
. The superscript here denotes the
i
th
trial evolution as part of an iterative scheme. If we update
θ
ð
i
Þ
ID
to
θ
ð
i
þ
1
Þ
ID
¼
θ
ð
i
Þ
ID
þ
δ
θ
ID
, a Taylor expansion yields
S
ð
θ
ð
i
Þ
ID
þ
δ
θ
ID
Þ¼
S
ð
θ
ð
i
Þ
ID
Þþ
δ
S
δ
θ
ID
δ
θ
ID
:
ð
2
Þ
Requiring that the updated evolution brings us to the target
parameters
S
ð
θ
ð
i
Þ
ID
þ
δ
θ
ID
Þ¼
θ
orb
;T
yields
δ
θ
ID
¼

δ
S
δ
θ
ID

1
h
θ
orb
;T
S
ð
θ
ð
i
Þ
ID
Þ
i
¼
δ
G
δ
θ
orb
h
θ
orb
;T
S
ð
θ
ð
i
Þ
ID
Þ
i
G
ð
θ
orb
;T
Þ
G
h
S
ð
θ
ð
i
Þ
ID
Þ
i
G
ð
θ
orb
;T
Þ
θ
ð
i
Þ
ID
;
ð
3
Þ
where
G
S
1
is the inverse operator. To go from the first
to the second line we have used the property of inverse
functions

δ
S
δ
θ
ID

1
¼
δ
G
δ
θ
orb
;
ð
4
Þ
while for the third line, we use the definition of the
directional derivative, assuming that
δ
θ
ID
is small.
Equation
(3)
provides a way to compute
δ
θ
ID
such that
θ
ð
i
þ
1
Þ
ID
¼
θ
ð
i
Þ
ID
þ
δ
θ
ID
brings us closer to the target orbital
parameters
θ
orb
;T
.
If we knew the operator
G
exactly, then the update would
converge in one iteration using the procedure above.
Unfortunately, the exact
G
is unknown as it corresponds
to the inverse operation of a full numerical simulation.
Therefore, the primary purpose of this paper is to derive
suitable approximations to
G
such that the iterative pro-
cedure converges in a reasonable number of iterations. In
Sec.
II B
we describe the iterative process, while in Sec.
II C
we use post-Newtonian Keplerian equations to approximate
the map between
θ
orb
and
θ
ID
.
B. Parameter control loop
As will be further discussed in Sec.
II C
, computing
G
includes a number of approximations. We therefore intro-
duce a parameter control loop that iteratively performs trial
simulations, extracts their parameters at
t
T
, updates the
initial parameters, and starts new simulations until achiev-
ing a predefined tolerance threshold.
The full workflow is presented in Fig.
2
.
(1) We begin with a user-selected set of target orbital
parameters
θ
orb
;T
and spins
χ
A;T
(box 1).
(2) We calculate the SpEC input parameters for the
i
¼
0
iteration
θ
ð
i
Þ
ID
based on Sec.
II F
and set
χ
ð
i
Þ
A
ð
t
¼
0
Þ¼
χ
A;T
. Unless otherwise indicated, by
default we start simulations at apastron (box 2).
(3) We numerically evolve the system up to a time
5
π
=
Ω
0
. In practice this corresponds to roughly 3 to 5
orbits, depending on the eccentricity of the orbit
(boxes 3
4).
(4) We fit the numerical data for the orbital parameters
θ
ð
i
Þ
orb
and read off the spins
̃
χ
ð
i
Þ
A
ð
t
¼
t
T
Þ
(box 5).
(5) If the orbital parameters
θ
ð
i
Þ
orb
are not equal to
the target parameters
θ
orb
;
T
within some tolerance
(box 6), we update the input parameters with
FIG. 2. Flow chart for the SpEC eccentricity and spin control
loop, an extension of the zero-eccentricity loop of Ref.
[25]
.
TAYLOR KNAPP
et al.
PHYS. REV. D
111,
024003 (2025)
024003-4
Eqs.
(11)
(13)
and the spins with Eq.
(14)
(box 7)
and start a new simulation (box 3).
(6) Once the eccentricity threshold is achieved (box 6),
we evolve the system to merger (box 8).
For simplicity, we adopt the eccentricity threshold
imposed by the initial implementation of the parameter
control loop in SpEC aiming at vanishing eccentricity
[25]
:
j
e
T
e
ð
i
Þ
j
7
×
10
4
. This tolerance is set just below the
eccentricity value for which spin-induced oscillations
dominate the eccentricity-induced oscillations in
̇
Ω
.
Therefore, it is a sufficient constraint for a parameter
control loop that relies on
̇
Ω
.
Even though we only set a tolerance threshold on the
eccentricity, in practice we find that all orbital and spin
parameters converge to their target value reasonably well.
Adding further tolerance thresholds on spins or further
orbital parameters is straightforward. Having established
the iterative parameter update procedure, we turn to the
post-Newtonian equations required for the fits and updates
of boxes 4 and 6.
C. Eccentric post-Newtonian dynamics
We fit the numerical data with the Keplerian equations of
eccentric motion augmented with post-Newtonian correc-
tions as given in Ref.
[39]
. We use these post-Newtonian
equations to extract the Keplerian parameters
θ
orb
from
numerical data and to compute the inverse operator
G
in Eq.
(3)
.
A generic Keplerian orbit at Newtonian order, 0PN, is
characterized by its semimajor axis
a
, eccentricity
e
, and
mean anomaly
l
defined at the epoch
t
¼
0
. These three
quantities are related via Kepler
s equations:
̄
Ω
t
¼
u
ð
t
Þ
e
sin
u
ð
t
Þ
l
and
̄
Ω
¼
ffiffiffiffiffiffiffiffiffiffiffiffi
M=a
3
p
. The eccentric anomaly,
u
ð
t
Þ
, is obtained by inverting Kepler
s equation, which in
the small-eccentricity limit can be approximated as:
u
ð
t
Þ¼
̄
Ω
t
þ
l
. To next, 1PN order, the Keplerian param-
eters and all equations obtain corrections. These can be
expressed in closed form and as a function of the
Newtonian orbital parameters
ð
a; e;
l
Þ
. We present these
equations and their derivation in Appendix. The pertinent
point is that there is a closed-form expressions for the
(derivative of the) orbital angular frequency as a function of
time,
̇
Ω
ð
t
;
a; e;
l
Þ
.
D. PN trajectory fitting model
We characterize the numerical trajectory through the
derivative of the orbital angular velocity
̇
Ω
ð
t
Þ
as shown in
Fig.
1
. This choice follows Ref.
[25]
, which states that
eccentricity-induced oscillations become much more pro-
nounced, and therefore easier to fit, in
̇
Ω
ð
t
Þ
than in
Ω
ð
t
Þ
.
In the post-Newtonian framework,
̇
Ω
ð
t
Þ
deviates from
zero (circular orbit) due to three effects
̇
Ω
ð
t
Þ¼
δ
̇
Ω
RR
ð
t
Þþ
δ
̇
Ω
SS
ð
t
Þþ
δ
̇
Ω
e
ð
t
Þ
:
ð
5
Þ
The first term in Eq.
(5)
encodes radiation reaction effects
which, including 1PN corrections, read
δ
̇
Ω
RR
ð
t
Þ¼
96
5
η
M
7
=
2

D
4
0
256
5
η
M
3
t

11
=
8
48
5
η
ð
3
η
Þ
M
9
=
2

D
4
0
256
5
η
M
3
t

13
=
8
;
ð
6
Þ
where
D
0
is the initial BH separation.
The second term in Eq.
(5)
encodes 2PN spin-spin
interactions
[25]
δ
̇
Ω
SS
ð
t
Þ¼
1
2
½ð
S
0
·
ˆ
n
0
Þ
2
þð
S
0
·
ˆ
λ
0
Þ
2

sin
ð
2
̄
Ω
t
þ
φ
Þ
;
ð
7
Þ
which induce frequency modulations at twice the orbital
frequency even for circular orbits. We use the definitions of
Ref.
[25]
, where
S
0
is the mass-weighted spin vector,
ˆ
n
0
is
the unit normal to the orbital plane, and
ˆ
λ
0
is orthogonal to
ˆ
n
0
and the angular orbital momentum vector,
ˆ
L
. Subscripts
indicate that these unit vectors are defined at
t
¼
0
. The
phase offset
φ
is a known function of
S
0
given by Eq. (50)
in Ref.
[25]
.
The final term in Eq.
(5)
encodes eccentric effects,
therefore it is a function of the orbital parameters. To 1PN
order
[39]
δ
̇
Ω
e
ð
t
Þ¼
A
sin
u
ð
t
Þð
1
̃
e
cos
u
ð
t
ÞÞ
ð
1
e
t
cos
u
ð
t
ÞÞ
3
ð
1
e
φ
cos
u
ð
t
ÞÞ
2
;
ð
8
Þ
where
A; e
t
;e
φ
;
̃
e; u
ð
t
Þ
are known functions of the
Keplerian orbital parameters
ð
a; e;
l
Þ
given in Appendix.
Appendix
A2
also presents a derivation of Eq.
(8)
based
on Ref.
[39]
. Crucially,
δ
̇
Ω
e
ð
t
Þ
at this order depends only on
ð
a; e;
l
Þ
, so we can use this expression to fit numerical data
and obtain estimates for the orbital parameters.
Motivated by the three terms of Eq.
(5)
, we define the
following model to fit the numerical data, shown by the
dashed curve in Fig.
1
̇
Ω
model
ð
t
Þ¼
C
1
ð
T
c
t
Þ
11
=
8
þ
C
2
ð
T
c
t
Þ
13
=
8
þ
C
3
cos
α
ð
t
Þþ
C
4
sin
α
ð
t
Þ
A
sin
u
ð
t
Þð
1
̃
e
cos
u
ð
t
ÞÞ
ð
1
e
t
cos
u
ð
t
ÞÞ
3
ð
1
e
φ
cos
u
ð
t
ÞÞ
2
:
ð
9
Þ
The first line of Eq.
(9)
captures radiation-reaction. It
depends on amplitude parameters
C
1
and
C
2
that we fit for.
The coalescence time,
T
c
, is calculated from the evolution
data using using a 0PN approximation for the time to
merger
PARAMETER CONTROL FOR ECCENTRIC, PRECESSING
...
PHYS. REV. D
111,
024003 (2025)
024003-5
T
c
¼
t
T
þ
5
M
2
256
η
ð
M
Ω
ð
t
T
ÞÞ
8
=
3
;
ð
10
Þ
where
Ω
ð
t
T
Þ
is the angular frequency at time
t
¼
t
T
, read
off from the numerical data.
The second line captures spin-spin interaction effects of
the functional form of Eq.
(7)
. It again depends on
amplitude parameters
C
3
and
C
4
that we fit for, while
α
ð
t
Þ
is a known function that is directly calculated from the
numerical data, see Eq. (50) in Ref.
[25]
. Finally, the third
line of Eq.
(9)
captures eccentricity-induced oscillations,
the primary objective of this work. All quantities appearing
in this term are known functions of the orbital parameters
a
,
e
,
l
that we fit for.
Formally, Eq.
(9)
should contain additional PN correc-
tions to each of the three effects beyond the PN order in
which they each appear. In practice, we find that this
expression is adequate for fitting numerical data with the
goal of obtaining updated input parameters.
E. Fitting the numerical data
We fit the numerical data with Eq.
(9)
to extract
the orbital parameters,
θ
orb
¼ð
a; e;
l
Þ
. This is a seven-
parameter nonlinear fit, but the values of four of the
parameters
ð
C
1
;C
2
;C
3
;C
4
Þ
are not used in any of the
subsequent formulas. First, we must choose a reference
time,
t
T
, the beginning of the fitting interval. In Fig.
1
,this
interval is indicated by the white region to the right of the
dotted line (
t
¼
t
T
). Our motivation for choosing
t
T
is
twofold. First, we only want to fit the trajectory after junk
radiation has ended. Second, we require
t
T
to be a multiple
of the orbital period such that the BHs have returned as
closely as possible to their starting positions. Both require-
ments are met with
t
T
set to 1 orbital period, in practice
occurring on the order of
300
M into the simulation
evolution. Currently, we cap
t
T
at 500 M but could update
the determination of the end of junk radiation using
Ref.
[23]
. We implement a cap in the event
5
π
=
Ω
0
yields
an inefficiently large
t
T
.
Next, we must decide how long to evolve the BBH after
t
T
. We motivate this choice with three arguments. First, the
fitting window must be short to mitigate computational cost
for the trial simulation whose main goal is to evaluate
Eq.
(3)
and provide improved initial parameters.
Second, the eccentricity, semimajor axis, and mean
anomaly are PN parameters that evolve due to radiation
reaction and become ill defined close to merger. The fitting
window must therefore both be short and as far away from
merger as possible. Third, since we are fitting harmonic
functions, at least a few orbits are required to achieve
enough precision. In practice, we evolve the trial simu-
lation for a little less than 5 orbits up to a final time
T
end
¼
10
π
=
Ω
0
, calculated based on the initial orbital
angular frequency. The PN fit then occurs between times
t
T
and
T
end
.
F. Updating the input parameters
At this stage, we have performed a trial simulation with
input parameters
θ
ð
i
Þ
ID
and spins
χ
ð
i
Þ
A
ð
t
¼
0
Þ
. Using Eq.
(9)
,
we have fitted the evolved trajectory to obtain orbital
parameters
θ
ð
i
Þ
orb
. Using the simulation data, we also read off
spins
̃
χ
ð
i
Þ
A
ð
t
¼
t
T
Þ
. The next step is to use this information to
compute updated initial-data parameters for the next
iteration. We introduce two independent update procedures
for the orbit and spins.
1. Orbit
Given the orbital parameters of the current iteration
θ
ð
i
Þ
orb
,
we calculate initial data parameters for the next iteration
θ
ð
i
þ
1
Þ
ID
. First, we need a mapping from
θ
orb
to
θ
ID
.
We evaluate
Ω
ð
θ
orb
;t
Þ
,
v
r
ð
θ
orb
;t
Þ
,and
D
ð
θ
orb
;t
Þ
using
Eqs.
(A18)
(A23)
, which are derived from 1PN equations
of motion in Appendix
A3
. Per Eq.
(3)
, the difference
between these functions evaluated for
θ
ð
i
Þ
orb
and
θ
orb
;
T
at
time
t
¼
0
is the amount we need to correct the current
iteration initial data parameters by, yielding
Ω
ð
i
þ
1
Þ
0
¼
Ω
ð
i
Þ
0
þ½
Ω
ð
θ
orb
;
T
;
0
Þ
Ω
ð
θ
ð
i
Þ
orb
;
0
Þ
;
ð
11
Þ
v
ð
i
þ
1
Þ
r;
0
¼
v
ð
i
Þ
r;
0
þ½
v
r
ð
θ
orb
;
T
;
0
Þ
v
r
ð
θ
ð
i
Þ
orb
;
0
Þ
;
ð
12
Þ
FIG. 3. Schematic for spin vector updates described in
Sec.
II F 2
. The bottom diagram shows the BBH at
t
¼
0
and
the top diagram shows the BBH at
t
¼
t
T
, one orbit later. The
vectors in the bottom panel represent the initial spin vector for the
current iteration,
χ
ð
i
Þ
ð
t
¼
0
Þ
, and the updated initial spin vector
for the next iteration,
χ
ð
i
þ
1
Þ
ð
t
¼
0
Þ
. The top panel shows the
target spin vector,
χ
T
ð
t
¼
t
T
Þ
, and the evolved spin vector at time
t
¼
t
T
, denoted
χ
ð
i
Þ
ð
t
¼
t
T
Þ
. The green and red arrows map
vectors between times
t
¼
0
and
t
T
via rotation matrix
R
0
t
T
and
its inverse, respectively. The spin arrows and text at times
t
¼
0
and
t
¼
T
correspond to the relevant rotation matrix and arrow.
Though the figure depicts only spin on one BH for clarity, we
apply the same procedure to both BHs
spins.
TAYLOR KNAPP
et al.
PHYS. REV. D
111,
024003 (2025)
024003-6
D
ð
i
þ
1
Þ
0
¼
D
ð
i
Þ
0
þ½
D
ð
θ
orb
;
T
;
0
Þ
D
ð
θ
ð
i
Þ
orb
;
0
Þ
:
ð
13
Þ
2. Spin
The spin update procedure is applied to both BHs
spins;
for clarity we discuss only one case and drop the BH
subscripts. We recall that the target spins
χ
T
are given in the
coorbiting BH frame at
t
¼
t
T
, whereas the numerical
evolution yields spin vectors in the inertial frame. Thus, the
first step is to rotate the target spins into the inertial frame.
We define this rotation matrix as
R
BH
in
ð
t
¼
t
T
Þ
, con-
structed using the inertial frame basis of
ˆ
n;
ˆ
λ
, and
ˆ
L
introduced in the context of Eq.
(7)
but defined at time
t
T
. Next, we correct for the change in the spin directions
during the evolution from
t
¼
0
to
t
¼
t
T
. We introduce
another rotation matrix, denoted
R
0
t
T
, shown in Fig.
3
,as
the rotation from
χ
ð
t
¼
0
Þ
to
χ
ð
t
¼
t
T
Þ
. We then apply the
inverse rotation
R
1
0
t
T
to
χ
T
to approximate the position of
that vector at the beginning of the simulation,
t
¼
0
.
Combining the two rotations yields the initial spin of the
next iteration
̃
χ
ð
i
þ
1
Þ
ð
t
¼
0
Þ¼
R
1
0
t
T
R
BH
in
ð
t
¼
t
T
Þ
χ
T
ð
t
¼
t
T
Þ
:
ð
14
Þ
The derivation and explicit forms of rotation matrices
R
1
0
t
T
and
R
BH
in
ð
t
¼
t
T
Þ
are given in Appendix
A4
.
III. ECCENTRIC, NONPRECESSING ORBITS
We begin with an example of the parameter control
process. Figure
4
shows results for an equal-mass binary,
target orbital parameters
θ
orb
;T
¼f
a
T
¼
15
M
;e
T
¼
0
.
2
;
l
T
¼
π
g
, and vanishing spins.
1
The top left panel shows
the numerical data for
̇
Ω
ð
t
Þ
for each of the three iterations
required to achieve the eccentricity threshold. Each iter-
ation follows the control process shown in Fig.
2
to produce
fitted orbital parameters for the numerical data and update
the initial data parameters of the next iteration accordingly.
The remaining panels show the difference between the
achieved (fitted from the NR data) and the target eccentricity
(top right), semimajor axis (bottom left), and mean anomaly
(bottom right) as a function of iteration. Here and in
subsequent figures,
a
ð
i
Þ
;e
ð
i
Þ
and
l
ð
i
Þ
refer to the fitted
orbital parameters in each iteration
i
of the control loop.
Even though we only impose a tolerance threshold on the
eccentricity, all orbital parameters converge to their target
values after 3 iterations. In fact, the mean anomaly,
l
,
converges almost immediately. Since this behavior is typical
FIG. 4. Application of the parameter control loop outlined in Fig.
2
for an equal-mass binary with vanishing spins and target orbital
parameters
θ
orb
;T
¼f
a
T
¼
15
M
;e
T
¼
0
.
2
;
l
T
¼
π
g
. The top left panel shows the NR data for each iteration of the control loop (top
subpanel) and the residuals between the PN fit and NR data for each iteration (bottom subpanel). The rightmost gray vertical line
indicates
t
¼
t
T
. The PN fitting window occurs between the two vertical gray lines. The difference between the achieved, i.e., fitted from
the NR data, orbital parameters in each iteration
i
and the target parameters is plotted in the remaining panels: eccentricity (top right),
semimajor axis (bottom left), and mean anomaly (bottom right). All parameters converge to their target value after 3 iterations.
1
Unless otherwise indicated, we set
l
T
to apastron as this
choice leads to more stable evolutions.
PARAMETER CONTROL FOR ECCENTRIC, PRECESSING
...
PHYS. REV. D
111,
024003 (2025)
024003-7
in subsequent explorations, we omit the mean anomaly from
now on. Moreover, the semimajor axis converges in 1 fewer
iteration than the eccentricity. This justifies continuing to
consider only the eccentricity in the parameter control
tolerance.
Convergence in all parameters is monotonic, as in every
iteration they are closer in absolute value to their target
value. However, while the eccentricity is larger than the
target value in the first iteration, it is slightly smaller in the
second iteration. This suggests that the updated parameters
from the first to the second iteration
overshoot
slightly.
This characteristic
oscillatory
convergence is typical for
the eccentricity parameter.
A. Impact of target eccentricity value
Having demonstrated the parameter control loop in
practice for one binary, we now explore the impact of
the target eccentricity value on convergence. We again use
equal masses, vanishing spins, and vary
e
T
between 0 and
0.65. In highly eccentric binaries, the periastron passage
might bring the BHs close enough that the post-Newtonian
approximation becomes inaccurate or the system directly
merges. The BH separation at periastron depends on
whether we increase the eccentricity at a fixed semimajor
axis
a
or apastron separation
r
a
¼
a
ð
1
þ
e
Þ
. For binaries
with
e
T
<
0
.
5
, we find that the former suffices and fix
a
T
¼
15
M. For binaries with
e
T
0
.
5
, we instead increase
the eccentricity at a fixed
r
a;T
¼
60
M.
Figure
5
shows parameter convergence at fixed
a
T
and
e
T
<
0
.
5
(top) and fixed
r
a;T
and
e
T
0
.
5
(bottom). We
show the difference between the target and fitted eccen-
tricity (left panels) and semimajor axis (right panels) as a
function of iteration and colored by the target eccentricity
value. The inlaid plots show the BH coordinate separation
over time, using the same color scheme. In the top row, the
separation reflects a constant
a
T
for each
e
T
<
0
.
5
.Asthe
eccentricity increases (darker purple), the average separa-
tion stays constant but the apastron and periastron passages
become more extreme. In the bottom row, all trajectories
begin at the same apastron and result in roughly similar
periastron values.
In all cases, i.e., for target eccentricty up to
e
T
¼
0
.
65
,
the eccentricity and semimajor axis converge to their target
value, with higher eccentricities requiring on average more
iterations. This is likely because higher eccentricity orbits
are more relativistic, thus the post-Newtonian approxima-
tion to the dynamics is less accurate. Such less accurate
fitting equations lead to less efficient initial data parameter
FIG. 5. Impact of target eccentricity value on parameter convergence. We plot the difference between the fitted and target eccentricity
(left panel) and semimajor axis (right panel) at each iteration for binaries with equal masses, vanishing spins,
a
T
¼
15
M (top) or
r
a;T
¼
60
M (bottom), and different values of target eccentricity
e
T
(color bar). We achieve parameter convergence for all the cases we
attempted, including eccentricities up to 0.65. For
e
T
>
0
.
65
, the BHs risk head-on collision, which terminates the control loop
prematurely before obtaining
θ
ID;T
. Thus we omit cases
e
T
>
0
.
65
in this work. As
e
T
increases to 0.65, more iterations are required for
parameter convergence. We also inlay separation over time plots based on the final (most converged) iteration of the control loop.
TAYLOR KNAPP
et al.
PHYS. REV. D
111,
024003 (2025)
024003-8
updates in each control loop iteration, thus requiring more
iterations to converge.
B. Aligned spins
Next, we explore the impact of aligned spins. Aligned
spins remain constant in both direction and magnitude
(modulo horizon absorption effects
[40]
). Therefore the
rotations of Fig.
3
reduce to the identity, and the spins are
not updated during the eccentricity-control iterations.
However, the spins still affect the binary dynamics and
rate of orbital decay and could therefore impact the orbital
parameter convergence. We study binaries with equal
masses,
a
T
¼
15
M,
e
T
¼f
0
.
1
;
0
.
2
;
0
.
3
g
, and vary the
amount of aligned spin
χ
z;A;T
¼
χ
z;B;T
in
½
0
.
7
;
0
.
7

.
Results are shown in Fig.
6
for the different target
eccentricities (top to bottom) and aligned spin values
(colorbar). In all cases, the eccentricity and semimajor
axis converge to their target value after
7
iterations.
The initial (iteration 1) evolution consistently under
(over)-shoots the target eccentricity for (anti)aligned spins.
Similar to the effect for larger target eccentricities, larger
absolute value aligned spin magnitudes require more
iterations for convergence.
IV. ECCENTRIC, PRECESSING ORBITS
In this section, we generalize to precessing binaries.
Going forward, we assume one spinning BH and no
aligned component for simplicity. In-plane spins precess
in the inertial frame, causing the BH frame to change
orientation. The parameter control loop of Fig.
2
and the
updating formulas in Eqs.
(11)
(13)
and
(14)
treat the
orbital parameters and spins independently for conven-
ience. Since the orbital parameters and the spins are
updated independently, we begin by exploring two param-
eter control methods:
(1)
Simultaneous.
We follow the process of Fig.
2
where
in the first iteration we simulate an eccentric,
precessing binary and then iteratively update the
orbital parameters and the spins, per Eqs.
(11)
(13)
and
(14)
.
(2)
Sequential.
We execute the parameter control loop
twice. First, we iterate while updating only the
orbital parameters; while doing so, we keep the
spins constant at only the aligned-spin component of
the target spins. Subsequently, once the target orbital
parameters have been achieved, we restart the loop
FIG. 6. Impact of aligned spins on parameter control. We plot the difference between the fitted and target eccentricity (left panel) and
semimajor axis (right panel) at each iteration for binaries with equal masses,
a
T
¼
15
M,
e
T
¼f
0
.
1
;
0
.
2
;
0
.
3
g
(top to bottom) and
different values of the aligned spin (color bar). We successfully converge to the target parameters for eccentric, aligned spin binaries,
with higher eccentricities again requiring more iterations.
PARAMETER CONTROL FOR ECCENTRIC, PRECESSING
...
PHYS. REV. D
111,
024003 (2025)
024003-9