# Effective field theory for compact object evolution in binary inspirals

Irvin Martínez  & Amanda Weltman

###### Abstract

Using the effective field theory framework for extended objects we describe
the evolution of spinning compact objects in the late inspiral of the
coalescence of a binary, before the plunge and merger, by including leading
order corrections due to spin, tides, dissipation and gravitational wave
radiation. Our implementation is of particular relevance for probing the
stellar structure of compact objects with gravitational wave observations. A
spinning compact object in the effective field theory framework is described
as a spinning point particle, with its finite size effects encoded in higher
order operators in the effective action, operators which have coefficients
that encapsulates the internal structure of the star. For the inspiral regime
described by non-relativistic general relativity, post-Newtonian corrections
to each term of the action can be obtained in a diagrammatic approach,
including gravitational radiation effects. Taking into account the
aforementioned ingredients of the effective theory, we solve for the dynamics
of the inspiral of binary systems using an algorithm for point particle
simulations. We extract the gravitational wave as a function of the orbital
frequency, input that is generated numerically and then evaluated in the
analytic function of the waveform. By performing illustrative numerical
experiments of systems that the LIGO-Virgo observatories have already
detected, we show the role of the stellar structure and its coefficients in
the phase evolution of the waveform, as well as the order in which they arise
and the sensitivity required for the gravitational wave observatories to
measure them. If these coefficients are to be measured, tight constraints on
fundamental physics can be placed.

##  1 Introduction

We live in a time with access to data from the most energetic "particle"
collisions in the universe. The astonishing first detection of gravitational
waves (GWs) by the LIGO-Virgo observatories from a binary black hole (BBH)
merger [1], the multi-messenger detection from a binary neutron star (BNS)
collision [2, 3, 4], and the recent detection of the coalescence of black hole
(BH) - neutron star (NS) binaries [5], provide us with ample opportunities to
test fundamental physics in the strong regime of gravity.

These detections of binary systems through GWs have already allowed us to
better understand the properties of compact objects [6, 7]. With more
sensitivity upgrades planned for the LIGO-Virgo observatories, the upcoming
third generation Einstein Telescope [8, 9] and the space based detector, LISA
[10], the era of precision gravity is upon us, and with it the potential for
great discovery using multi-messenger astronomy to better understand the
universe on all scales. On the theoretical side, the development of improved
tools to model the dynamics of any astrophysical source to extract the
gravitational wave signature, can only aid in this discovery.

One of the key potentials with GW observations is to test fundamental physics
by probing the internal structure of the compact objects. For example, the
equation of state of matter for NSs is still unknown, and it has been
suggested that it can be constrained by matching the Love numbers of the NSs
with GW observations [11, 12]. The Love numbers are dimensionless parameters
that measure the rigidity and tidal deformation of the stellar object, and
different values for the Love numbers correspond to different equations of
state of matter [13]. On the other hand, BHs within the theory of general
relativity are well constrained by the no hair theorem, which states that a BH
can be described by only three parameters, its mass, spin and charge, which
implies the vanishing of its Love numbers [14, 15]. Any finding of deviation
of this description or discovery of a new parameter that describes the BH must
show hints of a more complete underlying gravitational theory [16, 17].

Therefore, the modelling of compact objects and their interactions must take
into account different effects, such as the spin and the stellar structure of
the various stellar types. In this work we bring the tools of effective field
theory (EFT) for extended objects [18, 19, 20] to model BHs and NSs as point
particles, with additional effects encoded as higher order corrections in the
action. Although at first it might seem counter intuitive to model such
massive objects as point particles, there is in fact, a lot of similarity in
the description of compact objects and elementary particles [21]. For
instance, the no hair theorem implies that a BH effectively behaves as a point
particle, and it has been shown that a spinning BH is described
relativistically to all orders in the multipole expansion from minimal
coupling, matching on to the effective action of a spinning point particle
[22, 23].

We start with an effective theory for spinning extended objects derived from
the coset construction [20, 24], a general technique from the EFT framework
that allows us to construct effective actions [25]. In this approach, the
coefficients appearing in the effective action are treated as free parameters
to be fixed by observations. This theory is of particular relevance because it
leads to a more natural description of spinning extended objects without the
need for redundant degrees of freedom, as is the case for the EFT for spinning
extended objects used in [26, 27].

For the stellar structure, we build on [28], where in the context of non-
compact objects, static tidal effects and dissipation were considered, using
the tools introduced in [18], and in [19], for each effect respectively. In
this work we go further to consider compact objects and a response function
for dynamical tides of NSs derived in [29] using the tools from [18, 19], and
consider relativistic effects, such as GW radiation, determined by the PN
expansion [18]. Although it has been pointed out that this is a theory for
"slowly" spinning rigid objects due to the existence of a rotational frequency
at which the theory breaks down [20], our current observations of spinning
compact objects [6] suggest that we can safely consider that most of the
astrophysical objects spin slowly, and therefore are well described by our
effective theory.

The description of the coalescence of binary systems through GW radiation is
divided into four stages: the inspiral, plunge, merger and ringdown, which is
depicted in fig. 1. The inspiral phase is described using the PN expansion
[18, 30]. The transition from the inspiral to the plunge can be modelled
analytically [31], although it is usually modelled altogether with the merger
by solving the full Einstein’s equations numerically [32], which is
computationally expensive. After merging, if the final outcome is a BH, it
enters into a ringdown phase that can be modelled analytically using
perturbation theory [33].

Furthermore, the need to extract efficient and accurate waveforms for all of
the phases of the coalescence led to the development of the effective one body
(EOB) framework [34], which is a combination of analytical methods based on
the Hamiltonian, that uses the PN expansion as an input for the inspiral, and
non-perturbative methods and numerical relativity results to model the merger
semi-analytically. Since the introduction of the EOB many additional
improvements have been done to include the stellar structure, such as tidal
effects [35, 36, 37] for NSs, and BHs horizon GW absorption or dissipative
effects [38]. The EOB is the current framework used for gravitational wave
extraction and comparison to observations. In fig. 1, we show a waveform from
an equal mass BBH extracted using the EOB method from the LIGO library [39].

Figure 1: The coalescence of an equal mass BBH system, each with mass,
$M=20\>M_{\odot}$, and initial gravitational wave frequency,
$f_{GW}=50\mathrm{Hz}$. The left figure represents the dynamics of the binary
until the collision of the BHs using our numerical code for point particles,
with the purple and cyan line the trajectories followed by each of the BHs.
The orbital trajectory once the LSO is reached is shown in red until the radii
of the holes interact. The system decays due to GW radiation by including the
leading order 2.5 PN correction in the equations of motion. The figure on the
right is the gravitational waveform extracted from the merger of the same
binary but using the EOB method from the LIGO Algorithm Library [39]. The
coalescence is divided into four different phases: inspiral, plunge, merger
and ringdown.

The properties of the compact objects play an important role mostly in the
late inspiral phase of the coalescence, and thus we restrict ourselves to this
regime where the dynamics are described by the PN approximation, before the
plunge and merger of the binary. The PN expansion encodes non-linearities of
the theory of general relativity and its predictability has been tested by
comparing PN to numerical relativity simulations. It has been found that the
PN expansion reproduces the waveform with high accuracy up to a few orbits
before the merger [40, 41], where the PN approximation becomes less accurate
as one approaches the inner most circular orbit, or last stable orbit (LSO)
[42]. The LSO is the closest distance at which a point particle has a stable
orbit around the stellar object and it will set the limit of our theory. Once
the LSO has been reached, the binary enters into the plunge phase. In figure
1, we show the dynamics of the coalescence of an equal mass BBH in a quasi-
circular orbit made with our simulations, and show where our theory breaks
down.

The PN expansion in the EFT for extended objects can be obtained by computing
the corresponding Feynman diagrams of the interaction [18]. In this approach,
relativistic corrections are described as a perturbative series in terms of
the expansion parameter, $v/c<1$, with $v$ the relative velocity of the
binary. For each $n$-PN order, the expansion of the equations of motion is of
order $v^{2n}$. Although high order PN corrections are necessary for a precise
description of the interaction, in this paper we consider only the leading
order PN correction due to gravitational wave radiation, the $2.5-$PN term
[18], and neglect the leading order conservative $1$ and $2$-PN effects which
only contribute to a shift of the orbit and do not contribute to the decay of
the binary. With the purpose of building an intuitive understanding of our EFT
to perform numerical simulations, we skip the use of Feynman diagrams and take
the Newtonian limit of the effective action as in [28], which is contained in
the leading order correction of each term in the action. The PN expansion of
the effective action is discussed in [43], where the LO expansion to each of
the terms considered in this work is obtained. The state of the art of EFT and
the PN expansion can be found in [18, 26, 27, 30, 44, 45, 46, 47, 48, 49, 50,
51, 52, 53].

Taking into account the aforementioned ingredients, and with the point
particle approximation in mind, we obtain the equations of motion for the
effective action and introduce a known methodology for simulating stellar
dynamics as point particles [54] into the framework of EFT for compact
objects. Our implementation of the coalescence of a binary is tested with
analytical results, and the stellar structure effects are internally tested by
matching the coefficients of the theory during the simulation. Solving
numerically for the point particle as well as for its additional effects,
allows us to evolve the system with high accuracy for low computational cost.
We propose a hybrid method for extracting gravitational waveforms as a
function of the orbital frequency of the binary, $\omega_{s}$, input which is
generated numerically and then evaluated in the analytical formula of the
waveform. Our method reproduces results from known methodologies that
evaluates the coalescence time, $t_{c}$, of the binary.

We perform illustrative numerical "experiments" of the inspiral of the
different detected binaries to study the imprints of the stellar structure in
the dynamics. With these, we show the role of the stellar structure and its
coefficients in the waveform, given that the addition of the stellar
properties in the equations of motion generates changes in the orbital
frequency, which can be measured as a phase shift in the waveform. We show the
order at which the stellar structure effects arise in the waveform, and the
sensitivity that the current LIGO-Virgo observatories will need to detect
them. We suggest that the matching of the coefficients due to the internal
structure of spinning NSs, $n_{\Omega}$, is the first to be matched from
observations, given that it can be the leading order effect, making a wave
shift in the waveform of order up to, $\Delta\varphi\simeq O(10)$, given our
considered spinning objects. Moreover, we perform a systematic study of the
internal structure of the objects, and show that, at least at this level of
the implementation, in order to constrain the coefficients of all the leading
order effects due to the internal structure of the compact objects, current
and future GW observatories need to measure wave shifts with an accuracy of at
least, $\Delta\varphi\simeq O(10^{-4})$ radians.

The plausibility for LIGO-Virgo observatories and the Einstein telescope to
detect such effects, is of course dependent on the frequency where the
dephasing is accumulated in the signal relative to the most sensitive band of
the detectors. At current sensitivity, if the phase shift is of order
$O(10^{0})$ radians in the most sensitive band, then the effect could be
measurable and a systematic study needs to be carried out, i.e. using Fisher
matrix estimates [55]. With upcoming sensitivity upgrades, a systematic study
needs to be done to put bounds on the measurability of such effects, and
whether all the effects due to the internal structure of the compact object
can be constrained.

In section 2, we explain the effective action which systematically includes
spin/size, finite-size and dissipative effects. Then take the Newtonian limit,
derive the equations of motion and consider leading order GW radiation effects
to the point mass term of the action. In section 3, we introduce the hybrid
method for gravitational wave extraction and the signatures due to the stellar
structure. In section 4, we implement the formalism numerically, test it with
analytical results and measure the coefficients of the theory as an internal
tests of the implementation. We then perform illustrative numerical
experiments of the detected binaries by LIGO and quantify the effects in the
gravitational waveform. In section 5, we discuss the implications of our
results.

##  2 Effective field theory for compact objects

Within the framework of EFT for extended objects [18, 19], we build an
effective action for spinning compact objects that takes into account the
stellar structure. The key idea of modeling compact objects with EFT is that
we can treat them as point particles, with the stellar structure encoded in
higher order operators. This is represented by the effective action

$\displaystyle\mathcal{S}=\int\mathrm{d}\tau\left\\{-mc^{2}+\sum_{n}c_{n}\tilde{\mathcal{O}}_{n}\right\\},$
(2.1)

with the first term describing a point particle, and the second, the sum over
all possible higher order corrections. The coefficients, $c_{n}$, are the
Wilson coefficients of the effective theory that are to be matched from
experiments or from the full known theory, and $\mathcal{O}_{n}$, higher order
operators that are allowed by the symmetries. This sum is in principle
infinite, and can be cut-off at the desired accuracy. Some of the Wilson
coefficients are related to the Love numbers, making the EFT a very good
framework to work with. We refer to the Wilson coefficients that contain Love
numbers as Love coefficients. Furthermore, not all the coefficients are Love
coefficients, i.e. they instead can encode the capacity of a BH to absorb GWs,
or the viscosity of a NS which can generate an energy loss during the
interaction. We simply refer to this coefficient as the dissipation
coefficient.

The action for spinning objects in gravity derived from the coset construction
[20], taking into account the stellar structure [28, 24], has the schematic
form

$\mathcal{S}=\mathcal{S}_{PP}+\mathcal{S}_{\Omega}+\mathcal{S}_{Q}+\mathcal{S}_{\mathcal{D}},$
(2.2)

where, $\mathcal{S}_{PP}$, is the action for a spinning point particle,
$\mathcal{S}_{\Omega}$, gravitational corrections due to spin,
$\mathcal{S}_{Q}$, dynamical tides and $\mathcal{S}_{\mathcal{D}}$,
dissipative effects. We consider each of the leading order terms for these
effects in turn, to see their overall role. Thus, the effective action reads

$\displaystyle\begin{split}\mathcal{S}=&\int\mathrm{d}\tau\left\\{-mc^{2}+\frac{I}{2}\tilde{\Omega}_{ij}\tilde{\Omega}^{ij}+n_{\Omega}\tilde{\Omega}^{i}\tilde{\Omega}^{j}\tilde{C}_{0i0j}+\tilde{Q}^{ij}\tilde{C}_{0i0j}+\tilde{\mathcal{D}}^{ij}\tilde{C}_{0i0j}+.\>.\>.\right\\},\\\
\end{split}$ (2.3)

where $\tau$ is the proper time. The indices $i,j,k,...$, denote spatial
components of the Lorentz indices and the tildes indicates quantities in the
comoving frame. This action is Lorentz invariant and is compatible with all
possible symmetries.

The first two terms of equation (2.3) describes the leading order corrections
for a relativistic spinning particle with mass, $m$, moment of inertia, $I$,
and angular velocity,
$\tilde{\Omega}_{ij}=\Lambda^{\;\;k}_{i}D_{\tau}\Lambda_{kj}$, and
$\tilde{\Omega}_{ij}=\epsilon_{ijk}\tilde{\Omega}^{k}$. The third term is a
spin correction due to gravity, with $\tilde{C}_{0i0j}$, the traceless Weyl
tensor. The fourth term accounts for dynamical quadrupolar tidal effects with,
$\tilde{Q}_{E}=\tilde{F}(\tau)\tilde{C}_{0i0j}$ and $\tilde{F}(\tau)$, a
dynamical response function whose Fourier transformation can be Taylor
expanded around small frequencies, $\omega$, such that
$\mathcal{F}(\tilde{F})\propto n_{E}+n^{\prime}_{E}\omega^{2}$.111We have
explicitly separated dissipative effects from the response function in the
action compared to [29]. [29] The coefficient, $n_{E}$, is the dimensionful
Love number for static tides [29], or simply Love coefficient, which is
related to the dimensionless Love number $k$, via $n_{E}=2k\ell^{5}/3G$ [11,
56], with $\ell$, the radius of the star. Furthermore the coefficient,
$n^{\prime}_{E}$, is the dynamical Love coefficient. These quadrupolar effects
encode the fact that the object has finite size. Finally, the last term
contains dissipative effects due to the internal structure of the stellar
object, with $\tilde{\mathcal{D}}^{ij}$ a composite operator that represents
these leading order additional degrees of freedom in a model independent way
[19]. The ellipsis represents the tower of higher order operators not taken
into account in this work.

We can ensure that our theory will be predictive by expanding our action in
small parameters $v/c$, $\ell/r$ and $\Omega/\Omega_{0}$ to the desired
accuracy. The scales, $\ell$, the radius of the object and $\Omega_{0}$, the
typical frequency, do not appear explicitly in the action but determine the
characteristic size of the dimensionful coupling constants. The tower of
higher order spin terms is under control as long as the rotational velocity is
much less than the speed of sound $c_{s}$ of the material in the star, i.e.
$\Omega\ell\ll c_{s}$ [20]. For a star, this rotational frequency would be
such that the body would undergo large nonlinear stresses and order one
distortions, breaking down the theory.

###  2.1 Stellar dynamics

We perform a Lorentz transformation from the frame embedded in the rigid body
to the lab frame. We then take the Newtonian limit of our effective action,
which is contained in the lowest order PN expansion of each of the terms.
Therefore, in the lab frame, $\tau=t$,
$C_{0i0j}\approx\partial_{i}\partial_{j}\Phi$, with $\Phi$ the Newtonian
potential, and
$\Omega^{i}=\frac{1}{2}\epsilon^{ijk}\mathcal{R}_{jl}\partial_{t}\mathcal{R}_{k}^{\;l}$,
with $\mathcal{R}(\theta^{i})$ a rotation matrix, where the $\theta^{\prime}$s
are the Euler angles describing the orientation of the rigid body [28, 24].

The action then reads,

$\displaystyle\begin{split}\mathcal{S}=&\int\mathrm{d}t\left\\{\frac{mv^{2}}{2}-m\Phi+\frac{I}{2}\Omega_{ij}\Omega^{ij}+\frac{n_{\Omega}}{2}\Omega^{i}\Omega^{j}\partial_{i}\partial_{j}\Phi\right.\\\
&\left.+\frac{n_{E}}{4}\partial^{i}\partial^{j}\Phi\partial_{i}\partial_{j}\Phi+\frac{n^{\prime}_{E}}{4}\partial^{i}\partial^{j}\dot{\Phi}\partial_{i}\partial_{j}\dot{\Phi}+\frac{1}{2}\partial_{i}\partial_{j}\Phi\mathcal{R}^{i}_{\;k}\mathcal{R}^{j}_{\;l}\tilde{\mathcal{D}}^{kl}_{E}+.\;.\;.\right\\}.\end{split}$
(2.4)

Each term in the action has a simple physical interpretation as described in
[28]. The first three terms describe a non-relativistic spinning point
particle coupled to gravity. The fourth term describes the coupling between
gravity and the ensuing quadrupole. By expanding the gravitational potential
around some background value $\bar{\Phi}$, this term can be seen as a
deformation of the inertia tensor of the form $\delta I\propto
n_{\Omega}\partial^{2}\bar{\Phi}$. The first term in the second line describes
the coupling of the induced quadrupole $\delta Q\propto
n_{E}\partial^{2}\bar{\Phi}$, and the second term its dynamical part. The last
term contains dissipative effects, which are encoded in the composite operator
$\tilde{\mathcal{D}}^{ij}_{E}$.

As we are not interested in explicitly keeping the degrees of freedom encoded
in $\tilde{\mathcal{D}}^{ij}_{E}$, it is necessary to average over them in a
systematic way using the in-in formalism [57] for this context [18, 28], which
for dissipative systems, leads to obtainment of the equations of motion from a
modified variation [58],

$\delta\mathcal{S}+i\int\mathrm{d}t\mathrm{d}t^{\prime}\delta
J_{ij}(t)\tilde{G}^{ijkl}_{R}(t-t^{\prime})J_{kl}(t^{\prime})=0,$ (2.5)

where
$J_{kl}=\frac{1}{2}\partial_{i}\partial_{j}\Phi\mathcal{R}^{i}_{\;k}\mathcal{R}^{j}_{\;l}$,
and $\tilde{G}_{R}$ is the retarded correlation function of the operators
$\tilde{\mathcal{D}}^{ij}_{E}$ [28]. Then, by considering low frequencies,
from which we assume that the degrees of freedom from
$\tilde{\mathcal{D}}^{ij}_{E}$ are near equilibrium, the time ordered two
point correlation function imply that the Fourier transform $\tilde{G}_{R}$
must be an odd, analytic function of the frequency $\omega>0$ [19]. Thus, the
retarded correlation function $\tilde{G}_{R}$ reads

$\tilde{G}^{ijkl}_{R}(\omega)\simeq\eta_{E}\omega\left(\delta^{ik}\delta^{jl}+\delta^{il}\delta^{kj}-\frac{2}{3}\delta^{ij}\delta^{kl}\right),$
(2.6)

with the coefficient for dissipative effects, $\eta_{E}\geq 0$.

Then, by combining eq. (2.6) with the modified variation eq. (2.5), and
varying it with respect to $x^{i}$ and $\theta^{i}$, the equation of motion
reads

$\displaystyle\begin{split}m\dot{v}_{i}=&-m\partial_{i}\Phi+\frac{n_{\Omega}}{2}\Omega^{k}\Omega^{j}\partial_{i}\partial_{j}\partial_{k}\Phi\\\
&+\frac{n_{E}}{2}\partial_{i}\partial_{j}\partial_{k}\Phi\partial^{j}\partial^{k}\Phi+\frac{n^{\prime}_{E}}{2}\partial_{i}\partial_{j}\partial_{k}\dot{\Phi}\partial^{j}\partial^{k}\dot{\Phi}\\\
&-\frac{\eta_{E}}{2}\partial_{i}\partial_{j}\partial_{k}\Phi(\partial^{j}\partial^{k}\dot{\Phi}+2\partial^{j}\partial_{l}\Phi\epsilon^{klm}\Omega_{m}),\end{split}$
(2.7)

for the acceleration of the compact object, and

$\displaystyle\begin{split}&\partial_{t}(I\Omega_{i}+n_{\Omega}\Omega^{j}\partial_{i}\partial_{j}\Phi)=-n_{\Omega}\epsilon_{ijk}\Omega^{k}\Omega^{l}\partial^{j}\partial_{l}\Phi\\\
&+\eta_{E}\partial^{j}\partial^{k}\Phi(3\partial_{i}\partial_{j}\Phi\Omega_{k}-2\partial_{j}\partial_{k}\Phi\Omega_{i}+\epsilon_{ikl}\partial^{l}\partial_{j}\dot{\Phi}),\end{split}$
(2.8)

for the change of rotational angular momentum with contributions from
$n_{\Omega}$ and $\eta_{E}$. These equations of motion were originally
obtained in [28], with the exception of the dynamical tides term with
coefficient $n^{\prime}_{E}$, which we have derived.

By substituting the Newtonian potential, $\Phi$, and restricting the spin,
$\Omega^{i}$, to be aligned with the angular momentum of the binary, we obtain

$\displaystyle\begin{split}m_{1}\dot{\vec{v}}=&-\frac{Gm_{1}m_{2}\hat{r}}{r^{2}}-\frac{3n_{\Omega}Gm_{2}\Omega^{2}\hat{r}}{2r^{4}}-\frac{9n_{E}G^{2}m_{2}^{2}\hat{r}}{r^{7}}\\\
&-\frac{9\eta_{E}G^{2}m_{2}^{2}}{r^{8}}(\vec{r}\times\vec{\Omega}+\vec{v}+2\hat{r}(\hat{r}\cdot\vec{v}))\\\
&-\frac{18n^{\prime}_{E}G^{2}m_{2}^{2}}{r^{9}}\left(2v^{2}\hat{r}+5(\hat{r}\cdot\vec{v})^{2}\hat{r}-\vec{v}(\hat{r}\cdot\vec{v})\right),\\\
\end{split}$ (2.9)

and

$\displaystyle\begin{split}&\left(I+\frac{Gm_{2}n_{\Omega}}{r^{3}}\right)\partial_{t}\vec{\Omega}=\frac{3Gm_{2}n_{\Omega}(\hat{r}\cdot\vec{v})\vec{\Omega}}{r^{4}}+\frac{9\eta_{E}G^{2}m_{2}^{2}}{r^{7}}\left(\hat{r}\times\vec{v}-r\vec{\Omega}\right).\end{split}$
(2.10)

These equations of motion are to be integrated numerically. We will denote the
acceleration of the star as
$\vec{a}_{\star}=\vec{a}_{PP}+\vec{a}_{\Omega}+\vec{a}_{Q}+\vec{a}_{\mathcal{D}}$,
with $a_{PP}$, the point particle or Newtonian term, $a_{\Omega}$, the
gravitational correction due to spin, $a_{Q}$, tides both static and dynamic,
and $a_{\mathcal{D}}$ dissipative effects. Although we have restricted the
spin to be aligned for simplicity, in principle any spin can be considered.
Furthermore, these equations of motion are valid for any stellar object that
can be described by an equation of state of matter, with the only difference
on the values that the coefficients will take for different stellar types.

### Slowly spinning objects

As previously mentioned, our theory is valid as long as, $\Omega\ell\ll c_{s}$
holds, which can be tested once we specify the properties of the objects
below. First, it is important to note that instead of specifying values of the
angular velocity, we will work in terms of the dimensionless spin,
$\vec{\chi}$, of the object

$\chi=\frac{cJ}{GM^{2}}=\frac{cI\Omega}{GM^{2}},$ (2.11)

where, $J=I\Omega$, is the scalar value of the angular momentum of the star.
In order for our theory to work appropriately, we need $\chi$, to be slow and
non-negligible $0<\chi\ll 1$, such that we have a significant contribution in
the simulations and not break down the theory.

#### Neutron stars

We start by introducing the inner most circular orbit distance or LSO distance
$d_{LSO}$ for a NS. For a non-spinning NS, the spacetime outside it is
described by the Schwarzschild metric, which leads to $d_{LSO}=6Gm/c^{2}$,
with $m$ the mass of the star. For the case of a spinning NS, if the orbiting
object is in prograde motion, to first order in $\chi$, the LSO is
$d_{LSO,\Omega}=6Gm/c^{2}(1-0.54433\chi)$ [59], with $\chi$ being the
dimensionless spin parameter defined above. Nevertheless, for our simulations
we use the $d_{LSO}$ for the non-spinning star as the limit of the simulation.
The moment of inertia for a NS, can be approximated as [60]

$I=0.21\frac{m\ell^{2}}{1-2Gm/\ell c^{2}},$ (2.12)

where $\ell$ is the radius of the NS. We take the speed of sound inside the
star $c_{s}\lesssim c/\sqrt{3}$ [61].

To take into account tidal effects, we use the static Love coefficient,
$n_{E}=2k\ell^{5}/3G$, with $k=[0.449154,0.36260,0.259909]$, the dimensionless
quadrupolar ($k=k_{2}$) Love number for NSs, for each polytropic index,
$n=[0.5,0.7,1]$, respectively [11, 12, 62] and the spin coefficient,
$n_{\Omega}$, in the Newtonian limit $n_{\Omega}=n_{E}$ [63]. For the
dynamical component of the tidal effects, considering a NS with
$m=1.2M_{\odot}$, $\ell=8.89$ km and polytropic index, $n=1$, the
dimensionless dynamical Love number $k^{\prime}$, related to $n^{\prime}_{E}$,
is obtained from the response function of the NS in terms of the frequency
$\omega_{f}$, and dimensionful overlap integral, $\mathcal{I}_{f}$, of the
fundamental mode of the star [29], which describes to which extent an external
field excites the mode. $\mathcal{I}_{f}$ corresponds to $\mathcal{I}_{02}$ in
[29], where the subscript indicates the fundamental mode, $l=0$, and the
quadrupolar moment, $k=k_{2}$. The response function reads [29]

$\displaystyle\begin{split}\mathcal{F}(F)&=\frac{1}{2}\frac{\ell^{5}}{G}\frac{q_{f}^{2}}{\ell^{2}(\omega_{f}^{2}-\omega^{2})/c^{2}}=\frac{1}{2}\frac{\ell^{5}}{G}\frac{q_{f}^{2}}{\ell^{2}\omega_{f}^{2}/c^{2}}\left(1+\frac{\omega^{2}}{\omega_{f}^{2}}+.\;.\;.\right),\end{split}$
(2.13)

where we have expanded over $\omega/\omega_{f}$, and $q_{f}$, is the
dimensionless overlap integral which is related to the $\mathcal{I}_{f}$
through $q_{f}^{2}=G\mathcal{I}^{2}_{f}/\ell^{3}$. The combination,
$\ell\omega_{f}/c$, is dimensionless as well.

We identify the Love coefficient

$n_{E}=\frac{1}{2}\frac{\ell^{5}}{G}\frac{q_{f}^{2}}{\ell^{2}\omega^{2}_{f}/c^{2}},$
(2.14)

and comparing to the Newtonian tidal Love number, $n_{E}=2k\ell^{5}/3G$, we
identify the dimensionless Love number $k$

$k=\frac{3}{4}\frac{q_{f}^{2}}{\ell^{2}\omega^{2}_{f}/c^{2}},$ (2.15)

from which we find agreement with $k$ obtained using the Clairaut-Radau
equation [64]. From the above expansion of the response function, we can
extract the term $n^{\prime}_{E}$

$n^{\prime}_{E}=\frac{1}{2}\frac{\ell^{7}}{Gc^{2}}\frac{q_{f}^{2}}{\ell^{4}\omega_{f}^{4}/c^{4}}=\frac{1}{2}\frac{\ell^{7}}{Gc^{2}}k^{\prime},$
(2.16)

for a limit in which $\omega/\omega_{f}\ll 1$. We have defined $k^{\prime}$ as
the dimensionless part of $n^{\prime}_{E}$ and found $k^{\prime}=1.1311$, for
this particular model of NS with polytropic index $1$.

For all the considered NS simulations, we consider the mass and radius defined
above, even when considering different equations of state. From the considered
speed of sound inside the NS, and its radius, we can safely add spin,
$\chi\leq 0.2$, without breaking down the theory. Furthermore, is worth noting
that the coefficient regarding dissipative effects in NSs is unknown and is to
be matched from hydrodynamical numerical simulations.

#### Black holes

The radius of a Schwarzschild BH is $\ell_{\bullet}=2GM/c^{2}$. The LSO
distance for a nonspinning BH is $d_{LSO}=6Gm/c^{2}$. The Love numbers for BHs
within general relativity are expected to be zero [15, 56]. The extra
parameter that can describe a black hole besides its mass and charge, encodes
the capacity for the BH to absorb GWs through its horizon. Dissipative effects
for BHs due to GW radiation in a binary decay were first taken into account in
[14], and added into the EFT framework in [19]. In the latter, the graviton
absorption cross section for nonrotating black holes was matched from [65,
66], in which the response function containing dissipative effects was
obtained, and from which the dissipative coefficient can be read off. The
coefficient for dissipative effects for nonrotating BHs reads [19]

$\eta_{E}=\frac{16}{90}\frac{G^{5}M^{6}}{c^{13}}=\frac{\ell^{6}_{\bullet}}{360Gc},$
(2.17)

with $M$ the mass of the BH.

In the spinning case, there are more coefficients to account for dissipative
effects, including the coefficient, $\eta_{E}$, for the nonrotating black hole
[15]. The explicit expression of these coefficients has been shown in [24].
Given that our current numerical implementation does not allow us to take into
account these extra coefficients, we do not perform simulations of rotating
BHs.

###  2.2 Non-relativistic general relativity

In the EFT for extended objects [18], the PN expansion for the conservative
and radiative dynamics of the binary is organized and obtained by solving tree
level Feynman diagrams. This framework has allowed us to obtain higher order
PN corrections than those derived using different methods, reproducing the
well known results up to 3-PN order [18, 44], and deriving new results for the
4-PN [44] and recently for the 5-PN [67] corrections for the non-spinning
case. Spinning compact objects were introduced into the EFT in [26] to obtain
the PN contributions from the spin-orbit and spin-spin coupling, and further
corrections were obtained in [27, 46, 47, 48, 49, 50, 53]. The state of the
art for spinning objects is 4-PN order [68]. A review on the PN expansion
using EFTs can be found in [51, 52].

Our model for spinning objects has a different construction [20, 24], for
which the leading order PN expansion has been obtained in [43], reproducing
known results from [50]. The Newtonian action, eq. (2.4), is contained in the
leading order PN correction to each effect [43]. Nevertheless, our work
focuses in the internal structure of the compact objects, and the
implementation of the conservative PN corrections is left for a future work.
The leading order non-conservative correction to the point particle due to
gravitational wave radiation, can be extracted from the literature [18, 30].

#### Post-Newtonian dynamics

From the EFT perspective, we can obtain the PN corrections to each of the
terms in the action, eq. (2.4). We make use of the PN expansion to the point
particle term [18]. The PN corrections to the point mass can be expanded in a
series as

$\vec{a}_{PN}=\vec{a}_{PP}+c^{-2}\vec{a}_{2}+c^{-4}\vec{a}_{4}+c^{-5}\vec{a}_{5}+O(c^{-6}),$
(2.18)

where, $\vec{a}_{PP}=\vec{a}_{N}$, is the point particle acceleration or
Newtonian term. The terms $\vec{a}_{2}$ and $\vec{a}_{4}$ are the first and
second order PN contributions of the conservative sector which accounts for,
i.e. the periastron shift. The leading order term that accounts for radiation
of energy and momentum through GW radiation is the 2.5 PN term, $\vec{a}_{5}$,
which we denote as, $\vec{a}_{GW}=\vec{a}_{5}$. This correction reads [18, 30]

$\displaystyle\begin{split}\vec{a}_{GW}=&\frac{4}{5}\frac{G^{2}m_{1}m_{2}}{r^{3}}\left[\left(\frac{2Gm_{1}}{r}-\frac{8Gm_{2}}{r}-v^{2}\right)\vec{v}+(\hat{r}\cdot\vec{v})\left(\frac{52Gm_{2}}{3r}-\frac{6Gm_{1}}{r}+3v^{2}\right)\hat{r}\right].\end{split}$
(2.19)

For the rest of the paper we make use of the modified acceleration,
$\vec{a}=\vec{a}_{\star}+c^{-5}\vec{a}_{GW}$, where $a_{\star}$, is the
acceleration which includes the Newtonian term and the stellar structure of
the star that we have taken into account in eq. (2.9). We neglect the PN
contribution from the conservative sector given that it does not contribute to
the decay of the binary. The GW decay due to the $2.5$ PN term is shown in
Fig. 1 for an equal mass BBH.

##  3 Waveforms and observational signatures

###  3.1 Waveform extraction

Gravitational radiation is produced at lowest order by the time varying mass
quadrupole moment. An analytical expression for the quadrupole gravitational
waveform from the inspiral of a binary can be obtained in the setting of flat
background spacetime with linearized gravity. One can solve the linearized
Einstein’s equations in the presence of the binary as a matter source and
obtain a solution for the amplitude, project it into the transverse-traceless
gauge, and expand it as a multipole expansion. The solution is the amplitude
for each polarization mode, $h_{+}$ and $h_{\times}$ [69].

Binaries whose orbital distance is large enough, such that changes in the
orbital distance due to GW radiation over several periods are small, can be
considered as binaries with fixed orbits. In this scenario, the quadrupolar,
$+$ polarized, GW amplitude for an inspiral binary in a circular fixed orbit
reads [69]

$\displaystyle
h_{+}(\omega_{s},t)=\frac{2}{r_{o}}\frac{G^{5/3}}{c^{4}}\frac{m_{1}m_{2}}{m^{1/3}}\omega_{s}^{2/3}\cos(2\omega_{s}t)=\frac{2G\mu
x}{c^{2}r_{o}}\cos(2\omega_{s}t),$ (3.1)

where $r_{o}$, is the distance from the observer to the binary, $\mu$, the
reduced mass, $m$, the total mass of the binary, and $\omega_{s}$, the orbital
frequency. We have introduced the dimensionless variable,
$x\equiv\left(Gm\omega_{s}/c^{3}\right)^{2/3}$. Furthermore, we have chosen
the inclination angle, $\iota=\pi/2$, such that, the $\times$ polarized
amplitude, $h_{\times}$, vanishes.

For the late inspiral regime, where the changes in the orbital distance and
frequency are relevant, it is necessary to take into account the orbital decay
in the waveform. This can be done by obtaining the radiated power, the amount
of energy radiated through GWs per unit time, and relating it to the change of
energy of the system, $\dot{E}_{orb}$, as $P=-\dot{E}_{orb}$. From the latter,
an expression for the change of the orbital frequency in time,
$\dot{\omega}_{s}$, is obtained in terms of $\omega_{s}$. By integrating it
over the coalescence time measured by an observer, $t_{c}=T_{c}-t$, with
$T_{c}$, the time that the binary takes to merge, one obtains the relation
between orbital frequency and the time of coalescence,

$\displaystyle\omega_{s}(t_{c})=\frac{1}{8}\left(\frac{5c^{5}}{G^{5/3}}\frac{m^{1/3}}{m_{1}m_{2}}\right)^{3/8}t_{c}^{-3/8}.$
(3.2)

Then, in eq. (3.1), one replaces the phase,
$2\omega_{s}t\rightarrow\varphi(t)$, with the accumulated orbital phase,

$\displaystyle\varphi(t)=2\int_{t}\mathrm{d}t\;\omega_{s}(t)=-2\int_{t_{c}}\mathrm{d}t_{c}\,\omega_{s}(t_{c}),$
(3.3)

where we have used $\mathrm{d}t=-\mathrm{d}t_{c}$. By substituting eq. (3.2)
into last equation, solving the integral for the time of coalescence, and
expressing it in terms of the orbital frequency, one obtains the accumulated
orbital phase,

$\displaystyle\varphi(x)=-\frac{1}{32}\frac{m^{2}}{m_{1}m_{2}}x^{-5/2}+\varphi_{0}=-\frac{1}{32}\frac{m^{2}}{m_{1}m_{2}}x^{-5/2}-\frac{2Gm\omega_{s}}{c^{3}}\log\left(\frac{\omega_{s}}{\omega_{0}}\right),$
(3.4)

where $\varphi_{0}$ has been fixed from [30], and with $\omega_{0}$, a
constant frequency that is chosen as the frequency at which the binary enters
into the detectable band.

Thus, the gravitational waveform that takes into account the GW decay, reads

$\displaystyle h_{+}(x(\omega_{s}))=\frac{2G\mu x}{c^{2}r_{o}}\cos
2\varphi(x).$ (3.5)

with $\mu$ the reduced mass. We simply denote,
$h_{+}(x(\omega_{s}))=h_{+}(\omega_{s})$. This expression is the lowest order
or $0-$PN waveform in the PN approximation. By writing
$x=((Gm/r)(r\omega_{s}/c^{3}))^{2/3}$, we find its scaling given that,
$Gm/r\sim v^{2}$, and $r\omega_{s}\sim v$. Thus, every $x$ scales as
$x=O\left(v^{2}/c^{2}\right)$, meaning that corrections in powers of $v/c$ can
be expressed in powers of $x^{1/2}$. The PN expansion of the $+$ polarized
waveform in powers of $x$, reads [30]

$\displaystyle\begin{split}h_{+}(x)=\frac{2G\mu
x}{c^{2}r_{0}}\left(H^{0}_{+}+x^{1/2}H^{(1/2)}_{+}+xH^{(1)}_{+}+x^{3/2}H^{(3/2)}_{+}+x^{2}H^{(2)}_{+}+O(c^{-5})\right),\end{split}$
(3.6)

with $H^{0}_{+}=\cos 2\varphi$. The last expression shows the program to
extract high order PN waveforms that matches observations. Nevertheless, for
the purpose of this paper we consider only the amplitude in eq. (3.5).

###  3.2 Comparison to known results

Working with the amplitude formula as a function of the orbital frequency,
$\omega_{s}$, we need to generate the input numerically by solving the
equations of motion as shown in the numerical section. This is in contrast to
the amplitude as a function of the time of coalescence, $t_{c}$ measured by
the observer, which can be obtained via eq. (3.2) in eq. (3.5). In this case,
the waveform can be simply extracted if the coalescence time is known. The
analytical formula for the time of coalescence for two point masses orbiting
each other in a circular orbit, a good approximation for a binary black hole,
reads [70]

$T_{c}(d)=\frac{5}{256}\frac{d^{4}c^{5}}{G^{3}m_{1}m_{2}(m_{1}+m_{2})},$ (3.7)

with $d$, the orbital distance of the binary. Nevertheless, the stellar
structure changes the merger time, and deriving analytical expressions becomes
a challenge. Thus, to measure the imprints of the stellar structure in the
late inspiral of the binary, we extract the waveform as a function of the
orbital frequency, $\omega_{s}$.

Figure 2: Comparison of the leading order gravitational waveform as a function
of the time of coalescence, $t_{c}$, and of the orbital frequency,
$\omega_{s}$. Using the same equal mass BBH as in figure 1, we extract the
orbital frequency and time of coalescence to obtain the waveforms. The cyan
line is the waveform as a function of $\omega_{s}$, input generated with our
numerical simulations. The magenta line is the waveform as a function of the
$t_{c}$, which is obtained from using eq. (3.7) as the only input.

To compare our methodology with known results, we show the extraction of the
orbital frequency from the decay of a binary from the simulations. Consider a
binary with masses $m_{1}$ and $m_{2}$, and relative distance,
$\vec{r}=\vec{r}_{2}-\vec{r}_{1}$, in some inertial frame. Then, considering
only Newtonian gravity, the relative motion,
$\vec{a}=\vec{a}_{2}-\vec{a}_{1}$, reads

$\displaystyle\vec{a}_{N}=-G\frac{m}{r^{2}}\hat{r}.$ (3.8)

with $m=m_{1}+m_{2}$. The Newtonian orbital frequency, $\omega_{N}$, is simply
obtained from Kepler’s law, by equating the centripetal to the Newtonian
relative acceleration. The centripetal acceleration, $a_{c}=v^{2}/r$, is a
radial force that has the opposite direction to the gravitational radial
acceleration. Furthermore, the relative velocity, $v$, is related to the
orbital frequency as $v=\omega r$. Thus, the Newtonian orbital frequency is
simply given by

$\displaystyle\omega_{N}=\sqrt{\frac{Gm}{r^{3}}}.$ (3.9)

Using the above reasoning, we obtain the orbital frequency with the GW effects
included. This can be done by equating the total radial acceleration,
including the radial term in (2.19), to the orbital frequency, $\omega_{s}$.
Thus, the orbital frequency taking into account gravitational radiation, reads

$\displaystyle\begin{split}\omega_{s}=\omega_{N+GW}=&\sqrt{\frac{Gm}{r^{3}}-\frac{4}{5}\frac{G^{2}m_{1}m_{2}}{r^{4}}(\hat{r}\cdot\vec{v})\left(\frac{52Gm}{3r}-\frac{6Gm}{r}+6v^{2}\right)}.\end{split}$
(3.10)

In figure 2, we show the comparison of our method for GW extraction dependent
on the orbital frequency using eq. (3.5) and eq. (3.10), to leading order in
the amplitude as a function of the coalescence time, $t_{c}$ [69]. Given the
initial distance, $d=10{\ell}_{\bullet}$, of an equal mass binary,
$m_{\bullet}=20M_{\odot}$, we can obtain the coalescence time and extract the
waveform, $h_{+}(t_{c})$ [69]. To extract the waveform as a function of the
orbital frequency, $h_{+}(\omega_{s})$, we generate eq. (3.10) numerically
with our code described below, and then evaluate it in eq. (3.5). Both
waveforms match onto most of the inspiral, with some discrepancies near the
coalescence time, which might be due to cutting off the PN corrections. Recall
that we are only using the leading order GW radiation correction, or 2.5 PN
term. With these results, we can safely proceed to study, to leading order,
the effects due to the stellar structure.

###  3.3 Observational signatures

There are observable consequences due to the internal properties of the star
that can be quantified. The effective potential changes the acceleration of
each object, making the frequency and amplitude of the GW shift. To obtain the
frequency shift, we include the stellar structure and derive the orbital
frequency in the same way as for the GW decay. We only show the static tidal
effects for simplicity, but all other effects are incorporated in the same
manner. The contribution from static tidal effects is purely radial, such that
the orbital frequency now reads

$\displaystyle\begin{split}\omega_{N+GW+T}=&\sqrt{\frac{Gm}{r^{3}}-\frac{4}{5}\frac{G^{2}m_{1}m_{2}}{r^{4}}(\hat{r}\cdot\vec{v})\left(\frac{52Gm}{3r}-\frac{6Gm}{r}+6v^{2}\right)+\frac{9n_{E}G^{2}(m_{1}^{2}+m_{2}^{2})}{r^{8}}},\end{split}$
(3.11)

where the subscript, $T$, in this example, is referring to static tidal
effects.

To quantify the phase shift, $\Delta\varphi$, generated in the waveform (3.5),
we compare simulations of binaries including different effects, but the same
initial conditions. For instance, consider the orbital frequency from two
binaries, one with GW radiation only, $\omega_{N+GW}$, and the other with GW
radiation and static tidal effects, $\omega_{N+GW+T}$. From equation (3.5), we
find the change in the phase shift due to the additional effects,

$\displaystyle\Delta\varphi=(\varphi_{N+GW+T}-\varphi_{N+GW}).$ (3.12)

The phase shift is to be measured at each time step of the simulation. In
general, the phase shift due to the contribution of an additional effect is
measured as

$\Delta\varphi=(\varphi_{b}-\varphi_{a}),$ (3.13)

where the $a$ and $b$ refer to different effects in the simulations.

##  4 Numerical simulations

We implement a 4th order Hermite integrator [54, 71]. We briefly review the
basics to evolve a system of point particles numerically, in order to solve
for the dynamics of the compact objects in a binary system.

###  4.1 Point particle simulations

#### Position, velocity and acceleration

In the point particle approximation, the predicted values of the position,
velocity and acceleration for the next step are obtained from

$\displaystyle\begin{split}&\vec{r}_{i+1}=\vec{r}_{i}+\vec{v}_{i}\Delta
t+\frac{1}{2!}\vec{a}_{i}\Delta t^{2}+\frac{1}{3!}\vec{j}_{i}\Delta
t^{3}+.\;.\;.\,,\\\ &\vec{v}_{i+1}=\vec{v}_{i}+\vec{a}_{i}\Delta
t+\frac{1}{2!}\vec{j}_{i}\Delta t^{2}+.\;.\;.\,,\\\
&\vec{a}_{i+1}=\vec{a}_{i}+\vec{j}_{i}\Delta t+.\;.\;.\,,\end{split}$ (4.1)

where $\vec{j}=\dot{\vec{a}}$. The force calculation is computed as

$\displaystyle\begin{split}&\vec{a}_{i}=\frac{Gm_{2}\hat{r}}{r^{2}}+.\;.\;.\,,\\\
&\vec{j}_{i}=\frac{Gm_{2}}{r^{3}}\\{\vec{v}-3(\vec{v}\cdot\hat{r})\hat{r}\\}+.\;.\;.\,,\end{split}$
(4.2)

with the position and velocity from the previous step; the ellipses denote all
other effects we have taken into account: GW radiation, spin size corrections,
dynamical tides and dissipation. The force calculation, $\vec{a}_{i+1}$ and
$\vec{j}_{i+1}$, is done using the predicted values, $\vec{r}_{i+1}$ and
$\vec{v}_{i+1}$, to get the corrected position, velocity and acceleration

$\displaystyle\begin{split}&\vec{r}_{i+1,c}=\vec{r}_{i}+\frac{1}{2}(\vec{v}_{i}+\vec{v}_{i+1})\Delta
t+\frac{1}{12}(a_{i}-\vec{a}_{i+1})\Delta t^{2},\\\
&\vec{v}_{i+1,c}=\vec{v}_{i}+\frac{1}{2}(\vec{a}_{i}+\vec{a}_{i+1})\Delta
t+\frac{1}{12}(\vec{j}_{i}-\vec{j}_{i+1})\Delta t^{2},\\\
&\vec{a}_{i+1,c}=\vec{a}_{i}+\frac{1}{2}(\vec{j}_{i}+\vec{j}_{i+1})\Delta
t,\end{split}$ (4.3)

up to the jerk correction. With this very simple framework we are capable of
numerically evolving the equations of motion of our theory, obtaining the
desired accuracy by using a specific time-step $\Delta t$.

#### Spin evolution

We solve for the angular velocity, $\vec{\Omega}$, in a similar fashion as
above. We define $\vec{b}=\partial_{t}\vec{\Omega}$ and correct the angular
velocity as follows

$\begin{split}\vec{\Omega}_{i+1,c}=\vec{\Omega}_{i}+\frac{1}{2}(\vec{b}_{i}+\vec{b}_{i+1})\Delta
t+\frac{1}{12}(\dot{\vec{b}}_{i}-\dot{\vec{b}}_{i+1})\Delta t^{2}.\end{split}$
(4.4)

###  4.2 Numerical tests and the matching of coefficients

Figure 3: Simulation of the coalescence of a BH-NS binary and the measurement
of the dimensionless static tidal number $k$ of the NS during the interaction.
The binary is composed of a BH of mass $m_{\bullet}=5M_{\odot}$, and a NS of
mass $m_{\star}=1.2\>M_{\odot}$ and radius $\ell_{\star}=8.89$ km. The
simulation has an initial orbital distance of $d=8\>\ell_{\bullet}$ and it is
stopped at the LSO distance of the BH. The figure on the left illustrates the
trajectories followed by the binary until the LSO orbit is reached. The purple
and cyan line represents the trajectory followed by the NS and the BH
respectively. On the right figure the coefficient $k$ is measured at each
fixed time-step in the simulation. The theoretical value of the coefficient
for the NS used in this case is $k=2.59909$, which can be measured with high
accuracy with an appropriate time step, $dt$. For this simulation, $dt=0.1$ in
code units.

Many Newtonian and PN codes have been implemented to do numerical experiments.
Testing the Newtonian limit is straightforward. To test the PN corrections we
use the GW time decay eq. (3.7). In figure 2, we show that our implementation
reproduces the leading order waveform as a function of the coalescence time
$t_{c}$.

To test the rest of the effects we can match the coefficients of the theory
from the corrected position, velocity and acceleration at each time-step,
which serves as an internal test of the implementation and shows the accuracy
of the simulation. From an EFT perspective, we can follow the energy hierarchy
to measure different coefficients. Consider for instance a BH-NS binary
without any spin, and neglect for a moment the GW radiation effects. The
dynamics of the NS will change due to the static tidal effects, and the
matching of the dimensionless coefficient reads

$k=\frac{1}{6}\frac{m_{1}}{m_{2}}\bigg{(}\frac{r}{\ell}\bigg{)}^{5}\left(\frac{a_{\star,j}}{a_{N,j}}-1\right),$
(4.5)

where $a_{N,j}$ is one of the spatial components of the Newtonian term
($j=\\{x,y,z\\}$) without any other effect, and $a_{\star,j}=a_{N,j}+a_{T,j}$
is the total acceleration of the star, with $a_{T,j}$ being the acceleration
due to static tides. Thus, if we solve our system in the $x-y$ plane with
these coordinates, we must match coefficients for each spatial component. Then
we add GW radiation effects and measure the coefficient as

$k=\frac{1}{6}\frac{m_{1}}{m_{2}}\bigg{(}\frac{r}{\ell}\bigg{)}^{5}\left(\frac{a_{\star,j}}{a_{N,j}}-\frac{a_{GW,j}}{a_{N,j}}-1\right).$
(4.6)

Figure 3 shows the time evolution measurement of the dimensionless static Love
number, $k$, from a BH-NS coalescence. The matching of coefficients from
static tidal effects is the simplest, but we can add all other effects and go
order by order in the energy hierarchy to measure the coefficients and test
the implementation.

Although it might seem trivial to match the coefficients within the
simulation, we suggest that in a similar manner, the matching of coefficients
from simulations can be done, to compare, i.e. our state of the art
hydrodynamics codes to our EFT for compact objects, a first step towards
matching the coefficients from GW observations.

###  4.3 Numerical experiments and observations

The LIGO-Virgo observatories detect gravitational waves from distant
astrophysical sources in the frequency range, $f_{GW}\sim[10,10000]$ Hz [72],
with $f_{GW}=\omega_{GW}/2\pi=\omega_{s}/\pi$ as in [69], and $\omega_{s}$ the
orbital frequency. Thus, for the purpose of waveform extraction and
measurement of the different effects, we run illustrative simulations of the
coalescence of various systems with initial and final frequency inside the
LIGO-Virgo band. We systematically add the stellar structure effects and
extract the lowest order gravitational waveform for different systems, (3.5).
We choose an arbitrary distance to the system, $r_{o}$, that makes the system
detectable in the LIGO-Virgo band.

Each binary system is set in a circular orbit, with initial conditions
generated by a simulation with twice the distance of the shown examples. This
is done with the purpose of avoiding numerical errors that can be generated in
the first orbits of the system.

#### Neutron star - black hole interactions

In this interaction, all the effects we have considered play a role in the
dynamics: static and dynamical tides in the NS and dissipation in the BH. We
divide the BH-NS simulations into two sets, one set to show the role of each
effect due to the internal structure of the objects, and the other set to
compare the role of different static Love coefficients by considering
different equations of state of matter for the NS. In each simulation, we have
the same initial conditions of a NS-BH binary, with the mass of the BH,
$m_{\bullet}=5\>M_{\odot}$, and the mass and radius of the NS,
$m_{\star}=1.2\>M_{\odot}$ and $\ell_{\star}=8.89$ km. The simulations start
at an initial distance of $d=25\>\ell_{\bullet}$ and initial frequency of
$f_{GW}=40$ Hz. The simulations are finished when the NS reaches the $d_{LSO}$
of the BH. Furthermore, we have included slow spin to the NS,
$\Omega=[0.1,0.2]$, aligned to the orbital momentum of the binary.

Figure 4: Gravitational waveforms and phase shifts from the late inspiral of a
BH-NS coalescence with $m_{\bullet}=5M_{\odot}$, $m_{\star}=1.2\>M_{\odot}$
and $\ell_{\star}=8.89$ km, at an initial distance of $d=12\>\ell_{\bullet}$.
The time to the LSO is $T_{c}\approx 10.096$ sec. We show the very last part
of the simulation. On the top figure, each waveform color is extracted from
simulations with different effects. The green waveform is only GWs, the
magenta is GWs + dynamical tides (which includes the static part), and the
yellow is GWs + dissipation. On the plot below, the phase shift is measured.
The cyan line is the wave shift generated by including static tidal effects to
the GWs, the magenta line by taking into account static and dynamical tides,
the yellow line by dissipation, and the green line is the phase shift that
dynamical tides generate with respect to the static part.

First we consider the BH-NS system without any spin. In Fig. 4, we compare the
phase evolution in the waveforms and quantify the phase shifts from the
binaries with different effects: GWs, GWs $+$ static tides, GWs $+$ dynamical
tides, and GWs $+$ dissipation. The effects due to the internal structure are
small, as one can not distinguish a difference by eye in the waveforms of the
upper plot of fig 4. In the lower plot of fig. 4, we measure the difference of
the accumulated orbital phase, $\Delta\varphi$, due to the different effects.
As expected from eq. (2.7), the major contribution comes from the static tides
term, having a final phase shift, $\Delta\varphi_{ST}\simeq 0.2$, in our
simulations, which is depicted cyan color. The addition of dynamical tides
contributes, $\Delta\varphi_{DT}\simeq 0.009$, with respect to the static
tides, which is shown by the green line. Finally, the effects of dissipation
from the BH contributes the least, shifting the waveform by
$\Delta\varphi_{Diss}\simeq 0.0018$.

Then we add the dimensionless spin, $\chi=[0.1,0.2]$, to the NS in the same
BH-NS system. We refer to the spin, $\Omega_{1}$ and $\Omega_{2}$, for the
dimensionless spin respectively, $\chi=[0.1,0.2]$, that is obtained form the
relation (2.11). On the top figure of fig. 5, we show the extracted
gravitational waveform from simulations containing the different effects. The
green waveform is GWs + dynamical tides (DT), the cyan is GWs + DT +
$\Omega_{1}$ , and the magenta is GWs + DT + $\Omega_{2}$. On the plot below,
the phase shift is measured. The green line is the wave shift generated by DT,
the cyan line is the wave shift generated by the spin $\Omega_{1}$, while the
magenta line is the wave shift due to $\Omega_{2}$. The wave shift generated
by dynamical tides is the same as before. The spin, $\Omega_{1}$, generates a
wave shift of $\Delta\varphi\simeq 0.45$, while $\Omega_{2}$, generates a wave
shift of $\Delta\varphi\simeq 1.8$. Therefore, the corrections due to spin-
size effects, can be the leading order contribution from the stellar
structure. The spin evolution is negligible, having the same final value as
the initial, which is expected from the aligned spin case.

Figure 5: Gravitational waveforms and phase shifts from the late inspiral of
the same BH-NS interaction described in figure 4, but with dimensionless spin
$\chi=[0.1,0.2]$ on the NS. On the top figure, each waveform color is
extracted from simulations with different effects. The green waveform is GWs +
DT, the cyan is GWs + DT + $\Omega_{1}$ , and the magenta is GWs + DT +
$\Omega_{2}$. On the plot below, the phase shift is measured. The green line
is the wave shift generated by dynamical tides, cyan line is the wave shift
generated by the spin $\Omega_{1}$, while the magenta line is the wave shift
by $\Omega_{2}$.

Finally, by considering static effects, as well as spin-size effects, in
figure 6, we show the role of different Love coefficients in the waveform by
measuring the phase shifts from the late inspiral of the same BH-NS
interaction described in figure 4 and 5, but with different dimensionless Love
numbers, with corresponding polytropic index, $n$, for the NS [12]. In all the
figures, the line color refers to different polytropic indices for the NSs.
The cyan line for $n=0.5$, the magenta for $n=0.7$ and the yellow for $n=1$.
On the top figure we show the different phase shift generated by considering
only static tides given different equations of state, without any spin, which
generates a phase shift, $\Delta\varphi\simeq[0.346,0.278,0.200]$ radians, for
each $q=[0.5,0.7,1]$, respectively. The middle plot is the phase shift with
respect to the top figure due to the inclusion of spin/size effects, for which
we have added spin $\Omega_{1}$. On the plot below we show the phase
difference from the top plot and the inclusion of the spin $\Omega_{2}$. The
spin $\Omega_{1}$ generates a wave shift of
$\Delta\varphi\simeq[0.808,0.649,0.458]$ radians, while $\Omega_{2}$ generates
a wave shift of $\Delta\varphi\simeq[3.136,2.555,1.851]$ radians, for each
$n=[0.5,0.7,1]$. As expected, our results shows that the smaller the value of
$n$, the more compact the star, thus the stronger the signatures in the
signal.

Figure 6: Phase shifts from the late inspiral of the same BH-NS interaction
described in figure 4 and 5, but with different dimensionless Love numbers,
with corresponding polytropic index $n=[0.5,0.7,1]$. In all the figures, the
line color refers to different polytropic indices for the NSs. On the top
figure we show the different phase shift generated by static tides given
different equations of state. The middle plot is the phase shift with respect
to the top figure due to the inclusion of spin/size effects, for which we have
added spin $\Omega_{1}$. On the plot below we show the phase difference from
the top plot and the inclusion of the spin $\Omega_{2}$.

#### Binary black hole coalescence

Figure 7: The phase shift of the gravitational waveform from the late inspiral
of three different BBHs that contains dissipative effects. We have performed
three different simulations, each of them starting with initial $f_{gw}=10$
Hz, but differing in the mass ratio, $q=m_{2}/m_{1}$, of the binary. The cyan
line is for an equal mass case, $q=1$, with $m_{1}=20M_{\odot}$. The magenta
line is for the unequal mass case, $q=1/2$, and the green line for $q=1/4$,
both with $m_{1}=20M_{\odot}$.

We have performed simulations for three different non-spinning BBHs to study
the effects of dissipation in the waveform, which differ in the mass ratio of
the binary, $q=m_{2}/m_{1}$. All of the simulations are set with an initial
gravitational wave frequency, $f_{gw}=10$ Hz. The cyan line is for an equal
mass case, $q=1$, with $m_{1}=20M_{\odot}$. The magenta line is for the
unequal mass case, $q=1/2$, and the green line for $q=1/4$, both with
$m_{1}=20M_{\odot}$. The phase shift due to dissipative effects is,
$\Delta\varphi\simeq[0.0009,0.0013,0.0018]$, for each $q=[1,1/2,1/4]$
respectively. Our results shows that the effect of dissipation can increases
for more unequal cases, but this also depends on the mass ratios taken into
account. For instance, a simulation of the first detected binary black hole
merger, GW150914, shows that dissipative effects makes a phase shift of
$\Delta\varphi\simeq 0.0004$, which is the smallest phase shift of all.

#### Neutron star binaries

Figure 8: Phase shifts from simulations of a NS-NS coalescence with an initial
frequency of $f_{gw}$ . We add GWs, DT and then spin to the NSs, first only to
one member and then to both. The magenta color is the phase shift due to
static tides, the cyan for dynamical tides, and the green the difference
between both. The blue line is the phase shift due to one spinning star with
$\chi=0.2$, and the orange line for the case in which both stars are spinning
with same magnitude and direction (aligned spins).

We set an equal mass NS-NS binary with initial frequency of $f_{gw}=40$ Hz. In
order to illustrate different scenarios, we run four different simulations
that includes GWs, dynamical tides and aligned spin, $\chi=0.2$, to the NSs,
first only to one member and then to both. In figure 8, we show that for tidal
effects, the contribution to the wave shift due to the static term is of
$\Delta\varphi\simeq 1$ radians, while the dynamic part contributes to the
latter by $\Delta\varphi\simeq 0.1$ radians. By comparing to the wave shift of
the BH-NS case, we find that tidal effects in a NS-NS binary have a stronger
effect on modifying the waveform. The role of spin-size effects is also
enhanced in the NS-NS binary case, for which the waveform is shifted by,
$\Delta\varphi\simeq 8$ radians, with only one spinning star, and
$\Delta\varphi\simeq 15$ radians, from both spinning NS with same $\Omega_{2}$
aligned to the orbital momentum of the binary. Each of both spinning NS
contributes the same phase shift, and just like in the BH-NS system, the spin
evolution is negligible, having the same final value as initial.

##  5 Discussion

In this work we have implemented the numerical evolution of the late inspiral
of compact objects within the EFT framework using a high accuracy 4th order
Hermite integrator for point particle simulations. We have successfully
included in the dynamics, leading order spin corrections due to gravity,
dynamical tidal and dissipative effects, and the 2.5 PN corrections due to GW
radiation. We have tested our implementation with analytical results, and
matched the coefficients of the effective theory at each time step from the
numerical simulations as an internal test. We have extracted the leading order
gravitational waveform as a function of the orbital frequency of the binary,
and showed the role of the coefficients by studying the overall phase
evolution of the waveform due to the different effects.

Although our current implementation is inaccurate for reproducing high order
PN waveforms, what can be learned from our numerical experiments is the order
at which the effects of the stellar structure arise in the gravitational
waveform. We have found that in the BH-NS case, the coefficients due to static
tides can be measured if LIGO-Virgo observatories, and future detectors, are
sensitive enough to measure a shift in the waveform of, $\Delta\varphi\simeq
O(10^{-1})$ radians, which will allow us to constrain between different Love
coefficients. Furthermore, the stellar structure of spinning NS plays a role
on shifting the waveform at order, $\Delta\varphi\simeq O(1)$ radians, which
suggest that the spin-size effects, with coefficient $n_{\Omega}$, can be the
leading order effect due to the stellar structure, depending on the value of
the spin. On the other hand, to detect dynamical tides and better constrain
the EOS of matter, it is necessary to measure, $\Delta\varphi\simeq
O(10^{-3})$ radians, as well as for dissipative effects in the NS-BH case.

In the case of the NS-NS binary, the wave shift due to tidal effects is
enhanced, with the static tidal effects generating a shift on the wave of,
$\Delta\varphi\simeq O(1)$ radians, while the dynamical tides,
$\Delta\varphi\simeq O(10^{-1})$ radians, from the contribution of both stars.
Just as in the case of the BH-NS, the gravitational spin correction is leading
order, and could generate a wave shift of $\Delta\varphi\simeq O(10)$ radians,
being the first coefficient to match from observations. Nevertheless, the
dynamics of nonaligned spins with different magnitudes may complicate this
task, as well as binaries with members described by different equations of
state. Finally, from the BBH simulations including dissipation, we have shown
that their effects generate a shift in the waveform at order,
$\Delta\varphi\simeq O(10^{-3})$ radians.

Thus, we conclude that the first coefficient to match, given the order at
which it contributes and which requires the least sensitivity, is the
coefficient, $n_{\Omega}$, from the spin-size effects of spinning NSs, both in
BH-NS and NS-NS binaries. Then the static tidal coefficient in BH-NS binaries,
and the dynamical tides coefficient in NS-NS binaries. Finally dynamical tides
in BH-NS interactions, and dissipative effects in both BH-NS and BBH, which
play a role roughly at the same order, and which requires the most sensitivity
for the observatories to measure. We argue that the sensitivity required for
the LIGO-Virgo and future GWs observatories to constrain the stellar structure
of the compact objects with high precision, needs to measure wave shifts with
an accuracy of at least, $\Delta\varphi\simeq O(10^{-4})$ radians.

This is the first step towards constructing realistic stellar models and
extracting accurate GWs within the EFT framework for a wide class of
astrophysical objects in the region in which the dynamics can be described by
the post-Newtonian expansion. In this sense, this framework can allows us to
create a rich template bank of waveforms for different classes of systems,
such as the ones described in this work, and many others that can be detected
by ground and space based GW detectors.

## Acknowledgments

We are extremely grateful to R. Penco, W. Goldberger, T. Hinderer, J.
Steinhoff, L. Heisenberg, J. Samsing, M. Levi, R. Yarza, and H. S. Chia for
the many enlightening conversations. I.M. is particularly grateful to E.
Ramirez-Ruiz, S. Rosswog and S. Nissanke for their support in early stages of
the work. We gratefully acknowledge support from the University of Cape Town
Vice Chancellor’s Future Leaders 2030 Awards programme which has generously
funded this research and support from the South African Research Chairs
Initiative of the Department of Science and Technology and the NRF. This
research was supported in part by the National Science Foundation under Grant
No. NSF PHY-1748958.

## References

  * [1] B. Abbott et al., “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett., vol. 116, no. 6, p. 061102, 2016\.
  * [2] B. Abbott et al., “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett., vol. 119, no. 16, p. 161101, 2017.
  * [3] B. P. Abbott et al., “Multi-messenger Observations of a Binary Neutron Star Merger,” Astrophys. J. Lett., vol. 848, no. 2, p. L12, 2017.
  * [4] B. P. Abbott et al., “Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A,” Astrophys. J. Lett., vol. 848, no. 2, p. L13, 2017.
  * [5] R. Abbott et al., “Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences,” Astrophys. J. Lett., vol. 915, no. 1, p. L5, 2021.
  * [6] B. P. Abbott et al., “GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs,” Phys. Rev. X, vol. 9, no. 3, p. 031040, 2019.
  * [7] R. Abbott et al., “Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog,” Astrophys. J. Lett., vol. 913, no. 1, p. L7, 2021.
  * [8] M. Punturo et al., “The Einstein Telescope: A third-generation gravitational wave observatory,” Class. Quant. Grav., vol. 27, p. 194002, 2010.
  * [9] M. Maggiore et al., “Science Case for the Einstein Telescope,” JCAP, vol. 03, p. 050, 2020.
  * [10] E. Barausse et al., “Prospects for Fundamental Physics with LISA,” Gen. Rel. Grav., vol. 52, no. 8, p. 81, 2020.
  * [11] E. E. Flanagan and T. Hinderer, “Constraining neutron star tidal Love numbers with gravitational wave detectors,” Phys. Rev. D, vol. 77, p. 021502, 2008\.
  * [12] T. Hinderer, “Tidal Love numbers of neutron stars,” Astrophys. J., vol. 677, pp. 1216–1220, 2008.
  * [13] K. Yagi and N. Yunes, “I-Love-Q Relations in Neutron Stars and their Applications to Astrophysics, Gravitational Waves and Fundamental Physics,” Phys. Rev. D, vol. 88, no. 2, p. 023009, 2013.
  * [14] E. Poisson, “Absorption of mass and angular momentum by a black hole: Time-domain formalisms for gravitational perturbations, and the small-hole / slow-motion approximation,” Phys. Rev. D, vol. 70, p. 084044, 2004.
  * [15] H. S. Chia, “Tidal deformation and dissipation of rotating black holes,” Phys. Rev. D, vol. 104, no. 2, p. 024013, 2021.
  * [16] L. Heisenberg, “A systematic approach to generalisations of General Relativity and their cosmological implications,” Phys. Rept., vol. 796, pp. 1–113, 2019.
  * [17] S. B. Giddings, S. Koren, and G. Treviño, “Exploring strong-field deviations from general relativity via gravitational waves,” Phys. Rev. D, vol. 100, no. 4, p. 044005, 2019.
  * [18] W. D. Goldberger and I. Z. Rothstein, “An Effective field theory of gravity for extended objects,” Phys. Rev. D, vol. 73, p. 104029, 2006.
  * [19] W. D. Goldberger and I. Z. Rothstein, “Dissipative effects in the worldline approach to black hole dynamics,” Phys. Rev. D, vol. 73, p. 104030, 2006\.
  * [20] L. V. Delacrétaz, S. Endlich, A. Monin, R. Penco, and F. Riva, “(Re-)Inventing the Relativistic Wheel: Gravity, Cosets, and Spinning Objects,” JHEP, vol. 11, p. 008, 2014.
  * [21] N. Arkani-Hamed, Y.-t. Huang, and D. O’Connell, “Kerr black holes as elementary particles,” JHEP, vol. 01, p. 046, 2020.
  * [22] A. Guevara, A. Ochirov, and J. Vines, “Scattering of Spinning Black Holes from Exponentiated Soft Factors,” JHEP, vol. 09, p. 056, 2019.
  * [23] M.-Z. Chung, Y.-T. Huang, J.-W. Kim, and S. Lee, “The simplest massive S-matrix: from minimal coupling to Black Holes,” JHEP, vol. 04, p. 156, 2019.
  * [24] I. Martínez, “Effective actions for compact objects in an effective field theory of gravity,” 11 2021.  arXiv:2111.09070.
  * [25] S. R. Coleman, J. Wess, and B. Zumino, “Structure of phenomenological Lagrangians. 1.,” Phys. Rev., vol. 177, pp. 2239–2247, 1969.
  * [26] R. A. Porto, “Post-Newtonian corrections to the motion of spinning bodies in NRGR,” Phys. Rev. D, vol. 73, p. 104031, 2006.
  * [27] M. Levi, “Next to Leading Order gravitational Spin1-Spin2 coupling with Kaluza-Klein reduction,” Phys. Rev. D, vol. 82, p. 064029, 2010.
  * [28] S. Endlich and R. Penco, “Effective field theory approach to tidal dynamics of spinning astrophysical systems,” Phys. Rev. D, vol. 93, no. 6, p. 064021, 2016.
  * [29] S. Chakrabarti, T. Delsate, and J. Steinhoff, “Effective action and linear response of compact objects in Newtonian gravity,” Phys. Rev. D, vol. 88, p. 084038, 2013.
  * [30] L. Blanchet, “Gravitational radiation from post-Newtonian sources and inspiralling compact binaries,” Living Rev. Rel., vol. 9, p. 4, 2006.
  * [31] A. Buonanno and T. Damour, “Transition from inspiral to plunge in binary black hole coalescences,” Phys. Rev. D, vol. 62, p. 064015, 2000.
  * [32] F. Pretorius, “Evolution of binary black hole spacetimes,” Phys. Rev. Lett., vol. 95, p. 121101, 2005.
  * [33] E. W. Leaver, “An Analytic representation for the quasi normal modes of Kerr black holes,” Proc. Roy. Soc. Lond. A, vol. 402, pp. 285–298, 1985.
  * [34] A. Buonanno and T. Damour, “Effective one-body approach to general relativistic two-body dynamics,” Phys. Rev. D, vol. 59, p. 084006, 1999\.
  * [35] T. Damour and A. Nagar, “Effective One Body description of tidal effects in inspiralling compact binaries,” Phys. Rev. D, vol. 81, p. 084016, 2010\.
  * [36] T. Hinderer et al., “Effects of neutron-star dynamic tides on gravitational waveforms within the effective-one-body approach,” Phys. Rev. Lett., vol. 116, no. 18, p. 181101, 2016.
  * [37] J. Steinhoff, T. Hinderer, T. Dietrich, and F. Foucart, “Spin effects on neutron star fundamental-mode dynamical tides: Phenomenology and comparison to numerical simulations,” Phys. Rev. Res., vol. 3, no. 3, p. 033129, 2021\.
  * [38] S. Bernuzzi, A. Nagar, and A. Zenginoglu, “Horizon-absorption effects in coalescing black-hole binaries: An effective-one-body study of the non-spinning case,” Phys. Rev. D, vol. 86, p. 104038, 2012.
  * [39] LIGO Scientific Collaboration, “LIGO Algorithm Library - LALSuite.” free software (GPL), 2018.
  * [40] J. G. Baker, J. R. van Meter, S. T. McWilliams, J. Centrella, and B. J. Kelly, “Consistency of post-Newtonian waveforms with numerical relativity,” Phys. Rev. Lett., vol. 99, p. 181101, 2007.
  * [41] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, “High-accuracy comparison of numerical relativity simulations with post-Newtonian expansions,” Phys. Rev. D, vol. 76, p. 124038, 2007.
  * [42] A. Buonanno, B. R. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, “Comparison of post-newtonian templates for compact binary inspiral signals in gravitational-wave detectors,” Phys. Rev. D, vol. 80, p. 084043, Oct 2009.
  * [43] I. Martínez, “The post-Newtonian expansion in the effective theory of compact objects,” 1 2022.
  * [44] S. Foffa and R. Sturani, “Effective field theory calculation of conservative binary dynamics at third post-Newtonian order,” Phys. Rev. D, vol. 84, p. 044031, 2011.
  * [45] S. Foffa and R. Sturani, “Dynamics of the gravitational two-body problem at fourth post-Newtonian order and at quadratic order in the Newton constant,” Phys. Rev. D, vol. 87, no. 6, p. 064011, 2013.
  * [46] R. A. Porto and I. Z. Rothstein, “Next to Leading Order Spin(1)Spin(1) Effects in the Motion of Inspiralling Compact Binaries,” Phys. Rev. D, vol. 78, p. 044013, 2008.  [Erratum: Phys.Rev.D 81, 029905 (2010)].
  * [47] R. A. Porto and I. Z. Rothstein, “Spin(1)Spin(2) Effects in the Motion of Inspiralling Compact Binaries at Third Order in the Post-Newtonian Expansion,” Phys. Rev. D, vol. 78, p. 044012, 2008.  [Erratum: Phys.Rev.D 81, 029904 (2010)].
  * [48] M. Levi, “Next to Leading Order gravitational Spin-Orbit coupling in an Effective Field Theory approach,” Phys. Rev. D, vol. 82, p. 104004, 2010\.
  * [49] R. A. Porto, “Next to leading order spin-orbit effects in the motion of inspiralling compact binaries,” Class. Quant. Grav., vol. 27, p. 205001, 2010.
  * [50] M. Levi and J. Steinhoff, “Spinning gravitating objects in the effective field theory in the post-Newtonian scheme,” JHEP, vol. 09, p. 219, 2015\.
  * [51] R. A. Porto, “The effective field theorist’s approach to gravitational dynamics,” Physics Reports, vol. 633, p. 1–104, May 2016.
  * [52] M. Levi, “Effective Field Theories of Post-Newtonian Gravity: A comprehensive review,” Rept. Prog. Phys., vol. 83, no. 7, p. 075901, 2020.
  * [53] W. D. Goldberger, J. Li, and I. Z. Rothstein, “Non-conservative effects on spinning black holes from world-line effective field theory,” JHEP, vol. 06, p. 053, 2021.
  * [54] P. Hut, J. Makino, and S. McMillan, “Building a better leapfrog,” Astrophys. J. Lett., vol. 443, p. L93, 1995.
  * [55] C. Cutler and E. E. Flanagan, “Gravitational waves from merging compact binaries: How accurately can one extract the binary’s parameters from the inspiral wave form?,” Phys. Rev. D, vol. 49, pp. 2658–2697, 1994.
  * [56] T. Binnington and E. Poisson, “Relativistic theory of tidal Love numbers,” Phys. Rev. D, vol. 80, p. 084018, 2009.
  * [57] R. Jordan, “Effective Field Equations for Expectation Values,” Phys. Rev. D, vol. 33, pp. 444–454, 1986.
  * [58] C. R. Galley, “Classical Mechanics of Nonconservative Systems,” Phys. Rev. Lett., vol. 110, no. 17, p. 174301, 2013.
  * [59] M. C. Miller, F. K. Lamb, and G. B. Cook, “Effects of rapid stellar rotation on equation of state constraints derived from quasi-periodic brightness oscillations,” Astrophys. J., vol. 509, p. 793, 1998.
  * [60] D. G. Ravenhall and C. J. Pethick, “Neutron Star Moments of Inertia,” apj, vol. 424, p. 846, Apr. 1994.
  * [61] P. Bedaque and A. W. Steiner, “Sound velocity bound and neutron stars,” Phys. Rev. Lett., vol. 114, no. 3, p. 031103, 2015.
  * [62] E. Poisson and C. M. Will, Gravity: Newtonian, Post-Newtonian, Relativistic.  Cambridge University Press, 2014.
  * [63] K. Yagi and N. Yunes, “Approximate Universal Relations for Neutron Stars and Quark Stars,” Phys. Rept., vol. 681, pp. 1–72, 2017.
  * [64] E. Poisson and C. Will, Gravity: Newtonian, Post-Newtonian, Relativistic.  Cambridge University Press, 2014.
  * [65] A. A. Starobinskil and S. M. Churilov, “Amplification of electromagnetic and gravitational waves scattered by a rotating ”black hole”,” Sov. Phys. JETP, vol. 65, no. 1, pp. 1–5, 1974.
  * [66] D. N. Page, “Particle Emission Rates from a Black Hole. 2. Massless Particles from a Rotating Hole,” Phys. Rev. D, vol. 14, pp. 3260–3273, 1976.
  * [67] S. Foffa, P. Mastrolia, R. Sturani, C. Sturm, and W. J. Torres Bobadilla, “Static two-body potential at fifth post-Newtonian order,” Phys. Rev. Lett., vol. 122, no. 24, p. 241605, 2019.
  * [68] M. Levi and J. Steinhoff, “Complete conservative dynamics for inspiralling compact binaries with spins at the fourth post-Newtonian order,” JCAP, vol. 09, p. 029, 2021.
  * [69] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments.  Oxford Master Series in Physics, Oxford University Press, 2007.
  * [70] P. C. Peters, “Gravitational Radiation and the Motion of Two Point Masses,” Physical Review, vol. 136, pp. 1224–1232, Nov. 1964.
  * [71] J. Makino, M. Taiji, T. Ebisuzaki, and D. Sugimoto, “GRAPE-4: A Massively Parallel Special-Purpose Computer for Collisional N-Body Simulations,” apj, vol. 480, pp. 432–446, May 1997.
  * [72] B. P. Abbott et al., “Sensitivity of the Advanced LIGO detectors at the beginning of gravitational wave astronomy,” Phys. Rev. D, vol. 93, no. 11, p. 112004, 2016.  [Addendum: Phys.Rev.D 97, 059901 (2018)].