Gravitational waveforms of binary neutron star inspirals using PN Tidal Splicing
Kevin Barkett,
1
Yanbei Chen,
1
Mark A. Scheel,
1
and Vijay Varma
1
1
TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA
(Dated: November 26, 2019)
The tidal deformations of neutron stars within an inspiraling compact binary alter the orbital dy-
namics, imprinting a signature on the gravitational wave signal. Modeling this signal could be done
with numerical-relativity simulations, but these are too computationally expensive for many appli-
cations. Analytic post-Newtonian treatments are limited by unknown higher-order non-tidal terms.
This paper further builds upon the “Tidal Splicing” model in which post-Newtonian tidal terms are
“spliced” onto numerical relativity simulations of black-hole binaries. We improve on previous treat-
ments of tidal splicing by including spherical harmonic modes beyond the (2,2) mode, expanding the
Post-Newtonian expressions for tidal effects to 2.5 order, including dynamical tide corrections, and
adding a partial treatment of the spin-tidal dynamics. Furthermore, instead of numerical relativity
simulations, we use the spin-aligned BBH surrogate model “NRHybSur3dq8” to provide the BBH
waveforms that are input into the tidal slicing procedure. This allows us to construct spin-aligned,
inspiraling TaylorT2 and TaylorT4 splicing waveform models that can be evaluated quickly. These
models are tested against existing binary neutron star and black hole-neutron star simulations. We
implement the TaylorT2 splicing model as an extention to “NRHybSur3dq8”, creating a model that
we call “NRHybSur3dq8Tidal”.
I. INTRODUCTION
In August 2017, the network of interferometers con-
sisting of Advanced LIGO [1] and VIRGO [2] first ob-
served the gravitational radiation from the inspiral and
merger of a binary neutron star (BNS) [3], opening
the door to exploring extremely compact objects other
than binary black holes (BBH). Coincident detection
with the electromagnetic observational counterpart GRB
170817A [4], showed that these systems are the progen-
itors of short gamma-ray bursts, and herald the start
of multi-messenger astronomy. Additionally, this detec-
tion served as a probe of the neutron star equation of
state (EOS) [3, 5], and provided constraints on the grav-
itational wave speed [4]. As the detectors’ sensitivities
improve, observing more such systems will further con-
strain the governing physics [6, 7]. However, capturing
all of the information encoded in the measured signals
requires detailed, accurate templates that precisely de-
scribe waveforms from BNS or black hole-neutron star
(BHNS) systems.
A common approach to generate BHNS and BNS wave-
forms is to create analytic and phenomenological models
that capture the behavior of BHNS and BNS systems,
often in the form of additional corrections to BBH wave-
form models. Within the post-Newtonian (PN) formal-
ism, Ref [8] first computed the leading order tidal ef-
fects on the orbital evolution, characterized by the static
quadrupolar tidal deformability,
̄
λ
2
.
However recent
work suggests that the choice of the BBH model used as
the background for tidal corrections can impact parame-
ter estimation [9, 10], suggesting that errors in PN wave-
forms for BNS and BHNS systems might be dominated
by unknown higher-order vacuum terms rather than tidal
terms.
An alternative to the PN formalism is to implement
tidal corrections as an extension to the calibrated ef-
fective one-body (EOB) model SEOBNRv4 [11], includ-
ing additional static higher order effects [12–15]. This
leads to the time domain model SEOBNRv4T [14, 16].
A distinct EOB model, TEOBResumS [17], utilizes phe-
nomenological fits against NR results with tidal effects
that also include leading order self-spin effects.
Additional models have been constructed by calibrat-
ing to numerical simulations of BNS and BHNS binaries.
The LEA+ model [7] is a frequency domain waveform
model calibrated to a series of BHNS simulations with
mass ratios of
q >
2 [18] as an enhancement to the SEOB-
NRv2 model [19]. While LEA+ is a full waveform model,
it is valid only over the limited parameter space of the
simulations used to calibrate it,
i.e.
q >
2.
The
frequency
domain
models
SEOB-
NRv4
ROM
NRTidal, PhenomD
NRTidal, and Phe-
nomPv2
NRTidal [20] are built by combining the BBH
models SEOBNRv4
ROM [11], PhenomD [21, 22], and
PhenomPv2 [23, 24] with the NRTidal model [25].
NRTidal is a phenomenological fit of BNS/BHNS simu-
lation data to PN-like coefficients. While these models
cover a wide range of BNS parameter space, the fits are
calibrated to a limited number of waveforms and seem
to overestimate the tidal effects during the inspiral [20],
though improvements to this models are expected to
eliminate this problem [26].
The most accurate means of obtaining BNS and BHNS
waveform templates would be to run full numerical sim-
ulations. Running simulations that incorporate the rel-
evant matter physics for BHNS/BNS is a field of active
development [18, 25, 27–34]. However, the range of pos-
sible systems spans not just the masses and spins of the
components, but also includes all allowable EOSs. This
would require a large number of simulations to populate
such a high dimensional parameter space. Furthermore,
the large computational cost of such simulations makes
them impractical for parameter estimation purposes.
arXiv:1911.10440v1 [gr-qc] 24 Nov 2019
2
On the other hand, numerical simulations of BBH sys-
tems have made great strides over recent years, and there
are now public repositories of hundreds of simulations
for binaries with a variety of different initial masses and
spins [35–43]. Furthermore, surrogate models now allow
interpolation of numerical-relativity waveforms to desired
values of initial masses and spins [44–50]. Ref [49] showed
that surrogate models can robustly generate faithful rep-
resentations of binary black hole systems with spin mag-
nitudes
χ <
0
.
8 and masses low enough to be valid for
BNS systems (
M
≥
2
.
25
M
).
We build on a hybrid method called “tidal splic-
ing” [51], which computes inspiral waveforms for BNS
and BHNS by combining the accuracy of numerical BBH
simulations with the efficiency of PN models for tidally
deformable systems. This method does so by decompos-
ing the numerical BBH waveform in a manner akin to the
PN formalism and using this decomposition to replace all
orders of the vacuum terms in the analytic PN expansion
with their numerical equivalents. We combine these vac-
uum terms with the analytic tidal PN terms to build up
a waveform that models the inspiral of a BNS/BHNS.
In this paper, we continue the development of tidal
splicing beyond Ref. [51] by extending the method to
spinning systems and spherical harmonic modes beyond
the (2,2) mode. Using results from newly available EOB
models that incorporate higher PN effects [12–15], we
also extend the known higher order tidal effects from
EOB to the time domain PN approximants. Previous
explorations of tidal splicing [51] were based on particu-
lar individual numerical relativity BBH simulations, and
therefore could be tested only for the masses and spins
of those simulations. Here, instead of using numerical
relativity simulations directly, we will use the hybridized
surrogate model ‘NRHybSur3dq8’ [49] as our BBH base.
We organize this paper as follows: in Sec II we sum-
marize the current existing work on time domain tidal
waveforms in the PN framework; in Sec III we discuss
how we partially extend the PN tidal approximants to
2.5PN order and how we correct the dynamical tide ef-
fects for spinning NS; in Sec IV we explain our method
of tidal splicing; and in Sec V we compare tidal splicing
with some recent BHNS and BNS simulations.
Except where otherwise noted, we shall use the sub-
scripts
A,B
to refer to the individual NS or BH ob-
jects, the subscripts
`
= 2
,
3
,...
refer usually to the spe-
cific polar mode of the tidal effect in consideration (
i.e.
2=quadrupolar, 3=octopolar, etc.), while the
`,m
su-
perscripts will typically correspond to the spin-weighted
spherical harmonic modes
−
2
Y
`m
of the waveform. We
chose units of
G
=
c
= 1.
II. POST-NEWTONIAN THEORY
The PN approximation describes the binary’s orbital
behavior as series expansions that are valid in the slow-
moving, weak-field regime. The expansion parameter is
the characteristic velocity of the inspiraling objects,
v
(another common parameter is
x
=
v
2
). We denote an
expansion term of order
O
(
v
n
) by the label
n
2
PN (e.g.
2.5PN corresponds to
v
5
beyond leading order). A more
detailed summary of PN theory as it pertains to point-
particle systems can be found in [52].
We start with a quasi-circular binary system of a pair
compact objects with component masses
m
A
and
m
B
,
with total mass
M
, and spins of dimensionless magnitude
χ
A
and
χ
B
aligned with the orbital angular momentum.
Here
χ
A
=
S
A
/m
2
A
, where
S
A
is the spin angular mo-
mentum. We define the mass ratio
q
as the larger mass
over the smaller mass,
m
A
/m
B
, so that
q
≥
1. For con-
venience, we also define the mass fraction
X
A
=
m
A
/M
,
and symmetric mass ratio
ν
=
X
A
X
B
.
When the objects are not simply point particles, but
extended objects like neutron stars, each object responds
to the changing tidal fields. The leading tidal effects
are the result of the deformation of the NS due to the
tidal field generated by the other object in the binary.
This effect is characterized by the dimensionless
`
−
polar
tidal deformability parameter,
̄
λ
`
. Other commonly used
parameters are dimensionful tidal parameter
λ
`
or the
tidal love number
k
`
, and are related to
̄
λ
`
by the NS
radius
R
A
or compactness
C
A
=
m
A
/R
A
according to
̄
λ
`A
=
2
(2
`
−
1)!!
k
`A
C
2
`
+1
A
=
λ
`A
m
2
`
+1
A
.
(1)
As each object in the binary can have its own deforma-
bility, we add the subscript
A,B
to specify the particular
object.
In the following, we will often separate the PN expres-
sions into two parts: “BBH terms”, the terms that de-
scribe a BBH inspiral, and “tidal terms”, the terms that
describe corrections due to one or more of the objects be-
ing something other than a BH. This will be important
for tidal splicing, for which we replace the BBH terms
with numerical relativity (or a surrogate model thereof)
but we use the PN expressions for the tidal terms. The
tidal deformability
̄
λ
BH
of black holes is generally treated
as vanishing but is somewhat difficult to define [53]; here
we will set
̄
λ
BH
= 0, so all terms that depend on
̄
λ
`
are tidal terms. The BBH terms are identical to point-
particle terms up to 4PN for nonspinning BHs and 2.5PN
(with the exception of 2PN quadrupole moment terms)
for spinning BHs [54]. We include the 2PN correction for
spinning BHs in Section III A.
A. Orbital Evolution
Two equations govern the evolution of the quasi-
circular binary system in PN theory. The first relates
the orbital phase
φ
to
v
by a correspondence with the
orbital frequency,
ω
,
ω
=
dφ
dt
=
v
3
M
.
(2)
3
The other equation is the energy balance equation as
the emission of gravitational radiation drives the adia-
batic evolution by bleeding away the orbital energy,
E
(
v
).
If the energy flux is given by
F
(
v
), then this energy bal-
ance equation is
dE
(
v
)
dt
=
−
F
(
v
)
.
(3)
The quadrupolar tidal deformability
̄
λ
2
enters
E
(
v
)
and
F
(
v
) first at 0PN as a
O
(
v
10
) term , and through
1PN corrections at
O
(
v
12
) [8, 55]. While normally such
high order effects would be neglected, the relatively large
size of
̄
λ
2
∼ O
(1000) suggest the tidal deformations im-
pact the waveform earlier in the inspiral than expected
by their formal PN order.
Ref [8] provides the energy and flux expansions to 1PN,
E
(
v
) =
−
νv
2
2
[
1 +
(
−
3
4
−
ν
12
)
v
2
+
O
(
v
3
)
+
̄
λ
2
A
v
10
X
4
A
(
9(
−
1 +
X
A
) +
11
2
(
−
3 +
X
A
−
X
2
A
+ 3
X
3
A
)
v
2
+
O
(
v
3
)
)
+ (
A
→
B
)
]
,
(4)
F
(
v
) =
32
ν
2
v
10
5
[
1 +
(
−
1247
336
−
35
12
ν
)
v
2
+
O
(
v
3
)
+
̄
λ
2
A
v
10
X
4
A
(
6(3
−
2
X
A
) +
1
28
(
−
704
−
1803
X
A
+ 4501
X
2
A
−
2170
X
3
A
)
v
2
+
O
(
v
3
)
)
+ (
A
→
B
)
]
.
(5)
B. TaylorT Approximants
We can now insert the energy and flux expressions
into the orbital evolution, Eq. (2), and energy balance,
Eq. (3), equations and solve these equations to describe
the binary’s evolution. These equations are expanded in
powers of
v
, and then truncated at a particular order.
There are many choices of how to do this truncation,
and these choices give rise to different families of PN ap-
proximants. All these families agree to the same formal
PN order but have different higher order terms in
v
. The
two approximants that we will examine here are usually
referred to as TaylorT4 and TaylorT2.
Because the tidal terms are formally proportional to
at least
v
10
, naively truncating an expansion at a given
power of
v
would eliminate these terms until we reached
v
10
, where the point-particle terms are unknown. To en-
sure that tidal terms are included, and in light of the
fact that
̄
λ
2
is large, we handle the leading tidal terms
as if they were the same order as the leading PN terms:
O
(1)
∼O
(
̄
λ
2
v
10
) [8]. (See Appendix E for further discus-
sion regarding this correspondence in PN orders.) In the
expansions below we then keep all terms through 1PN
beyond leading order effects.
1. TaylorT4
The TaylorT4 [56] method generates the orbital evolu-
tion by rewriting the energy balance equation as
dv
dt
=
−
F
(
v
)
M
dE
(
v
)
dv
,
(6)
then expanding the ratio on the right hand side as a
power series in
v
and truncating at the appropriate order,
so that
dv
dt
=
F
BBH
(
v
) +
F
Tid
(
v
)
.
(7)
Here, we have broken the series into 2 parts: the terms
corresponding to a BBH system in
F
BBH
(
i.e.
̄
λ
2
= 0)
and the terms corresponding to the tidal correction in
F
Tid
. We don’t reproduce
F
BBH
here (as it is not needed
for our methods), but
F
Tid
to 1PN order is [8],
F
Tid
(
v
) =
32
νv
9
5
M
[
̄
λ
2
A
X
4
A
v
10
(
72
−
66
X
A
+
(
4421
56
−
12263
X
A
56
+
1893
X
2
A
4
−
661
X
3
A
2
)
v
2
)
+(
A
→
B
)
]
.
(8)
The TaylorT4 method computes the quantity
v
(
t
) by
integrating Eq. (7), then computes the orbital phase by
integrating
dφ
dt
=
v
3
/M.
(9)
The two constants arising from integrating both equa-
tions correspond to the inherent freedom to choose the
initial time and phase of the waveform.
2. TaylorT2
The TaylorT2 [57] expansion begins at the same point
as TaylorT4 with the PN energy equation and definition
of
v
, except the equations are rearranged to get a pair of
integral expressions parametric in
v
,
t
(
v
) =
t
0
+
M
∫
dE
(
v
)
dv
F
(
v
)
dv,
(10)
φ
(
v
) =
φ
0
+
∫
v
3
dE
(
v
)
dv
F
(
v
)
dv.
(11)
The integration constants
t
0
and
φ
0
are both freely speci-
fiable, and can be used to set the initial time and phase
of the resulting waveform.
The above integrands are expanded as a power series,
truncated to the appropriate order, and then integrated
to get series expressions for both the time and phase,
4
which we break into a part corresponding to a BBH sys-
tem and a part comprised of all the additional tidal ef-
fects,
t
(
v
) =
t
0
+
T
BBH
(
v
) +
T
Tid
(
v
)
,
(12)
φ
(
v
) =
φ
0
+
P
BBH
(
v
) +
P
Tid
(
v
)
.
(13)
As before, we do not reproduce expressions for
T
BBH
or
P
BBH
, but the 1PN tidal terms are
T
Tid
(
v
) =
−
5
M
256
νv
8
[
̄
λ
2
A
X
4
A
v
10
(
288
−
264
X
A
+
(
3179
4
−
919
X
A
4
−
1143
X
2
A
2
+ 65
X
3
A
)
v
2
)
+(
A
→
B
)
]
,
(14)
P
Tid
(
v
) =
−
1
32
νv
5
[
̄
λ
2
A
X
4
A
v
10
(
72
−
66
X
A
+
(
15895
56
−
4595
X
A
56
−
5715
X
2
A
28
+
325
X
3
A
14
)
v
2
)
+(
A
→
B
)
]
.
(15)
C. BBH Strain Modes
The gravitational radiation emission pattern for dis-
tant observers can be represented via a decomposition
into spin-weighted spherical harmonics. Following the
PN formalism used in [58], we express the strain from
compact objects inspiraling in quasi-circular orbits as
h
`m
BBH
(
v
) =
2
νv
2
M
r
√
16
π
5
H
`m
(
v
)
e
−
im
Ψ(
v
)
,
(16)
where
r
is the distance from the source to the detector.
The various terms of
H
`m
(
v
) are complex series expan-
sions of the individual modes and are distinct from the
series expansions for the energy and flux from above [58].
Ψ(
v
) is the tail-distorted orbital phase variable [59, 60]
Ψ(
v
) =
φ
(
v
)
−
2
Mω
ln
(
ω
ω
0
)
.
(17)
The constant
ω
0
is the reference frequency, often chosen
to be the frequency the waveform enters the detector’s
frequency band.
We will now rewrite Eq. (16) in a simpler and more
convenient form. The first simplification arises because
φ
(
v
) is proportional to
v
−
5
to leading order, and because
ω
is proportional to
v
3
. This means that the correction
in Eq. (17) is of 4PN order (i.e.
v
8
beyond leading order),
higher order than we consider in this paper, so we neglect
it and set Ψ(
v
) =
φ
(
v
).
The second simplification is to rewrite
H
`m
(
v
) in
Eq. (16), which is complex, in terms of real quantities.
We do this by treating the imaginary part of
H
`m
(
v
) as
a phase correction as is done in [61]. Thus we write
h
`m
BBH
(
v
) =
A
`m
BBH
(
v
)
e
i
(
ψ
`m
BBH
(
v
)
−
mφ
(
v
)
)
,
(18)
where
A
`m
BBH
(
v
) and
ψ
`m
BBH
(
v
) are real. This is the ex-
pression we will use for the waveform amplitudes in sub-
sequent sections.
We make one further simplification for the special case
of the (2
,
2) mode: for that mode,
ψ
2
,
2
BBH
can be neglected.
To see why, we examine the expression for
H
22
(
v
) [58]
and find that the first imaginary terms enter at 2.5PN
order. We then can write
H
22
(
v
)
e
−
2
iφ
(
v
)
=
A
22
BBH
(
v
)(1 +
iv
5
δ
+
O
(
v
6
))
e
−
2
iφ
(
v
)
≈
A
22
BBH
(
v
)
e
iv
5
δ
e
−
2
iφ
(
v
)
,
(19)
where
δ
is the imaginary 2.5PN coefficient of
H
22
(
v
).
Because
φ
(
v
) is proportional to
v
−
5
to leading order,
ψ
2
,
2
BBH
=
v
5
δ
is a 5PN phase correction (i.e. a correc-
tion to
φ
(
v
) that is
v
10
beyond leading order), so we set
ψ
2
,
2
BBH
= 0. Therefore
h
22
BBH
(
v
) =
A
22
BBH
(
v
)
e
−
2
iφ
(
v
)
.
(20)
We will see later that in the tidal splicing procedure that
replaces BBH terms with numerical relativity, Eq. (20)
allows us to extract
φ
(
v
) as the phase of the (2,2) mode.
D. Tidal Correction to Strain
Ref [13] computed the leading order PN tidal correc-
tions to the strain modes (these are explicitly written out
in the form we use here in Eqs. (A14-A17) of [62]). There
are no corrections to the phase of the individual modes
at leading order, (
i.e.
ψ
`m
Tid
(
v
) = 0), so the strain modes
for systems with tidally deformed objects are then
h
`m
Tid
(
t
) =
(
A
`m
BBH
(
v
) +
A
`m
Tid
(
v
)
)
e
i
(
ψ
`m
BBH
(
v
)
−
mφ
(
v
)
)
.
(21)
The additive corrections to the strain amplitudes are then
given by
A
22
Tid
(
v
) =
∣
∣
∣
∣
∣
24
√
π
5
v
12
̄
λ
2
A
X
5
A
(
3
−
5
X
A
+ 2
X
2
A
)
×
(
1 +
α
22
2
A
v
2
+
α
22
4
A
v
4
)
+ (
A
→
B
)
∣
∣
∣
∣
∣
,
A
21
Tid
(
v
) =
∣
∣
∣
∣
∣
8
√
π
5
v
13
̄
λ
2
A
X
5
A
(
9
2
−
15
X
A
+
33
X
2
A
2
−
6
X
3
A
)
×
(
1 +
α
21
2
A
v
2
)
−
(
A
→
B
)
∣
∣
∣
∣
∣
,
A
33
Tid
(
v
) =
∣
∣
∣
∣
∣
108
√
3
π
14
v
13
̄
λ
2
A
X
5
A
(
1
−
2
X
A
+
X
2
A
)
5
×
(
1 +
α
33
2
A
v
2
)
−
(
A
→
B
)
∣
∣
∣
∣
∣
,
A
31
Tid
(
v
) =
∣
∣
∣
∣
∣
12
√
π
70
v
13
̄
λ
2
A
X
5
A
(
1
−
2
X
A
+
X
2
A
)
×
(
1 +
α
31
2
A
v
2
)
−
(
A
→
B
)
∣
∣
∣
∣
∣
,
(22)
where
α
`m
i
are coefficients that depend on the masses.
The term
α
22
2
is
α
22
2
A
=
−
202 + 560
X
A
−
340
X
2
A
+ 45
X
3
A
42(3
−
2
X
A
)
,
(23)
while the rest of the
α
`m
i
are currently not known. The
corrections arising from
A
`m
Tid
not listed in Eq. (22) enter
at higher PN orders and so we ignore them here.
For
m
= odd modes, note that the (
A
→
B
) terms in
Eqs. (22) appear with an overall minus sign. This can
be understood by considering that for identical objects
A
and
B
, the
m
=odd modes must vanish because of
symmetry under an azimuthal rotation by
π
.
E. Dynamical Tides
The tidal corrections considered so far are based on a
tidal deformability parameter, which describes the defor-
mation of a stationary object in the presence of a sta-
tionary tidal field. For dynamical objects in a binary
system, this amounts to treating the tidal deformation
as proportional to the instantaneous tidal field of the
companion. However, the objects also have internal
f
-
modes described by the resonant frequencies
ω
f`
. In late
inspiral, as the orbital frequency approaches
ω
f`
, it is
no longer appropriate to treat tides as stationary, and
dynamical tidal effects must be considered. Ref [14, 16]
analyzed how these dynamical tides affect the orbital mo-
tion for nonspinning systems. The approximate solution
they derive treats the dynamical tidal deformabilities as
frequency-dependent scalings of their static values. We
summarize their results here.
In effect, the dynamical tides serve primarily to am-
plify the static deformability during the evolution, peak-
ing when the orbital frequency is on resonance with one
of the object’s internal
f
−
modes, with frequency
ω
f`A
.
We shall also make use of the dimensionless
f
−
mode res-
onance frequency,
̄
ω
f`A
=
Mω
f`A
.
(24)
Note that our choice of defining ̄
ω
f`A
by scaling it as
the binary’s total mass
M
, rather than the object mass
m
A
as other works often use, is for our convenience when
comparing with the dimensionless orbital frequency, ̄
ω
=
Mω
.
In the nonspinning case, we denote the characteristic
parameter governing the resonance as
γ
`A
=
`ω
ω
f`A
=
`v
3
̄
ω
f`A
.
(25)
This parameter characterizes how close the system is to
resonance.
Dynamical tidal effects result in a multiplicative cor-
rection factor for the deformability. There are two such
correction factors: one that appears in the orbital evolu-
tion equations and one that appears in the strain ampli-
tudes.
For the orbital phase evolution, Ref [14] expressed the
effective enhancement factor
κ
`A
(
v
) on the static de-
formability of object
A
by
κ
`A
(
v
) =
a
`
+
b
`
[
1
1
−
γ
2
`A
+
4
3
√
`
ˆ
t
`
γ
2
`A
+
√
π
3
`
Q
`
γ
2
`A
]
,
(26)
where
`
=
256
νω
5
/
3
f`A
5
`
5
/
3
,
ˆ
t
`
=
8
5
√
`
(
1
−
γ
−
5
/
3
`A
)
,
Q
`
= cos
(
3
ˆ
t
`
8
)
[
1 + 2F
S
(
√
3
2
√
π
ˆ
t
`
)]
−
sin
(
3
ˆ
t
`
8
)
[
1 + 2F
C
(
√
3
2
√
π
ˆ
t
`
)]
,
(27)
with F
S
and F
C
as the Fresnel sine and cosine integrals.
The coefficients (
a
`
,b
`
) are given by (
a
2
,b
2
) = (1
/
4
,
3
/
4)
and (
a
3
,b
3
) = (3
/
8
,
5
/
8). The PN orbital evolution equa-
tions, Eqs. (8, 14), and (15), can be modified to incorpo-
rate the dynamical tides by taking
̄
λ
`A
→
̄
λ
`A
κ
`A
(
v
).
The correction to the strain amplitudes is computed in
Ref. [15]:
ˆ
κ
2
A
(
v
) =
(
κ
2
A
(
v
)
−
1)
ω
2
f
2
A
+ 6(1
−
X
A
)
κ
2
A
ω
2
(9
−
6
X
A
)
ω
2
.
(28)
As in the case of the deformability amplification for the
orbital evolution, the deformability amplification for the
dynamical strain can be incorporated into Eq. (22) by
the substitution
̄
λ
2
A
→
̄
λ
2
A
ˆ
κ
2
A
(
v
).
Note that these dynamical tide corrections do not pro-
duce proper power series expansions of
v
, because the
Fresnel integrals do not have a well defined power series
expansion about
v
= 0 (even though those terms vanish
as
v
→
0, they also oscillate infinitely fast). So including
κ
`A
(
v
) and ˆ
κ
2
A
(
v
) into the evolution and strain formula
means the tidal effects cannot be represented to a for-
mal PN order. While this is not a problem for the actual
generation of the waveforms, it does mean that the differ-
ences between the different families of PN approximants
no longer diverge at a well defined PN order.
6
III. ADDITIONAL TIDAL CORRECTIONS
A. Partial 2.5PN Tidal Terms
The static tidal corrections to the orbital evolution
discussed so far have been those derived from Eqs. (4)
and (5), which are 1PN order expressions computed
by Ref. [8]. Since then, higher-order corrections have
been computed, but only within the EOB formalism and
not for standard PN approximants. These higher-order
terms include not only nonspinning
̄
λ
2
tidal corrections to
the energy (through the EOB Hamiltonian) up through
2.5PN order, but also corrections according to the static
octopolar deformability
̄
λ
3
[12] through 2.5PN order.
These 2.5PN tidal effects are already included both for
SEOBNRv4T [14, 16], and for the frequency domain ap-
proximant TaylorF2 [62]. We use the EOB results to
obtain the time domain Taylor approximants here.
The details of how to convert the
̄
λ
2
and
̄
λ
3
2.5PN
tidal terms from the EOB Hamiltonian to corrections
to Eq. (4) are given in Appendix A. That computation
yields
E
̄
λ
2
(
v
) =
−
νv
2
2
[
̄
λ
2
A
v
10
X
4
A
(
9(
−
1 +
X
A
)
+
11
2
(
−
3 +
X
A
−
X
2
A
+ 3
X
3
A
)
v
2
+
13
2
(
−
51
4
−
15
4
X
A
+
361
42
X
2
A
+
47
21
X
3
A
+
47
12
X
4
A
+
7
4
X
5
A
)
v
4
+
O
(
v
6
)
)
+ (
A
→
B
)
]
,
(29)
E
̄
λ
3
(
v
) =
−
νv
2
2
[
̄
λ
3
A
v
14
X
6
A
(
(
−
65 + 65
X
A
)
+
(
75
2
−
875
2
X
A
+
475
2
X
2
A
+
325
2
X
3
A
)
v
2
+
(
−
6205
24
−
4165
8
X
A
+
425
36
X
2
A
−
1955
4
X
3
A
+
26095
24
X
4
A
+
12155
72
X
5
A
)
v
4
+
O
(
v
6
)
)
+(
A
→
B
)
]
.
(30)
From these expressions we can see that the leading order
octopolar deformability terms enter at
O
(
v
14
). Numeri-
cally,
̄
λ
3
is typically larger than
̄
λ
2
, with
̄
λ
3
∼O
(1000
−
10000). Thus, as we did for
̄
λ
2
, we treat the leading order
̄
λ
3
terms as being the same formal order as the leading or-
der PN terms,
i.e.
O
(1)
∼O
(
̄
λ
2
v
10
)
∼O
(
̄
λ
3
v
14
). Terms
for both deformabilties are computed up to 2.5PN order.
Unfortunately, unlike the energy, the fluxes are not
known to 2.5PN order. We shall introduce undefined
coefficients for the missing terms,
α
4
,β
0
,β
2
,β
4
so that
we may expand the PN expressions to a consistent order.
We write the flux terms as
F
̄
λ
2
(
v
) =
32
ν
2
v
10
5
[
̄
λ
2
A
v
10
X
4
A
(
6(3
−
2
X
A
)
+
1
28
(
−
704
−
1803
X
A
+ 4501
X
2
A
−
2170
X
3
A
)
v
2
+
α
4
v
4
+
O
(
v
6
)
)
+ (
A
→
B
)
]
,
(31)
F
̄
λ
3
(
v
) =
32
ν
2
v
10
5
[
̄
λ
3
A
v
14
X
6
A
(
β
0
+
β
2
v
2
+
β
4
v
4
+
O
(
v
6
)
)
+ (
A
→
B
)
]
.
(32)
We
keep
track
of
the
undefined
coefficients
α
4
,β
0
,β
2
,β
4
for the purposes of completeness; when the
value of those coefficients are computed in some future
work, those values can simply substituted into the flux
expression and the Taylor expansions of the orbital
evolution below. For our results in Section V, we set
these coefficients to zero.
For aligned spin systems, there will also be terms cor-
responding to spin-tidal terms, interactions in the Hamil-
tonian between object spins and the tidal deformations
of the NS. The spin-tidal connection terms are not in-
cluded in our expansions because of discrepencies in the
literature regarding the calculation of the leading order
1.5PN coefficients; a discussion of those coefficients, in-
cluding how to add them when those discrepencies are
resolved, can be found in Appendix D.
Because we have added tidal terms up to 2.5PN or-
der, for a consistent expansion we also introduce non-
tidal 2.5PN spin-orbit, spin-spin, and rotationally in-
duced quadrupolar moment effects [63–72] to both the
orbital energy and the flux expressions. The non-tidal
parts of these expressions are:
E
BBH
(
v
) =
−
νv
2
2
[
1 +
(
−
3
4
−
ν
12
)
v
2
+
(
2
χ
A
X
A
(
1 +
X
A
3
)
−
2
χ
B
X
B
(
1 +
X
B
3
))
v
3
+
(
1
8
(
−
27 + 19
ν
−
ν
2
3
)
+ 2
χ
A
χ
B
ν
−
(
̄
Q
A
+ 1)
χ
2
A
X
2
A
−
(
̄
Q
B
+ 1)
χ
2
B
X
2
B
)
v
4
+
(
χ
A
X
A
(
3 +
5
X
A
3
+
29
X
2
A
9
+
X
3
A
9
)
+
χ
B
X
B
(
3 +
5
X
B
3
+
29
X
2
B
9
+
X
3
B
9
))
v
5
+
O
(
v
6
)
]
,
(33)
7
F
BBH
(
v
) =
32
ν
2
v
10
5
[
1 +
(
−
1247
336
−
35
12
ν
)
v
2
+
(
4
π
+
χ
A
X
A
(
−
5
4
−
3
X
A
2
)
+
χ
B
X
B
(
−
5
4
−
3
X
B
2
))
v
3
+
((
−
44711
9072
+
9271
ν
504
+
65
ν
2
18
)
+
31
8
χ
A
χ
B
ν
+
(
33
16
+ 2
̄
Q
A
)
χ
2
A
X
2
A
+
(
33
16
+ 2
̄
Q
B
)
χ
2
B
X
2
B
)
v
4
+
((
−
8191
672
−
583
ν
25
)
π
+
χ
A
X
A
(
−
13
16
+
63
X
A
8
−
73
X
2
A
36
−
157
X
3
A
18
)
+
χ
B
X
B
(
−
13
16
+
63
X
B
8
−
73
X
2
B
36
−
157
X
3
B
18
))
v
5
+
O
(
v
6
)
]
.
(34)
In the above expressions,
̄
Q
A
is the dimensionless
quadrupole moment, which is unity for black holes,
̄
Q
BH
= 1, and is related to the dimensionful quadrupole
moment
Q
A
by
̄
Q
A
=
−
Q
A
m
3
A
χ
2
A
.
(35)
Therefore the full 2.5PN energy and flux expressions
(replacing the 1PN versions given by Eqs. (4) and (5))
are simply
E
(
v
) =
E
BBH
(
v
) +
E
̄
λ
2
(
v
) +
E
̄
λ
3
(
v
)
,
(36)
F
(
v
) =
F
BBH
(
v
) +
F
̄
λ
2
(
v
) +
F
̄
λ
3
(
v
)
.
(37)
At this point, we can repeat the expansion procedure
from Sections II B 1 and II B 2 to generate the TaylorT4
and TaylorT2 approximants. Again, we treat terms like
O
(1)
∼ O
(
̄
λ
2
v
10
)
∼ O
(
̄
λ
3
v
14
) as leading order and ex-
pand through 2.5PN order.
1. TaylorT4
With the 2.5PN expressions for the energy and flux,
we can now recompute the TaylorT4 approximant from
Eq. (7) up to 2.5PN order. As before, the terms in
dv/dt
corresponding to a BBH system (
̄
λ
2
=
̄
λ
3
= 0
,
̄
Q
= 1)
are denoted
F
BBH
(
v
) (which we do not reproduce here),
and the terms describing tidal corrections are denoted
F
Tid
(
v
). Here
F
Tid
(
v
) =
32
νv
9
5
M
[
(
5(
̄
Q
A
−
1)
χ
2
A
X
2
A
)
v
4
+
̄
λ
2
A
X
4
A
v
10
(
5
∑
i
=0
F
2
A,i
v
i
)
+
̄
λ
3
A
X
6
A
v
14
(
5
∑
i
=0
F
3
A,i
v
i
)
+ (
A
↔
B
)
]
.
(38)
The coefficients
F
2
A,i
and
F
3
A,i
are given in Appendix B,
Eq. (B1).
The
̄
λ
2
×
χ
A,B
cross terms appearing in these coeffi-
cients are not due to spin-tidal interaction terms in the
Hamiltonian, but instead are a consequence of the se-
ries expansion power counting. The (
̄
Q
A
−
1) in the
v
4
term arises from the fact that we have incorporated the
BBH part of
̄
Q
A
into
F
BBH
(
v
) (for which
̄
Q
BH
= 1) and
included the remainder here.
2. TaylorT2
Similarly for TaylorT2, the updated time and phase
expressions corresponding to Eqs. (14–15) are
T
Tid
(
v
) =
−
5
M
256
νv
8
[
−
(
10(
̄
Q
A
−
1)
χ
2
A
X
2
A
)
v
4
+
̄
λ
2
A
X
4
A
v
10
(
5
∑
i
=0
T
2
A,i
v
i
)
+
̄
λ
3
A
X
6
A
v
14
(
5
∑
i
=0
T
3
A,i
v
i
)
+ (
A
↔
B
)
]
,
(39)
P
Tid
(
v
) =
−
1
32
νv
5
[
−
(
25(
̄
Q
A
−
1)
χ
2
A
X
2
A
)
v
4
+
̄
λ
2
A
X
4
A
v
10
(
5
∑
i
=0
P
2
A,i
v
i
)
+
̄
λ
3
A
X
6
A
v
14
(
5
∑
i
=0
P
3
A,i
v
i
)
+ (
A
↔
B
)
]
.
(40)
The individual coefficients
T
2
A,i
,
T
3
A,i
,
P
2
A,i
,
P
3
A,i
are
given in Appendix B, Eqs. (B2–B3).
B. Spinning Dynamical Tides
The dynamical tidal corrections discussed earlier in
Section II E do not take into account the spin of the NS.
The dynamical tides are caused by the changing tidal
field due to the orbital motion interacting with the in-
ternal
f
−
modes of the deformable object. So when the
object is also spinning, then its internal modes will effec-
tively experience a driving frequency equal to the orbital
frequency shifted by the object’s rotational frequency.
We characterize this frequency shift by making a slight
correction to the characteristic parameter
γ
`A
in Eq. (25).
Given an aligned or antialigned spin of
χ
A
and a mo-
ment of inertia
I
A
, then we can compute the rotation
8
FIG. 1. The dynamical tide amplification
κ
`
from Eq. (26) modified by Eq. (43), as a function of orbital frequency, for the
nonspinning case (black) and the spin aligned/antialigned (blue/red) cases for both
`
= 2 (solid) and
`
= 3 (dashed) tidal
deformabilities. The vertical lines represent the resonance frequencies of both modes in every case. The parameters correspond
to an NS with
̄
λ
2
≈
800 in an equal mass system with total mass 2
.
8
M
. For this NS, a spin magnitude of
|
χ
NS
|
=
.
2
corresponds to a rotational frequency of
f
NS
≈
312Hz.
frequency of the deformable object as
̄
ω
A
=
Mω
A
=
χ
A
X
A
̄
I
A
,
(41)
where the dimensionless moment of inertia is
̄
I
A
=
I
A
m
3
A
.
(42)
Thus, the effective orbital frequency the resonant
modes will experience is simply the difference between
the orbital frequency and the rotational frequency of the
object. With that in mind, we rewrite
γ
`A
as
γ
`A
=
∣
∣
∣
∣
∣
`
(
v
3
−
̄
ω
A
)
̄
ω
f`A
∣
∣
∣
∣
∣
.
(43)
We take the absolute value because Eq. (25) is undefined
for negative values of
γ
`A
, which can occur at low orbital
frequencies with aligned spin objects. This corresponds
to saying that the resonance modes of the object only care
about the magnitude of the frequency of the changing
tidal field. Within Eq. (28), we also need to make the
change
ω
→
∣
∣
v
3
−
̄
ω
A
∣
∣
/M
.
To show how spin affects the profile of the dynamical
tide correction, in Fig 1 we plot the profile of
κ
`
from
Eq. (26) as a function of orbital frequency. We assume
an NS with
̄
λ
2
≈
800 in an equal mass system and use the
universal relations (see Table II) to compute the other rel-
evant tidal parameters. We compare the nonspinning NS
against both aligned and antialigned spinning NS with
magnitude
|
χ
NS
|
=
.
2, which corresponds to a rotational
frequency of
f
NS
≈
312Hz. The upper frequency termina-
tion point is at an orbital frequency of
Mω
ISCO
= 6
−
3
/
2
,
which has been used as an approximate BNS inspiral ter-
mination criterion [73].
From Fig 1 we can see that aligned spins push the res-
onance peak later in the inspiral while antialigned spins
move the peak to smaller frequencies. In fact, with large
enough aligned spins it is possible that the peak of the
resonance is never reached before the system enters the
merger/ringdown phase.
At low frequencies, the nonspinning system behaves
the same as for static tides (
κ
`
= 1). However, this is not
true for the spinning systems, both of which asymptote
to a slightly different value. Physically, these differences
are due to the deformations of the object experiencing a
driving frequency not from the orbital frequency (which
is vanishingly small), but from its own rotation along its
axis. We expect this difference to result in a negligible
contribution to the waveform as the relative size of the
tidal effects already vanishes (
O
(
v
10
)) at low frequencies.
For the aligned spin object, note that there is a point in
the evolution where the rotational and orbital frequencies
match and the effective dynamical tidal field vanishes.
One concern is that for antialigned spins,
κ
`
becomes
negative before
Mω
ISCO
. We recognize that Eq. (26)
is derived assuming that the driving frequency is not
much larger than the resonance peak [14]. While for
nonspinning systems and spin aligned systems the res-
9
onance peak occurs near the end of the inspiral or after
merger/ringdown and is thus within the range of valid-
ity, that condition does not necessarily apply in the an-
tialigned case. If the antialigned spin is large enough, as
seems to be the case in Fig 1, the resonance frequency
occurs early enough in the evolution that this approxi-
mate formalism potentially breaks down while still in the
inspiral, necessitating a more delicate handling of the dy-
namical tides.
Until such a formalism is developed for antialigned NS
spins, we instead assume for antialigned spin that the
object is nonspinning (
i.e.
we set
ω
A
= 0) for the pur-
poses of Eqs. (26) and (28); the aligned spin case will use
Eq. (43) as expected.
IV. TIDAL SPLICING
Putting together the results of the previous sections,
we can write the PN equations of motion in the form
dφ
dt
=
v
3
M
,
(44)
dv
dt
=
F
BBH
(
v
) +
F
Tid
(
v
)
,
(45)
h
`m
(
t
) =
(
A
`m
BBH
(
v
) +
A
`m
Tid
(
v
)
)
e
i
(
ψ
`m
BBH
(
v
)
−
mφ
)
,
(46)
where we have expressed the equations within the Tay-
lorT4 framework. The expression for
F
Tid
(
v
) is given by
Eq. (38),
F
BBH
(
v
) can be computed with the mathemat-
ica notebook attached to this paper, while
A
`m
BBH
(
v
) and
ψ
`m
BBH
(
v
) can be extracted from the expansions of Eq. (16)
given in [58] using the procedure we describe in Sec II C.
However,
F
BBH
(
v
),
A
`m
BBH
(
v
) and
ψ
`m
BBH
(
v
) are all unim-
portant for our purposes and so we do not present those
formulae here.
Equivalently, the PN equations of motion can be writ-
ten within the TaylorT2 framework as
t
(
v
) =
t
0
+
T
BBH
(
v
) +
T
Tid
(
v
)
,
(47)
φ
(
v
) =
φ
0
+
P
BBH
(
v
) +
P
Tid
(
v
)
,
(48)
h
`m
(
t
) =
(
A
`m
BBH
(
v
) +
A
`m
Tid
(
v
)
)
e
i
(
ψ
`m
BBH
(
v
)
−
mφ
)
.
(49)
The expression for
T
Tid
(
v
) and
P
Tid
(
v
) are given by
Eqs. (39–40).
Consider the case of a BBH system. In principle, if
we knew the PN expansions for the expressions of
F
BBH
(or
P
BBH
and
T
BBH
),
A
`m
BBH
and
ψ
`m
BBH
up through ar-
bitrarily large order, then we could perfectly reproduce
the gravitational waveforms of those inspiraling systems.
But unfortunately these terms are known only to a lim-
ited order.
However, numerical simulations of BBH systems are
able to accurately solve the full Einstein equations. Thus,
if we can represent these numerical waveforms in a form
akin to the systems of equations in either Eqs. (44–
46) or Eqs. (48–49) (with vanishing tidal terms), they
would provide perfect representations of the PN ex-
pressions up to numerical resolution error. This is the
main idea of tidal splicing: We use numerical relativ-
ity BBH waveforms to effectively obtain the functions
F
BBH
(
v
),
P
BBH
(
v
),
T
BBH
(
v
),
A
`m
BBH
(
v
) and
ψ
`m
BBH
(
v
) to all
orders in
v
. Then we add the analytic expressions for the
tidal terms described above, and we integrate the PN
equations of motion to generate waveforms correspond-
ing to the inspirals of BHNS and BNS systems.
A. Decomposition of NR waveforms
We start with a numerical BBH waveform correspond-
ing to a system with a particular mass ratio
q
and spins
χ
A
and
χ
B
. The spins are assumed here to be aligned
or antialigned with the orbital angular momentum. We
decompose the (2
,
2) mode of the BBH waveform in the
form
h
22
NR
(
t
) =
A
22
NR
(
t
)
e
−
2
iφ
NR
(
t
)
,
(50)
where
A
22
NR
(
t
) is real. This is the same decomposition as
Eq. (20), so we interpret
φ
NR
(
t
) as the NR orbital phase.
We then compute the effective PN expansion parameter
v
NR
using Eq. (9):
v
NR
(
t
) =
3
√
M
dφ
NR
(
t
)
dt
.
(51)
Because the numerical waveform is known at a finite
set of time samples, we use a 6th order finite difference
scheme to compute the
d/dt
derivatives numerically.
With the orbital phase and effective PN parameters in
hand, we decompose each mode from the waveform into
an amplitude and phase,
h
`m
NR
(
t
) =
A
`m
NR
(
t
)
e
i
Φ
`m
NR
(
t
)
.
(52)
Comparing with the expression for strain from Eq. (18),
we break up the phase as
Φ
`m
NR
(
t
) =
ψ
`m
NR
(
t
)
−
mφ
NR
(
t
)
.
(53)
Since we know
φ
NR
we can compute
ψ
`m
NR
from the total
phase by rearranging Eq. (53),
ψ
`m
NR
(
t
) = Φ
`m
NR
(
t
) +
mφ
NR
(
t
)
.
(54)
Up to this point, we have been treating
t
NR
as the
independent variable for the purposes of decomposition.
Since the PN formalism considers the frequency expan-
sion parameter
v
as the independent variable, we invert
v
NR
(
t
) to get
t
NR
(
v
) as a function of
v
. This inversion is
straightforward numerically, since we have a finite num-
ber of time samples
t
NR
. Thus we can represent all of the
individual parts of our waveform as functions of
v
, e.g.
φ
NR
(
t
NR
(
v
)) =
φ
NR
(
v
). Then we can write the strain for
each mode in the form corresponding to Eq. (18):
h
`m
NR
(
v
) =
A
`m
NR
(
v
)
e
i
(
ψ
`m
NR
(
v
)
−
mφ
NR
(
v
))
.
(55)
10
We now have numerical equivalents for the various PN
expansions for BBH systems; these are correct up to an
arbitrary PN order and limited only by the errors from
the simulations themselves.
B. Tidal parameters
With the numerical decomposition in hand, we will
need the tidal parameters for the particular BHNS or
BNS system under consideration.
A review of the
different tidal effects explored in the previous section
show there are 6 different parameters that character-
ize the tidal behavior of each object: the dimension-
less quadrupolar and octopolar tidal deformabilities
̄
λ
2
and
̄
λ
3
, their corresponding
f
−
mode resonance frequen-
cies
ω
f
2
and
ω
f
3
, the dimensionless rotationally-induced
quadrupole moment
̄
Q
, and the dimensionless moment
of inertia
̄
I
. For our model, we choose only
̄
λ
2
, and we
compute the values of the other parameters from
̄
λ
2
using
the universal relations, which are approximate relations
between
̄
λ
2
and the other tidal parameters; the details
are given in Appendix C. The choice of
̄
λ
2
depends on
the physical properties of the deformable object in con-
sideration. Once we have chosen
̄
λ
2
(and thus the other
parameters via the universal relations), the next step in
tidal splicing is the recomputation of the orbital evolu-
tion.
C. Splicing of the orbital evolution equations
The details of how to specifically splice the PN tidal
information into the orbital evolution depends on the
specific Taylor expansion considered. We shall discuss
the details of tidal splicing with TaylorT4 and TaylorT2.
This method can also be performed for TaylorT3 [57]
which involves expanding about an intermediate dimen-
sionless time variable. However, TaylorT3 is known to do
a poor job in general of reproducing the results of BBH
numerical simulations even in the equal mass, nonspin-
ning case [56, 74], so we ignore that method here. We
also ignore TaylorT1 [57] because we do not know how
to compute the BBH contribution to the PN energy and
flux separately using only the BBH waveform. While
Ref. [51] discussed tidal splicing under a TaylorF2 [75]
framework, in this paper we only examine time domain
approximants.
1. TaylorT4 Splicing
Splicing with TaylorT4, originally introduced within
[51], begins by examining how the tidal terms manifest
in the TaylorT4 framework. The evolution of the PN
parameter as seen in Eq. (45) for a BBH system is
dv
dt
=
F
BBH
(
v
)
.
(56)
With
t
NR
(
v
) in hand from the simulation, we can com-
pute a numerically accurate version of
F
BBH
(
v
) which we
shall call
F
NR
(
v
),
F
NR
(
v
) =
(
dt
NR
(
v
)
dv
)
−
1
(57)
The tidal terms,
F
Tid
(
v
NR
), represent the sum of the
additional tidal effects in the evolution and we compute
them according to Eq. (B1). We incorporate the dy-
namical tides by scaling the deformabilities according to
Eq. (26),
i.e.
̄
λ
`
→
̄
λ
`
κ
`
(
v
).
We introduce a new spliced time coordinate,
t
Spl
(
v
),
which we compute by integrating the differential equation
dt
Spl
dv
=
1
F
NR
(
v
) +
F
Tid
(
v
)
,
(58)
which is the inverse of Eq. (45). Once we have
t
Spl
, we
find the orbital phase of this new waveform by integrating
Eq. (44):
φ
Spl
(
v
) =
1
M
∫
v
(
t
Spl
)
3
dt
Spl
.
(59)
We use Simpson’s method for integration. We choose the
integration constants in Eqs. (58) and (59) to align the
waveform to the numerical waveform at the initial time.
2. TaylorT2 Splicing
For BBH systems Eqs. (47–48) take the form
t
(
v
) =
t
0
+
T
BBH
(
v
)
,
φ
(
v
) =
φ
0
+
P
BBH
(
v
)
.
(60)
The constants
t
0
and
φ
0
correspond simply to the start-
ing time and phase of the waveform. Our numerically
corrected versions of
T
BBH
and
P
BBH
are simply
T
NR
(
v
) =
t
NR
(
v
)
,
P
NR
(
v
) =
φ
NR
(
v
)
.
(61)
We compute
T
Tid
(
v
) according to Eq. (B2) and
P
Tid
(
v
)
according to Eq. (B3), incorporating the dynamical tides
by making the frequency dependent adjustment to
̄
λ
`
from Eq. (26).
The spliced waveform’s time
t
Spl
and phase
φ
Spl
are
then given by examining Eqs. (47–48) and making the
appropriate substitutions,
t
Spl
(
v
) =
t
0
+
t
NR
(
v
) +
T
Tid
(
v
)
,
φ
Spl
(
v
) =
φ
0
+
φ
NR
(
v
) +
P
Tid
(
v
)
.
(62)
We use the freedom inherent in choosing
t
0
and
φ
0
to
align the spliced waveform to the numerical waveform at
the initial time
As the waveform nears the merger phase of the evolu-
tion, the effect from
T
Tid
(
v
) might grow larger than that
of
t
NR
(
v
). This may cause
t
Spl
(
v
) to become nonmono-
tonic at some
v
; if this happens, we end the waveform at
that value of
v
.
11
D. Waveform Reconstruction
Once
t
Spl
(
v
) and
φ
Spl
(
v
) are computed, then the fi-
nal step is reconstructing the spliced waveform from
Eqs. (46) or (49) (which are the same equation). For
A
`m
BBH
(
v
) and
ψ
`m
BBH
(
v
) we use
A
`m
NR
(
v
) and
ψ
`m
NR
(
v
) com-
puted from Eq. (55). For
A
`m
Tid
(
v
) we use the expressions
in Eq. (22), with the dynamical tides accounted for by
the replacement rule
̄
λ
2
→
̄
λ
2
A
ˆ
κ
2
A
(
v
) from Eq. (28). We
thus arrive at the final formula for the spliced waveform
modes:
h
`m
Spl
(
v
) =
(
A
`m
NR
(
v
) +
A
`m
Tid
(
v
)
)
e
i
(
ψ
`m
NR
(
v
)
−
mφ
Spl
(
v
))
.
(63)
To get a time-domain waveform, we invert the function
t
Spl
(
v
). Because
v
is known only at discrete values, we
interpolate the amplitudes and phases of the waveforms
onto a set of uniformly-spaced values of
t
Spl
using a cubic
spline.
V. RESULTS
A. Models for Comparison
To measure the accuracy of the tidal splicing method,
we compare our spliced waveforms against numerical sim-
ulations of BHNS/BNS inspirals. In particular, we use
some of the recent numerical simulations from [33]. In
all simulations we compare against, the NSs were gen-
erated according to an EOS of Γ = 2 polytrope with a
mass
M
ADM
= 1
.
4
M
and compactness of
C
NS
= 0
.
1444
so that the quadrupolar tidal deformability is
̄
λ
2
∼
800.
Comparing against one specific EOS at a single NS mass
is a small slice of the full possible BHNS/BNS parame-
ter space, but should give some idea for how well tidal
splicing can perform. See Table I for the full list of the
simulations we consider here. In particular, there are 1
BNS and 5 BHNS runs, and two of the BHNS runs have
a small anti-aligned spin on the NS while the rest of the
runs have zero spins.
We generated our tidally spliced waveforms for each
of these cases using the hybridized surrogate model
‘NRHybSur3dq8’ [49] to compute the underlying BBH
signal and making use of the universal relations to ob-
tain the other tidal parameters from
̄
λ
2
; the details are
given in Appendix C. The TaylorT2 splicing model has
also been implemented as a tidal extension of ‘NRHyb-
Sur3dq8’ under the name “NRHybSur3dq8Tidal”.
To provide an additional point of comparison, we also
test another waveform model, SEOBNRv4T, which is
the time domain model SEOBNRv4 [11] augmented with
most of the same effects we have used here, including
higher order corrections to the static tides in the EOB po-
tential [12], strain corrections [13], and dynamical tides
(without the resonance frequency correction for spinning
NSs) [16]. These effects correspond to our Eq. (A2),
Eq. (22), and Eqs. (26) and (28).
Type
q
χ
NS
f
0
(Hz)
f
1
(Hz)
N
22
cyc
BHNS
1
0
218
578
19.9
BNS
1
0
211
629
20.8
BHNS
1
-.2
217
505
17.0
BHNS
1.5
0
154
537
28.9
BHNS
2
0
156
505
21.0
BHNS
2
-.2
156
485
19.8
TABLE I. List of parameters for numerical simulations
from [33] considered in this paper, namely the mass ratio
q
, the dimensionless spin of the neutron star
χ
NS
, the lower
and upper orbital frequency cutoff of the waveform
f
0
and
f
1
,
and the number of cycles in the (2,2) mode of the waveform
between the cutoff frequencies.
B. Waveforms
The numerical waveforms we compare against include
all modes through
`
= 5, while the surrogate and spliced
waveforms model only the modes [(2, 0), (2, 1), (2, 2),
(3, 2), (3, 0), (3, 1), (3, 3), (4, 2), (4, 3), (4, 4), (5, 5)],
and SEOBNRv4T models only the (2,2) mode. Since
all systems we consider have spins parallel to the orbital
angular momentum, we need only
m
≥
0 modes; the
m <
0 modes are obtained from symmetry. The strain
measured along a particular direction in the sky can then
be written as
h
+
(
ι,φ
0
)
−
ih
×
(
ι,φ
0
) =
∑
`m
h
`m
−
2
Y
`m
(
ι,φ
0
)
(64)
Here
−
2
Y
`m
are spin-weighted spherical harmonics, and
the angles (
ι,φ
0
) are defined so that
ι
is the inclination
angle between the binary angular momentum and the
line of sight to the observer while
φ
0
is the binary orbital
phase when it enters the detector sensitivity band.
We choose the beginning of the waveform to be at a
time after the initial burst of junk radiation,
t
= 200
M
.
That time also sets the starting orbital frequency of the
waveform (chosen as half the time derivative of the phase
of the (2,2) mode). To prevent the starting frequency
f
Initial
from being contaminated by residual junk radia-
tion in the imperfect BHNS/BNS initial data, we use a
quadratic fit of the simulations’ frequency against time
over the interval
t
∈
(200
M,
700
M
) to estimate the pre-
cise starting frequency. We window the waveform with a
Planck-Taper window [76] over that early 500
M
region
of the waveform. We label the orbital frequency at the
end of this window as
f
0
; this frequency will serve as the
initial frequency considered in our mismatches below.
At late times, we set the upper frequency cutoff
ω
Cutoff
by the frequency attained by the simulation at its peak
power, and window the waveform (again with Planck-
Taper) over the times from
f
1
=
.
85
f
Cutoff
to
f
Cutoff
.