# Observability of dynamical tides in merging eccentric neutron star binaries

János Takátsy<EMAIL_ADDRESS>Institute for Particle and
Nuclear Physics, HUN-REN Wigner Research Centre for Physics, Konkoly–Thege
Miklós út 29-33, 1121 Budapest, Hungary  Bence Kocsis  Rudolf Peierls Centre
for Theoretical Physics, University of Oxford, Clarendon Laboratory, Parks
Road, Oxford, OX1 3PU, UK  St Hugh’s College, University of Oxford, St
Margaret’s Rd, Oxford, OX2 6LE, UK  Péter Kovács  Institute for Particle and
Nuclear Physics, HUN-REN Wigner Research Centre for Physics, Konkoly–Thege
Miklós út 29-33, 1121 Budapest, Hungary

###### Abstract

While dynamical tides only become relevant during the last couple of orbits
for circular inspirals, orbital eccentricity can increase their impact during
earlier phases of the inspiral by exciting tidal oscillations at each close
encounter. We investigate the effect of dynamical tides on the orbital
evolution of eccentric neutron star binaries using post-Newtonian numerical
simulations and constructing an analytic stochastic model. Our study reveals a
strong dependence of dynamical tides on the pericenter distance, with the
energy transferred to dynamical tides over that dissipated in gravitational
waves (GWs) exceeding $\sim 1\%$ at separations $r_{\mathrm{p}}\lesssim 50$ km
for large eccentricities. We demonstrate that the effect of dynamical tides on
orbital evolution can manifest as a phase shift in the GW signal. We show that
the signal-to-noise ratio of the GW phase shift can reach the detectability
threshold of 8 with a single aLIGO detector at design densitivity for
eccentric neutron star binaries at a distance of $40$ Mpc. This requires a
pericenter distance of $r_{\mathrm{p0}}\lesssim 68$ km
($r_{\mathrm{p0}}\lesssim 76$ km) at binary formation with eccentricity close
to 1 for a reasonable tidal deformability and f-mode frequency of 500 and
$1.73$ kHz (700 and $1.61$ kHz), respectively. The observation of the phase
shift will enable measuring the f-mode frequency of neutron stars
independently from their tidal deformability, providing significant insights
into neutron star seismology and the properties of the equation of state. We
also explore the potential of distinguishing between equal-radius and twin-
star binaries, which could provide an opportunity to reveal strong first-order
phase transitions in the nuclear equation of state.

##  I Introduction

Since the first direct observation of a binary black hole (BH) inspiral with
gravitational waves (GWs) in September 2015 Abbott  _et al._ (2016), the
LIGO–Virgo–Kagra (LVK) Collaboration announced the detection of almost 100 GW
signals produced by inspiraling compact object binaries Abbott  _et al._
(2023), and the ongoing observing run already produced more than 117
additional significant candidate GW events111see
https://gracedb.ligo.org/superevents/public/O4/ for the list of candidates.
These events are dominated by binary BH mergers, however, some of them are
attributed to mergers of binary neutron star (NS) and NS–BH systems.
Particularly interesting was the GW event of a binary NS merger, detected in
August 2017, GW170817 Abbott  _et al._ (2017a). Together with the associated
gamma-ray burst and kilonova signals they were the first example of multi-
messenger astronomy with GWs Abbott  _et al._ (2017b, c). This multi-messenger
observation provided a rich laboratory for standard siren measurement of the
cosmic expansion history Abbott  _et al._ (2017d), the investigation of the
equation of state (EoS) of strongly interacting matter through tidal
deformability measurement Abbott  _et al._ (2019), as well as the inference on
the origin of r-process elements Kasen  _et al._ (2017).

In particular, the measurement of the tidal deformability provided strong
constraints on the nuclear EoS, since it excluded the stiffest equations of
state by putting an upper bound on the radius of NSs with a mass of
$1.4~{}\mathrm{M}_{\odot}$ Annala  _et al._ (2018); Most  _et al._ (2018); De
_et al._ (2018); Tews  _et al._ (2018a). NSs, in general, are vital sources of
information about the EoS of strongly interacting matter at low temperature
and intermediate densities. Under such conditions, none of our microscopic
models are predictive. At low densities, below nuclear saturation, the EoS is
well established. Here, two- and three-body nuclear interactions are
determined by experimental data mostly based on nucleon-nucleon scattering and
the properties of light nuclei (e.g. Akmal  _et al._ (1998); Wiringa  _et al._
(1988)). Hadron resonance gas models and Chiral Effective Field Theory provide
different but robust ways to determine the low-density EoS and to estimate
uncertainties Huovinen and Petreczky (2010); Bazavov  _et al._ (2012); Lynn
_et al._ (2016); Tews  _et al._ (2018b). According to state-of-the-art
calculations, the uncertainties of the nuclear EoS above $\sim 1.1\,n_{0}$
become increasingly significant, with $n_{0}$ being the baryon number density
at nuclear saturation Tews  _et al._ (2013, 2019). On the opposite side of the
phase diagram, at very high densities, due to the asymptotic freedom, one can
resort to perturbative QCD methods Kurkela  _et al._ (2010); Mogliacci  _et
al._ (2013). This method gives reliable results at $\mu_{B}\gtrsim 2.5$ GeV,
or equivalently $n_{B}\gtrsim 40\,n_{0}$. In the intermediate density region
the EoS is highly uncertain. While effective models based on the approximate
global chiral symmetry of QCD can be used to estimate the EoS and the position
of phase transitions, these transitions can virtually happen at any density in
this uncertain region. Therefore, observations of NSs and the measurement of
their properties can be important for constraining the EoS in this region.
(e.g. Refs. Annala  _et al._ (2018); Dietrich  _et al._ (2020); Huth  _et al._
(2022); Marczenko  _et al._ (2023); Takatsy  _et al._ (2023)).

Since tidal deformations are sensitive to the mass, radius and internal
structure of NSs, measuring the tidal deformability can be used to constrain
the EoS. During a binary NS’s inspiral the signature of tidal deformation in
the GW signal mainly appears through a change in the phase of the waveform
Flanagan and Hinderer (2008). The dominant effect of the binary evolution is
caused by adiabatic tides, where the distorted NS remains in equilibrium
according to the tidal force exerted by its companion. These adiabatic tides
can be characterized by a single constant for each multipole moment, the tidal
deformability Love (1909); Damour and Nagar (2009); Flanagan and Hinderer
(2008). Adiabatic tides provide a good approximation in the regime where the
eigenfrequencies of the tidal modes are much higher than the frequency of the
acting force. However, finite frequency effects, also known as dynamical
tides, can become relevant when approaching the merger of the binary system
Hinderer  _et al._ (2016); Steinhoff  _et al._ (2016); Schmidt and Hinderer
(2019); Andersson and Pnigouras (2021). The effects of dynamical tides caused
by the most dominant, fundamental-mode (f-mode) tides have been extensively
studied for circular binaries using Newtonian gravity Shibata (1994); Lai
(1994); Kokkotas and Schaefer (1995), and within a time-domain effective-one-
body approach Steinhoff  _et al._ (2016); Schmidt and Hinderer (2019). GW
models that describe the waveform in the frequency domain also provide a more
efficient description from a practical point of view (see e.g. Refs. Ajith
_et al._ (2008); Khan  _et al._ (2016); Field  _et al._ (2014); Bohé  _et al._
(2017)).

While circular compact binary inspirals are expected to be more common,
dynamical tides can be even more important in highly eccentric NS binary
systems due to the broad spectrum of the acting tidal forces. During close
encounters of the NSs, these tidal forces can excite the lowest-frequency
f-mode oscillations of the NSs, which in turn influences the orbital
evolution. Due to the broad spectrum of excitations, unlike in circular
binaries, dynamical tides in eccentric binaries might become relevant well
before the merger of the system, providing an opportunity to be observed Yang
_et al._ (2018); Yang (2019); Vick and Lai (2019).

While field binaries can be born with high eccentricities, GW radiation
efficiently damps eccentricities of binaries inspiraling from an initially
large periapsis distance Peters (1964), leaving nearly circular binaries by
the time they reach the sensitive band of current GW detectors. However,
measurable eccentricities can be produced for compact object binaries in dense
stellar clusters during dynamical interactions, such as single-single
gravitational-wave capture O’Leary  _et al._ (2009); Kocsis  _et al._ (2011);
Gondán  _et al._ (2018a); Takátsy  _et al._ (2019); Rasskazov and Kocsis
(2019); Gondán and Kocsis (2021); Zevin  _et al._ (2021), binary–single
Samsing  _et al._ (2014); Samsing (2018) or binary–binary interactions Zevin
_et al._ (2019); Arca-Sedda  _et al._ (2021), the secular evolution of
hierarchical triple systems Wen (2003); Antonini and Perets (2012); Antognini
_et al._ (2014); Silsbee and Tremaine (2017); Hoang  _et al._ (2018);
Petrovich and Antonini (2017); Rodriguez and Antonini (2018); Fragione and
Kocsis (2019); Fragione  _et al._ (2019), or interactions in gaseous disks
Tagawa  _et al._ (2021); Samsing  _et al._ (2022); Fabj and Samsing (2024); Li
_et al._ (2023); Rowan  _et al._ (2023); Whitehead  _et al._ (2023). While
theoretical predictions to form eccentric compact binaries are quite robust,
their detection and eccentricity measurement are difficult, since template
banks for matched-filtering GW searches typically assume quasi-circular
binaries (e.g. Ref. Messick  _et al._ (2017); Davies  _et al._ (2020); Aubin
_et al._ (2021)). Searches for burst-like GW signals might be effective in
detecting highly eccentric sources Tiwari  _et al._ (2016); Ramos-Buades  _et
al._ (2020); Cornish  _et al._ (2021); Dálya  _et al._ (2021), although their
sensitivity is more difficult to quantify Klimenko  _et al._ (2016). On the
other hand, recent studies have shown that eccentricities could be accurately
measured using matched filtering templates if $e_{\mathrm{10Hz}}\gtrsim 0.1$
Gondán  _et al._ (2018b); Gondán and Kocsis (2019). While the effects of
finite eccentricity and spin precession on the GW signal are difficult to
distinguish Romero-Shaw  _et al._ (2023), recent studies even suggest that
some of the binary BH inspirals observed by the LVK detectors show significant
support for non-negligible eccentricities Gamba  _et al._ (2023); Gayathri
_et al._ (2022); Romero-Shaw  _et al._ (2023, 2022); Gupte  _et al._
(2024).222Note that Ref. Iglesias  _et al._ (2022) found in their analysis
neglecting a possible spin-orbit misalignment that the incorporation of
eccentricity was not preferred compared to the quasi-circular case, however
eccentricty may be preferred for at least one observed system with spin-orbit
misalignment. Once observed, eccentricities can be used as efficient
discriminators between isolated and dynamical binary formation channels (e.g.
Ref. Zevin  _et al._ (2021)). While the merger rate of eccentric binary
systems is highly uncertain, recent studies suggest that it may be significant
Hoang  _et al._ (2018); Rodriguez  _et al._ (2018); Samsing  _et al._ (2018a),
with both numerical and semi-analytical models indicating that they may
account for up to $\sim 10\%$ of the total binary BH merger rate in globular
clusters Samsing (2018); Samsing  _et al._ (2018b); Rodriguez  _et al._
(2018); Zevin  _et al._ (2019); Antonini and Gieles (2020); Kremer  _et al._
(2020); Zevin  _et al._ (2021).

Several previous studies also investigated various aspects of adiabatic and
dynamical tides in binary systems with non-negligible eccentricities and the
corresponding effects on GWs. These include analytic expressions for the
energy and angular momentum transfer, as well as numerical assessments of the
effect of dynamical tides on the orbital evolution Parisi and Sturani (2018);
Yang  _et al._ (2018); Yang (2019); Vick and Lai (2019); Wang and Lai (2020);
Bernaldez and Datta (2023). Resonant excitations of tidal modes in circular
and eccentric binary white dwarfs have also been the subject of multiple
studies Rathore  _et al._ (2005); Fuller and Lai (2011); Burkart  _et al._
(2014). Refs. Yang  _et al._ (2018); Yang (2019) have previously derived the
energy transfer between stellar oscillations and the orbital motion, which we
rederive using a more convenient parameterization of the NS tides. Those
papers have also calculated the GW phase shift caused by dynamical tides in
the small eccentricity limit, finding that it is of the order of
$\mathcal{O}(0.5)$ radians. Vick & Lai Vick and Lai (2019) approached this
problem from a different perspective, solving the equations of orbital motion
numerically for high-eccentricity orbits. Their results suggest these phase
shifts can become an order of magnitude larger for these eccentricities.

In this paper we investigate the effect of dynamical tides on the evolution of
eccentric NS binary systems in more detail using realistic NS parameters and
focus on the detectability of the induced GW phase shift. We run post-
Newtonian numerical simulations and construct an analytic stochastic model to
interpret the results for arbitrary eccentricities. In particular, our
analytical model is based on the Lagrangian of Ref. Flanagan and Hinderer
(2008), which is used to calculate the transmitted energy to tidal oscillation
modes and the corresponding orbital phase shift generated during the dynamical
interactions of binary NSs during each close approach. While resonant
excitations of dynamical tides in NS binaries have been studied before, to our
knowledge this is the first study to investigate the stochastic evolution of
dynamical tides from orbit-to-orbit and the corresponding phase shift. We
confirm that the analytical results and numerical calculations are in good
agreement. Then, we also investigate if the effect of dynamical tides could be
observed with current GW detectors by estimating the signal-to-noise ratio of
the perturbation to the GW. While calculations suggest that quasi-universal
relations hold between the tidal deformability and the f-mode frequency of NSs
Chan  _et al._ (2014); Godzieba  _et al._ (2021); Pradhan  _et al._ (2022,
2023), by measuring the effect of dynamical tides we should be able to measure
these two quantities independently. Therefore, we also assess the effect of
deviating from these relations. Another useful aspect of independent
measurements on tidal deformabilities and f-mode frequencies is the
possibility of observing inspirals of the so-called twin star binaries, i.e.
NS-NS binaries with the same mass but different stellar radii. While adiabatic
tides are known to be only sensitive to the mass-weighted tidal deformability
of the binary in the circular case Flanagan and Hinderer (2008) with finite
eccentricities only negligibly modifying these tides Bernaldez and Datta
(2023), we show that dynamical tides could be used to break the degeneracy and
observationally differentiate between the equal-radius and the twin star
scenarios.

This paper is organised as follows. In Sec. II, we introduce our model for
eccentric binary evolution with tidal interactions, and summarize our
assumptions and the analytical methods. In Sec. III we present our results for
the phase shift caused by dynamical tides, for the comparison with numerical
calculations, as well as for the detectability of the phase shifts. We also
investigate the measurability of strong first-order phase transitions in the
EoS using dynamical tides. Finally, we conclude in Sec. IV.

##  II Tidal model for eccentric binaries

In the simplest approximation, tidal modes can be treated as forced damped
oscillators Lai (1994); Flanagan and Hinderer (2008). The dominant source of
dissipation in the inspiral phase is expected to be gravitational radiation
(see e.g. Refs. Flanagan and Hinderer (2008); Yang  _et al._ (2018)), which we
include in our model. While some of the damping forces caused by bulk or shear
viscosities may be important during the merger and post-merger phases Alford
and Harris (2019); Alford  _et al._ (2020); Most  _et al._ (2022), they can be
neglected during the inspiral phase due to the small viscosities of low-
temperature dense matter Cutler  _et al._ (1990); Kochanek (1992); Bildsten
and Cutler (1992). In the center-of-mass frame dipole moments vanish, while
the coupling of external fields to higher-order multipole moments drops
quickly. Considering the dominant $l=2$ quadrupole oscillations, the driving
force comes from the external quadrupolar tidal field $\mathcal{E}_{ij}$. In
this section we introduce our dynamical model by first examining the problem
with only classical interactions, then complementing it with general
relativistic effects. Throughout the paper we use $G=c=1$.

###  II.1 The Newtonian solution

Let us now consider a binary system with Newtonian interactions, consisting of
a NS and another compact object, which can be either a BH or a NS. The
components have masses $m_{1}$ and $m_{2}$ and Love numbers $\lambda_{1}$ and
$\lambda_{2}$. For simplicity, we only present the equations for
$\lambda_{2}=0$, however, note that our results correspond to binaries, where
both NSs have a finite tidal deformability. The tidal Love number and induced
quadrupole moment of object $1$ can be expressed as a sum of contributions
from $l=2$ modes with various $n\geq 0$ radial nodes,
$\lambda_{1}=\sum_{n}\lambda_{1,n}$ and $Q_{ij}=\sum_{n}Q_{ij}^{n}$, where
$Q_{ij}$ are the components of the trace-free quadrupole moment tensor, which
are related to the density perturbation $Q_{ij}=\int\delta\rho\
(x_{i}x_{j}-\frac{1}{3}\delta_{ij})d^{3}x$ for a Newtonian star. Generally,
however, one can neglect contributions from oscillations with $n>0$ radial
nodes, since they are subdominant ($\lambda_{1,0}\gg\lambda_{1,n}$ and
$\omega_{n}\geq\omega_{0}$, for $n>0$) Flanagan and Hinderer (2008). Hence, we
use $\lambda_{1}\approx\lambda_{1,0}$ and $Q_{ij}\approx Q_{ij}^{0}$.

The Lagrangian of the binary system then reads333see Appendix A for an
explanation of the interaction term:

$\displaystyle L$ $\displaystyle=$ $\displaystyle\left[{1\over
2}\mu{\dot{r}}^{2}+{1\over 2}\mu
r^{2}{\dot{\varphi}}^{2}+{m_{\mathrm{tot}}\mu\over r}\right]-{1\over
2}Q_{ij}{\cal E}_{ij}$ (1) $\displaystyle+{1\over
4\lambda_{1}\omega_{0}^{2}}\left[{\dot{Q}}_{ij}{\dot{Q}}_{ij}-\omega_{0}^{2}Q_{ij}Q_{ij}\right],$

where $m_{\mathrm{tot}}=m_{1}+m_{2}$ and $\mu=m_{1}m_{2}/m_{\mathrm{tot}}$ are
the total and reduced masses, $r$ and $\varphi$ are the relative distance and
true anomaly, respectively, and $\omega_{0}$ is the fundamental f-mode
frequency. Note that the static solution for the tidal deformation reads
$Q_{ij}=-\lambda\mathcal{E}_{ij}$.

The external tidal field can be expressed as

$\displaystyle{\cal E}_{ij}$ $\displaystyle=$ $\displaystyle-
m_{2}\partial_{i}\partial_{j}\left({1\over r}\right)=-{m_{2}\over
r^{3}}(3n^{i}n^{j}-\delta^{ij})=$ (2) $\displaystyle=$
$\displaystyle-{m_{2}\over r^{3}}\begin{pmatrix}{1\over 2}+{3\over 2}\cos
2\varphi&{3\over 2}\sin 2\varphi&0\\\ {3\over 2}\sin 2\varphi&{1\over
2}-{3\over 2}\cos 2\varphi&0\\\ 0&0&-1\end{pmatrix}.$

Although there are nine elements of the above tensor, it can only have five
independent components, due to its symmetrical and traceless nature.
Spherical-harmonic tensors, $\bm{\mathcal{Y}}_{lm}$, give a natural way of
representing such objects, since they form a basis for symmetric traceless
tensors. They are also related to spherical harmonics through the relation

$Y_{lm}(\bm{n})=\sqrt{\tfrac{(2l+1)!!}{4\pi\,l!}}\
\bm{\mathcal{Y}}_{lm}\cdot\bm{n}^{\otimes l}\ ,$ (3)

where $\bm{n}$ is a radial unit vector, $\bm{n}^{\otimes
l}=\overbrace{\bm{n}\otimes\ldots\otimes\bm{n}}^{l}$ is the $l$-fold tensor
product of $\bm{n}$, and the inner product is defined for two arbitrary
tensors $\bm{A}$ and $\bm{B}$ of order $k$ as
$\bm{A}\cdot\bm{B}=A_{a_{1}\ldots a_{k}}B_{a_{1}\ldots a_{k}}$. Another useful
property of spherical-harmonic tensors is that they are orthonormal with the
inner product:

$\bm{\mathcal{Y}}_{lm}\cdot\bm{\mathcal{Y}}_{lm^{\prime}}^{*}=\delta_{mm^{\prime}}\
,$ (4)

while $\bm{\mathcal{Y}}_{l(-m)}=(-1)^{m}\bm{\mathcal{Y}}_{lm}^{*}$. One can
formulate these tensors directly by considering the Cartesian form of the
spherical harmonics, replacing each component $n_{x},n_{y},n_{z}$ of the
radial unit vector $\bm{n}$ with Cartesian unit vectors
$\bm{e}_{x},\bm{e}_{y},\bm{e}_{z}$, the multiplications with tensor products,
and normalizing according to Eq. (4) (which means factoring in the constant in
Eq. (3), see Table 1 for explicit examples).

$lm$ | $Y_{lm}(\bm{n})$ | $\bm{\mathcal{Y}}_{lm}$
---|---|---
$2\,0$ | $\sqrt{\frac{5}{4\pi}}(n_{z}^{2}-n_{\text{\scalebox{0.8}{$\uparrow$}}}n_{\text{\scalebox{0.8}{$\downarrow$}}})$ | $\sqrt{\frac{2}{3}}(\bm{e}_{z}\otimes\bm{e}_{z}-\bm{e}_{\text{\scalebox{0.8}{$\uparrow$}}}\otimes\bm{e}_{\text{\scalebox{0.8}{$\downarrow$}}})$
$2\,1$ | $-\sqrt{\frac{15}{4\pi}}n_{\text{\scalebox{0.8}{$\uparrow$}}}n_{z}$ | $-\sqrt{2}\bm{e}_{\text{\scalebox{0.8}{$\uparrow$}}}\otimes\bm{e}_{z}$
$2\,2$ | $\sqrt{\frac{15}{8\pi}}n_{\text{\scalebox{0.8}{$\uparrow$}}}n_{\text{\scalebox{0.8}{$\uparrow$}}}$ | $\bm{e}_{\text{\scalebox{0.8}{$\uparrow$}}}\otimes\bm{e}_{\text{\scalebox{0.8}{$\uparrow$}}}$
Table 1: Spherical harmonics with Cartesian coordinates and the corresponding
spherical harmonic tensors for $l=2$.
$n_{\text{\scalebox{0.8}{$\uparrow$}}}=\frac{1}{\sqrt{2}}(n_{x}+in_{y})$,
$n_{\text{\scalebox{0.8}{$\downarrow$}}}=\frac{1}{\sqrt{2}}(n_{x}-in_{y})$,
and similarly
$\bm{e}_{\text{\scalebox{0.8}{$\uparrow$}}}=\frac{1}{\sqrt{2}}(\bm{e}_{x}+i\bm{e}_{y})$,
$\bm{e}_{\text{\scalebox{0.8}{$\downarrow$}}}=\frac{1}{\sqrt{2}}(\bm{e}_{x}-i\bm{e}_{y})$.
The normalization factor from Eq. (3) is $\sqrt{8\pi/15}$. Spherical harmonic
tensors for $m=-1,-2$ can be obtained from the identity below Eq. (4).

The tensors $Q_{ij}$ and $\mathcal{E}_{ij}$ are expanded in this basis as

$\bm{Q}=\sum_{m}Q_{m}\bm{\mathcal{Y}}_{2m}\,,\quad\bm{\mathcal{E}}=\sum_{m}\mathcal{E}_{m}\bm{\mathcal{Y}}_{2m}\,,\quad$
(5)

with coefficients

$\mathcal{E}_{0}=-\sqrt{3\over 2}{m_{2}\over r^{3}}\,,\quad\mathcal{E}_{\pm
2}=-{3\over 2}{m_{2}\over r^{3}}e^{\mp i2\varphi}\,,\quad\mathcal{E}_{\pm
1}=0\ ,$ (6)

From the Lagrangian in Eq. (1) the equations of motion for the independent
quadrupole components are

$\displaystyle{\ddot{Q}}_{m}+2\gamma_{0}{\dot{Q}}_{m}+\omega_{0}^{2}Q_{m}$
$\displaystyle=$ $\displaystyle-\lambda_{1}\omega_{0}^{2}\mathcal{E}_{m}\ ,$
(7)

where we have added a finite but small dissipative term with coefficient
$\gamma_{0}$. We identify this term with gravitational radiation caused by
f-mode oscillations. Assuming that $\gamma_{0}$ is small compared to
$\omega_{0}$ one can calculate $\gamma_{0}$ using the radiation power for a
single oscillation mode

$\frac{\mathrm{d}E_{\mathrm{GW,Q}}}{\mathrm{d}t}=-\frac{1}{5}\left\langle{\dddot{Q}_{ij}\dddot{Q}_{ij}}\right\rangle.$
(8)

Then for $\gamma_{0}$ we get

$\gamma_{0}=\frac{1}{5}{\lambda\omega_{0}^{6}}.$ (9)

For realistic NSs, $\gamma_{0}\ll\omega_{0}$ is satisfied.444In our examples
below, we typically assume $\lambda=\Lambda m_{\rm NS}^{5}$ where
$400\leq\Lambda\leq 700$ for which
$\gamma_{0}/\omega_{0}=\frac{1}{5}\Lambda(m_{\rm NS}\omega_{0})^{5}\sim
3\times 10^{-12}$ if $\omega_{0}=2\pi f_{0}$ and $1.61\,\mathrm{kHz}\leq
f_{0}\leq 1.81\,\mathrm{kHz}$, see Sec. II.3.

One should in principle also account for the spins of the NSs, as they can
strengthen dynamical tides Lai (1994); Ho and Lai (1999); Foucart  _et al._
(2019); Ma  _et al._ (2020); Steinhoff  _et al._ (2021). If the NS is rotating
in the retrograde direction with respect to the orbit, the rotation drags the
oscillation mode frequencies to lower values in the inertial frame. This
allows a resonance to occur at a lower orbital frequency. Studies
investigating the effect of spins on dynamical tides have shown that low-
frequency g-mode oscillations, which are weakly coupled to the external tides,
can be enhanced if the NS rotates faster than $100$ Hz, while f-mode resonance
can also significantly increase the phase shift if the spin frequency is
higher than $300-500$ Hz, depending on the EoS Ho and Lai (1999); Lai and Wu
(2006). However, most studies, as well as recent observations suggest that a
high spin rate is unlikely for NS binaries at the time they enter the aLIGO
band, due to spin-down allowed by the long timescale of the binary evolution.
For example, both GW events associated with binary NS inspirals, GW170817 and
GW190425, were found to be consistent with low spin configurations Abbott  _et
al._ (2019, 2020), while the fastest known spinning pulsar in a binary NS
system spins with $44$ Hz Lyne  _et al._ (2004). In less relativistic systems,
such as e.g. white dwarf binaries, external torques acting on dynamical tides
can spin up these objects in a way that they become tidally locked to the
orbit Burkart  _et al._ (2012, 2013, 2014). However, for binary NSs the spin-
up rate is very low compared to the orbital frequency evolution and tidal
locking fails to occur by many orders of magnitude at a wide range of orbital
parameters Burkart  _et al._ (2014); Lai (1994). Additionally, the spins of
isolated field binaries are expected to be aligned with their orbital angular
momentum (e.g. Ref. Farr  _et al._ (2017)). Nevertheless, since high
eccentricities and misaligned spins are both characteristic of dynamical
formation channels, while the inspiral timescale of eccentric binaries are
also much shorter than for circular ones, spin effects might still be
important for these systems. For the time being, we neglect the effect caused
by the rotation of NSs, however, this might be investigated in a follow-up
study.

As a first approximation in describing dynamical tides, we neglect the back-
reaction effect of the tidal deformations on the orbital parameters. Then the
orbital motion will be periodic and hence can be expanded as a Fourier-series
in time, or equivalently, in the mean anomaly, $M$.555To avoid confusion note
that we are using $m_{1,2}$ and $m_{\rm tot}=m_{1}+m_{2}$ to describe mass
parameters; $M$ is used exclusively to denote mean anomaly. For this purpose,
we utilize the so-called Hansen-coefficients, $X_{k}^{q,m}$, which are defined
as:

$\left(\frac{r}{a}\right)^{q}\exp(im\varphi)=\sum\limits_{k=-\infty}^{\infty}X_{k}^{q,m}(e)\
\mathrm{exp}(ikM)\ ,$ (10)

where $q$ is an arbitrary integer and $e$ is the eccentricity of the orbit.
These coefficients are calculated using the following relation (assuming $M=0$
at periapsis):

$X_{k}^{q,m}(e)=\frac{1}{2\pi}\int\limits_{-\pi}^{\pi}\left(\frac{r}{a}\right)^{q}\cos(m\varphi-
kM)\mathrm{d}M\ ,$ (11)

Therefore, the right hand side of Eq. (7) can be expressed as

$-\lambda_{1}\omega_{0}^{2}\mathcal{E}_{m}(t)=C_{m}\lambda_{1}\omega_{0}^{2}\frac{m_{2}}{a^{3}}\sum\limits_{k=-\infty}^{\infty}X_{k}^{-3,-m}e^{ikM}\
,$ (12)

where $C_{m}=3/2$ for $m=\pm 2$ and $C_{m}=\sqrt{3/2}$ for $m=0$. Hence, the
time domain solutions for the tidal components are:

$\displaystyle Q_{m}(t)$
$\displaystyle=C_{m}\lambda_{1}\omega_{0}^{2}\frac{m_{2}}{a^{3}}\times$
$\displaystyle\times\sum\limits_{k=-\infty}^{\infty}\frac{X_{k}^{-3,-m}}{\omega_{0}^{2}-(k\omega_{\mathrm{orb}})^{2}+2i\gamma_{0}k\omega_{\mathrm{orb}}}e^{ikM}\
,$ (13)

where $\omega_{\mathrm{orb}}=\sqrt{m_{\mathrm{tot}}/a^{3}}$ is the orbital
angular frequency. Note, that the implicit time dependency is in $M\equiv
M(t)=\omega_{\mathrm{orb}}t$. Practically, due to the exponential suppression
of the spectrum at large values of $k$ (see Fig. 1), we only need to sum a
finite number of contributions to obtain $Q(t)$ up to a given precision.

Figure 1: Dimensionless Fourier-coefficients of the quadrupolar deformation
$Q_{-2}(t)$ of object $1$, i.e. with respect to the $\bm{\mathcal{Y}}_{lm}$
mode with $l=2$ and $m=-2$, as a function of the angular frequency
$k=\omega/\omega_{\mathrm{orb}}$ (top) and the time dependence of $Q_{-2}(t)$
(bottom) in a binary NS system with $m_{1}=m_{2}=1.4~{}\mathrm{M}_{\odot}$.
The orbital parameters are fixed at $a=600$ km, $e=0.9$, while the tidal
parameters are $\Lambda_{1}=\lambda_{1}/m_{1}^{5}=400$ and $f_{0}=1800$ Hz
(based on Ref. Pradhan  _et al._ (2022)). The orbital angular frequency is
then $\omega_{\mathrm{orb}}\approx 41.5~{}s^{-1}$, which means that the
resonance peak is at $k_{\mathrm{res}}=2\pi f_{0}/\omega_{\mathrm{orb}}\approx
272.7$. The values of the discrete Fourier spectrum are connected to guide the
eye.

Fig. 1 shows an example for the Fourier components of $Q_{-2}$, while the time
dependence of $Q_{-2}$, obtained from Eq. (13) is also shown in the bottom
panel. The Fourier components exhibit a broad plateau due to the high
eccentricity, while there is also a narrow peak corresponding to
$f_{0}=\omega_{0}/2\pi$. We can also see the necessity of an eccentric orbit
for the excitation of dynamical tides. Without an eccentric orbit, the
spectrum of the tidal field would consist of a single peak (at
$-mf_{\text{orb}}$), hence, modes with a frequency of $f_{0}$ would not be
excited unless the system is close to resonance ($f_{\text{orb}}=f_{0}/2$ for
$|m|=2$). However, since $f_{0}\gg f_{\mathrm{orb}}$ during most of the
inspiral of binary NSs, the system would not reach a resonant state before
approaching its merger. From the bottom panel of Fig. 1 one can easily
understand what is happening during eccentric encounters. The compact objects
spend most of their orbit far from each other, hence, tidal forces are weak
for most of the orbit. When the NSs approach pericenter, tidal forces suddenly
increase and then decrease immediately. Therefore, the neutron stars act as
bells which are struck once during each orbit. Tidal modes get a short impulse
during the pericenter passage and then oscillate freely for the rest of the
orbit, while also slowly losing energy due to GW radiation.

###  II.2 General relativistic effects

In addition to nonrelativistic tidal forces, general relativistic effects are
also relevant during this process, such as energy loss due to gravitational
radiation. The orbit of an eccentric binary system with semi-major axis $a$
and eccentricity $e$ and without tidal effects is modified due to
gravitational radiation according to the following equations in the 2.5 Post-
Newtonian (PN) leading order Peters (1964):

$\frac{\mathrm{d}a}{\mathrm{d}t}=-\frac{64}{5}\frac{m_{1}m_{2}m_{\mathrm{tot}}}{a^{3}(1-e^{2})^{7/2}}\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right),$
(14)
$\frac{\mathrm{d}e}{\mathrm{d}t}=-\frac{304}{15}\frac{m_{1}m_{2}m_{\mathrm{tot}}e}{a^{4}(1-e^{2})^{5/2}}\left(1+\frac{121}{304}e^{2}\right).$
(15)

It is worth introducing the dimensionless pericenter distance
$\rho_{\mathrm{p}}=a(1-e)/m_{\mathrm{tot}}$. Eqs. (14) and (15) can be solved
analytically and yield the following relation for $\rho_{\mathrm{p}}(e)$
Peters (1964); O’Leary  _et al._ (2009):

$\rho_{\mathrm{p}}(e)=\frac{c_{0}e^{12/19}}{(1+e)}\left(1+\frac{121}{304}e^{2}\right)^{\frac{870}{2299}},$
(16)

where $c_{0}$ is a dimensionless constant that can be expressed from the
initial parameters, setting $e=e_{0}$ and
$\rho_{\mathrm{p}}=\rho_{\mathrm{p0}}$ in Eq. (16).

The inspiral ends either at physical collision or when reaching the last
stable orbit (LSO) determined byCutler  _et al._ (1994):

$\rho_{\mathrm{p}}(e_{\mathrm{LSO}})=\frac{6+2e_{\mathrm{LSO}}}{1+e_{\mathrm{LSO}}}.$
(17)

The orbit-averaged power radiated through GWs is

$\frac{\mathrm{d}E_{\mathrm{GW}}}{\mathrm{d}t}=-\frac{32}{5}\frac{\mu^{2}m_{\mathrm{tot}}^{3}}{a^{5}(1-e^{2})^{7/2}}\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right),$
(18)

therefore, the energy loss during one orbit is $\Delta
E_{\mathrm{GW}}=\dot{E}_{\mathrm{GW}}\cdot 2\pi/\omega_{\mathrm{orb}}$.

These formulas assume the leading order quadrupole formula for radiation-
reaction where the objects follow Newtonian trajectories with slowly decaying
$a$ and $e$. We examine the impact of other terms in a post-Newtonian (PN)
expansion in powers of $v/c$ which further modify the elliptical orbit and
examine the corresponding GW phase shift in Sec. III.3 (see also Lincoln and
Will (1990); Vick and Lai (2019)). In particular, we account for the
conservative 1PN and 2PN general relativistic corrections to the orbital
motion and the dissipative $2.5$PN GW radiation back-reaction terms, and also
include the backreaction of tidally excited $f$-mode stellar oscillations in
our numerical simulations, given explicitly in Appendix B.

Note that for a proper treatment of the post-Newtonian interactions we would
need to account for additional effects. First of all, in order to correctly
determine the radiation back-reaction terms, we would need to use the
derivatives of the total quadrupole moment instead of calculating the back-
reaction terms separately for the orbit and the tides (see e.g. Refs. Flanagan
and Hinderer (2008); Blanchet (2014)). This would induce cross-terms between
the orbital and tidal evolution. The impact of the post-Newtonian potential on
the tides is another possibly important effect to take into account Vines and
Flanagan (2013); Vines  _et al._ (2011); Henry  _et al._ (2020). In this paper
we neglect these terms and concentrate on the leading-order evolution,
however, these effects could be taken into consideration in a follow-up
analysis.

###  II.3 Modeling NS properties

Both the tidal deformability $\lambda$ and the f-mode frequency $f_{0}$ are
dependent on the EoS describing the matter inside NSs. However, various quasi-
universal relations can be drawn between different macroscopic properties of
NSs Yagi and Yunes (2013). The tidal deformability of $m_{\mathrm{NS}}=1.4$
$\mathrm{M}_{\odot}$ NSs has been recently constrained to
$\Lambda(1.4~{}\mathrm{M}_{\odot})<720$ by the observation of GW170817 Abbott
_et al._ (2019), where $\Lambda=\lambda/m_{\mathrm{NS}}^{5}$. We expect the
$f_{0}$ frequency of $1.4$ $\mathrm{M}_{\odot}$ NSs to fall between 1.5 kHz
and 2.5 kHz Pradhan  _et al._ (2022). In this paper, instead of calculating
tidal deformabilities and f-mode frequencies individually for specific
equations of state, we utilize quasi-universal relations between various
properties of NSs. We also examine possible deviations from these relations
due to first-order phase transitions or quark star solutions.

Most of our calculations are done using binary NSs with two
$1.4~{}\mathrm{M}_{\odot}$ components. In this case we only need to find the
connection between $\Lambda(1.4~{}\mathrm{M}_{\odot})$ and
$f_{0}(1.4~{}\mathrm{M}_{\odot})$. This is done by utilizing the quasi-
universal relation derived in Ref. Pradhan  _et al._ (2022) for NSs containing
nuclear or hyperonic matter. This relation gives us much more accurate results
than using general relations between these and other properties of NSs.

However, in order to be able to calculate f-mode frequencies and tidal
deformabilities for NSs with arbitrary masses, we need to apply a two-step
process. After specifying the mass and the radius of the NS, first we
calculate the tidal deformability from the compactness, $C\equiv M/R$,
inverting the quasi-universal relation presented in Ref. Godzieba  _et al._
(2021). Then, we calculate the f-mode frequency as well, utilizing the
relation derived in Ref. Pradhan  _et al._ (2023) between $C$ and
$M\omega_{0}$. Other parameterizations also exist for these quasi-universal
relations (see e.g. Refs. Pradhan  _et al._ (2022); Yagi and Yunes (2017);
Maselli  _et al._ (2013); Urbanec  _et al._ (2013)), however, we chose these
ones due to their generality in the choice of constructing the EoS, allowing
first-order phase transitions as well. Both of these relations apply within an
error of $\sim 5-10\%$. Crosschecking this two-step process using the quasi-
universal relation between $\Lambda(1.4~{}\mathrm{M}_{\odot})$ and
$f_{0}(1.4~{}\mathrm{M}_{\odot})$, we find that it overestimates the f-mode
frequency by $\sim 5\%$. Since the values of
$\Lambda(1.4~{}\mathrm{M}_{\odot})$ and $f_{0}(1.4~{}\mathrm{M}_{\odot})$ show
a tighter correlation, we assume that $f_{0}(1.4~{}\mathrm{M}_{\odot})$ can be
estimated more accurately from the quasi-universal relation of Ref. Pradhan
_et al._ (2022).

##  III Results

Dynamical tides excited by the eccentric orbital motion may be directly
observed through their GW emission. However, as we will see later, these GW
signals alone are typically below the detection threshold of current GW
detectors. Another way to observe these dynamical tides with GW instruments is
indirectly through their effect on the orbital evolution. While the amplitude
of the GW signals emitted by the binary system is not expected to change much,
the effect of these dynamical tides on the phase shift may be much more
pronounced.

In this section we provide analytical formulas for the energy transferred to
the tidal modes and the phase shift caused by this energy transfer. We also
compare the analytical results with our numerical simulation, calculating the
signal-to-noise ratio (SNR) of the phase shift, and examine the possible
observational implications for the EoS of NSs.

First we construct simple leading-order typically Newtonian analytic models to
interpret the results of the numerical simulations. We investigate post-
Newtonian effects in Sec. III.3.

###  III.1 Energy transmitted to tidal modes

To calculate the energy transferred from the orbital motion to tidal modes,
let us consider the binary system again without GR effects and energy loss due
to GW radiation. Then the change in the energy of a single oscillation mode is

$\displaystyle\Delta E$
$\displaystyle=\Delta\left[\frac{1}{4\lambda_{1}\omega_{0}^{2}}\left(\dot{Q}_{ij}\dot{Q}_{ij}+\omega_{0}^{2}Q_{ij}Q_{ij}\right)\right]=$
$\displaystyle=\frac{1}{2\lambda_{1}\omega_{0}^{2}}\int\limits_{t_{1}}^{t_{2}}\mathrm{d}t\
\dot{Q}_{ij}\left(\ddot{Q}_{ij}+\omega_{0}^{2}Q_{ij}\right)=$
$\displaystyle=-\frac{1}{2}\int\limits_{t_{1}}^{t_{2}}\mathrm{d}t\
\mathcal{E}_{ij}\dot{Q}_{ij}\ .$ (19)

Following Press and Teukolsky (1977), let $t_{1}\rightarrow-\infty$ and
$t_{2}\rightarrow\infty$. Then, using the Fourier-transforms of $\mathcal{E}$
and $Q$ defined as

$f(t)=\int\limits_{-\infty}^{\infty}\mathrm{d}\omega\ e^{-i\omega
t}\tilde{f}(\omega)\ ,$ (20)

we get Press and Teukolsky (1977):

$\displaystyle\Delta
E=\pi\lambda_{1}\omega_{0}^{2}\int\limits_{-\infty}^{\infty}\mathrm{d}\omega\
\frac{-i\omega|\tilde{\mathcal{E}}_{ij}(\omega)|^{2}}{\omega_{0}^{2}-\omega^{2}-i\omega\epsilon}=$
$\displaystyle=\pi^{2}\lambda_{1}\omega_{0}^{2}|\tilde{\mathcal{E}}_{ij}(\omega_{0})|^{2}=\sum_{m}\pi^{2}\lambda_{1}\omega_{0}^{2}|\tilde{\mathcal{E}}_{m}(\omega_{0})|^{2}\
,$ (21)

where in the periodic case from Eq. (12),

$\tilde{\mathcal{E}}_{m}(\omega)=-\sum\limits_{k=-\infty}^{\infty}C_{m}\frac{m_{2}}{a^{3}}X_{-k}^{-3,-m}\delta(\omega-k\omega_{\mathrm{orb}})\
.$ (22)

Hence, in order to get the energy transferred to tidal modes, one needs to
evaluate the Fourier-transform of the tidal field at the resonance frequency.
Note that in case we have a periodic excitation where none of the harmonics
coincide with the resonance frequency, the total energy transmitted is zero.
In contrast, the transmitted energy becomes infinity as these two coincide,
which corresponds to a resonant excitation.

However, we are interested in the energy exchange during a single passage.
Then, the Fourier transform should be limited to only a single orbital period,
while keeping the tidal field zero elsewhere. However, in many cases,
especially close to the merger, where the tidal field at the apocenter is not
negligible, the jump at the boundary results in an oscillatory spectrum, and
therefore the energy exchange will also include this artificial oscillation.
In fact, it is not trivial how to uniquely define the single-passage scenario
without considering the whole evolution of the system. We find that a
convenient way to define the energy transfer is to smoothly interpolate
between the harmonics of the discrete Fourier spectrum. This way, we get a
continuous and smooth spectrum, while this smooth interpolation also creates a
natural smoothing at the edges of the tidal field in the time domain. This way
we get the coefficient in Eq. (21) for all values of $\omega_{0}$ (for which
$\omega_{0}=k\omega_{\rm orb}$ does not necessarily hold) by interpolating the
discrete function:

$\tilde{\mathcal{E}}_{m}(\omega=k\omega_{\mathrm{orb}})=-C_{m}\frac{m_{2}}{\omega_{\mathrm{orb}}a^{3}}X_{-k}^{-3,-m}\>,$
(23)

where the extra factor $1/\omega_{\mathrm{orb}}$ comes from smearing the
$\delta$ functions separated from each other by an $\omega_{\mathrm{orb}}$
wide interval. The change in the orbital energy due to tidal excitations is
then:

$\Delta E_{Q}=-\Delta E\>.$ (24)

Since the tidal oscillations of the two NSs are independent of (not coupled
to) each other, when taking into account the quadrupole deformation of both
objects, the energies extracted from the orbit simply add up.

To understand the amount of energy transmitted to tidal modes during a single
orbit, $\Delta E_{Q}$, it is useful to compare it to the amount of energy
radiated away during a single orbit, $\Delta
E_{\mathrm{GW}}=\dot{E}_{\mathrm{GW}}T_{\rm orb}$ (Eq. 18) where
$T_{\mathrm{orb}}=2\pi/\omega_{\mathrm{orb}}$ is the orbital period. The solid
black lines in Figs. 2 and 3 show the logarithm of $\Delta E_{Q}/\Delta
E_{\mathrm{GW}}$ as a function of $\rho_{\mathrm{p}}$ and $e$ for a binary of
two $1.4~{}\mathrm{M}_{\odot}$ NSs and one with such a NS and a
$10~{}\mathrm{M}_{\odot}$ BH, respectively.

One can observe that there is an exponential-like dependency of $\Delta
E_{Q}/\Delta E_{\mathrm{GW}}$ as a function of $\rho_{\mathrm{p}}$. This can
be attributed to the fact that the Hansen coefficients have an exponential
tail for large values of $k$, and increasing $\rho_{\mathrm{p}}$ means
increasing $\omega_{0}/\omega_{\mathrm{orb}}$ as well. Looking at the
numerical simulation results, the tidal interaction becomes visible above
$\Delta E_{Q}/\Delta E_{\mathrm{GW}}\gtrsim 10^{-3}$. Another observation can
be made by looking at the values of $T_{\mathrm{orb}}/\tau$, where
$\tau=1/\gamma_{0}$ is the damping time. For most of the parameter space
$T_{\mathrm{orb}}/\tau\ll 1$, which means that tidal oscillations do not have
time to attenuate from one pericenter passage to the next. However, the
oscillation may also be cancelled by the subsequent pulse-like excitation if
it comes at the right phase. Due to the continuously changing orbital period,
the phase at which the tidal pulse arrives is mostly random. Thus, tidal
excitations follow a random walk and so will the energy extracted from the
orbit. Similar tidal evolutions were found in Refs. Vick and Lai (2019);
Samsing  _et al._ (2017). The orange line in Fig. 2 shows the evolution when
tidal effects are taken into account which deviate from simulations without
tidal effects for small $\rho_{p}$. Note that this figure corresponds to NSs
with tidal parameters $\Lambda=400$ and $f_{0}=1800$ Hz, and the excitation
energy would be increased if $\Lambda$ was higher or $f_{0}$ lower.

Fig. 3 shows the same for binaries consisting of a $1.4~{}\mathrm{M}_{\odot}$
NS and a $10~{}\mathrm{M}_{\odot}$ BH. The message of this figure is clear.
Dynamical tides play a much less prominent role for NS-BH binaries. The reason
is twofold. There is only one NS, which can go through substantial tidal
oscillations, and the orbital frequency of the system for given orbital
parameters become lower, which, due to the exponential dependence of the
Hansen coefficients on $k$ results in much lower excitation energies. Fig. 3
also shows that NSs are not expected to be tidally disrupted in a NS–BH binary
unless the BH is highly spinning and the EoS of matter inside the NS is stiff
at low densities – which would result in a larger radius and $\Lambda$, and a
lower f-mode frequency.

Figure 2: The logarithm of $\Delta E_{Q}/\Delta E_{\mathrm{GW}}$ as a function
of orbital parameters $\rho_{\mathrm{p}}$ and $e$ for two $1.4$
$\mathrm{M}_{\odot}$ NSs (solid black lines). The dashed red lines show the
orbit-averaged evolution of binaries with selected initial conditions
$e_{0}\approx 1$ and $\rho_{\mathrm{p0}}=(10,15,20)$, driven by Eqs. 14-15.
Black circles show the evolution of compact objects along the same line of
evolution without any tidal interactions. These come from our numerical
simulation, and each circle corresponds to an apocenter passage. The orange,
light green and dark green dashed lines with circles correspond to compact
binaries with the same initial parameters but with tidal interactions also
taken into account. The tidal deformability and the f-mode frequency was
chosen as $\Lambda=400$ and $f_{0}=\omega_{0}/2\pi=1800$ Hz. The blue line
represents parameters at the last stable orbit, while the vertical orange line
corresponds to direct collision of the two NSs, where we used the very
conservative estimate $R_{\mathrm{NS}}(1.4~{}\mathrm{M}_{\odot})=14$ km.
Dashed black lines show different values of $T_{\mathrm{orb}}/\tau$. For each
selected orbital evolution we display the sky position and orientation
averaged detection range with a single aLIGO detector, $d_{\rm L,det}$,
calculated from the orbit-averaged solution O’Leary  _et al._ (2009). Figure
3: Same as Fig. 2 but for a binary consisting of a NS with
$m_{1}=1.4~{}\mathrm{M}_{\odot}$ and a BH with
$m_{2}=10~{}\mathrm{M}_{\odot}$. The three selected orbital evolutions here
have $e_{0}\approx 1$ and $\rho_{\mathrm{p0}}=(5,6,8)$. Black circles
corresponding to an evolution without any tidal interaction were omitted,
since they would overlap with the colored circles. The tidal parameters of the
NS are $\Lambda=400$ and $f_{0}=1800$ Hz.

###  III.2 Phase shift caused by tidal interactions

Since the energy transmitted to tidal modes remain low during most of the
orbital evolution and can only exceed $\Delta E_{\mathrm{GW}}/10$ during the
last couple of orbits, and since the frequency of the tidal modes falls in the
kHz band where aLIGO is less sensitive, GWs produced by these oscillations
have a low SNR. However, due to back-reaction, the orbital evolution can be
modified substantially. GW detectors are particularly sensitive to changes in
the GW phase. We derive analytical formulas for the phase shift caused by the
back-reaction from tidal excitations and examine their detectability by GW
detectors. A common rule of thumb is to check if the phase shift exceeds 1
radian at any time during the detectable part of the evolution. If it does,
and the binary is detectable, then the phenomenon causing the phase shift is
expected to be detectable.

A common way to quantify the phase shift is to determine its dependence on the
orbital frequency Flanagan and Hinderer (2008); Kocsis  _et al._ (2011);
Chirenti  _et al._ (2017); Yang  _et al._ (2018); Yang (2019). Some other
studies investigated the effect of g-mode excitations in eccentric binaries or
resonant g-mode excitations in circular ones Lai (1994); Xu and Lai (2017);
Andersson and Ho (2018); Kuan  _et al._ (2021). Here we follow a similar
method. Let us consider a binary system whose orbital energy is decreasing due
to GW emission and some energy is transferred to tidal modes. Then, the change
in the mean anomaly during an infinitesimally short period of time is

$\displaystyle\mathrm{d}M$
$\displaystyle=\omega_{\mathrm{orb}}\>\mathrm{d}t=\omega_{\mathrm{orb}}\frac{\mathrm{d}E_{\mathrm{orb}}}{\dot{E}_{\mathrm{orb}}}=\omega_{\mathrm{orb}}\frac{\mathrm{d}E_{\mathrm{orb}}}{\dot{E}_{\mathrm{GW}}+\dot{E}_{Q}}\approx$
$\displaystyle\approx\omega_{\mathrm{orb}}\frac{\mathrm{d}E_{\mathrm{orb}}}{\dot{E}_{\mathrm{GW}}}\left(1-\frac{\dot{E}_{Q}}{\dot{E}_{\mathrm{GW}}}\right)\>.$
(25)

The phase difference between systems with and without tidal excitations is:

$\mathrm{d}(\delta
M)\approx-\omega_{\mathrm{orb}}\frac{\mathrm{d}E_{\mathrm{orb}}}{\dot{E}_{\mathrm{GW}}}\frac{\dot{E}_{Q}}{\dot{E}_{\mathrm{GW}}}\approx-\omega_{\mathrm{orb}}\frac{\mathrm{d}E_{Q}}{\dot{E}_{\mathrm{GW}}}\>,$
(26)

where $\mathrm{d}E_{Q}$ is the change in the energy due to tidal oscillations
while the orbital energy changes by $\mathrm{d}E_{\mathrm{orb}}$. Note that
the integration variable is $E_{\mathrm{orb}}$ in both cases, hence, this
formula applies for orbital phases compared at the same orbital energies, or
equivalently, the same orbital frequency. This will become important later,
when comparing this formula to numerical results. Also note that since both
$\mathrm{d}E_{Q}$ and $\dot{E}_{\mathrm{GW}}$ are negative, the phase
difference will be negative as well, i.e. with dynamical tides it takes fewer
orbits to get to a certain orbital frequency.

During a single orbit, the phase shift then becomes

$\delta M_{1}\approx-\omega_{\mathrm{orb}}\frac{\Delta
E_{Q}}{\dot{E}_{\mathrm{GW}}}=-2\pi\frac{\Delta E_{Q}}{\Delta
E_{\mathrm{GW}}}\>,$ (27)

where $\Delta E_{Q}/\Delta E_{\mathrm{GW}}$ is the quantity derived in the
previous section.

Now we also need to take into account, that the amplitude of tidal excitations
follows a random walk due to the random phase at which the short tidal pulses
arrive. The $j^{\mathrm{th}}$ pulse can be characterized by the complex number
$\Delta q_{j}=A_{j}\mathrm{exp}(i\varphi_{j})$, where $\varphi_{i}$ is a
random phase in the $[0,2\pi]$ interval. Thus, the expectation value of the
complex amplitude after the $j^{\mathrm{th}}$ pulse is

$\langle q_{j}\rangle=\sum\limits_{k=1}^{j}\langle\Delta q_{k}\rangle=0.$ (28)

Since the energy of tidal modes is quadratic in the amplitude, its expectation
value over the incoherent random realizations of the tidal pulses is

$\displaystyle\langle E_{Q,j}\rangle$
$\displaystyle=C\left\langle\left|\sum\limits_{k=1}^{j}\Delta
q_{k}\right|^{2}\right\rangle=C\sum\limits_{k=1}^{j}\left\langle\left|\Delta
q_{k}\right|^{2}\right\rangle=$
$\displaystyle=C\sum\limits_{k=1}^{j}A_{k}^{2}=\sum\limits_{k=1}^{j}\Delta
E_{Q,k}\>,$ (29)

where we have used the fact that the expectation value of cross terms gives
zero. $\Delta E_{Q,k}$ is the energy transmitted to tidal modes with zero
initial amplitude. Utilizing Eq. (27),

$\langle\delta M\rangle=\sum\limits_{k=1}^{j}\delta
M_{1,k}\approx\sum\limits_{k=1}^{j}\left(-2\pi\frac{\Delta E_{Q}}{\Delta
E_{\mathrm{GW}}}\right)_{k}\>.$ (30)

Hence, to get the expectation value of the phase shift, we must add up the
phase shift contributions of each orbit. Note that this result would be the
same if the tidal modes were strongly attenuated.

Let us also utilize the orbit-averaged evolution of the binary without tidal
interactions in Eq. (25) to get an approximate integral formula as

$\displaystyle\mathrm{d}(\delta M)$
$\displaystyle=-\omega_{\mathrm{orb}}\frac{dE_{\rm orb}}{df_{\rm orb}}df_{\rm
orb}\frac{\dot{E}_{Q}}{\dot{E}_{\mathrm{GW}}^{2}}+\mathcal{O}\left(\frac{\dot{E}_{Q}^{2}}{\dot{E}_{\mathrm{GW}}^{3}}\right)$
$\displaystyle=\frac{2\pi}{3}\mathcal{M}^{5/3}\omega_{\rm orb}^{2/3}df_{\rm
orb}\frac{\dot{E}_{Q}}{\dot{E}_{\mathrm{GW}}^{2}}+\mathcal{O}\left(\frac{\dot{E}_{Q}^{2}}{\dot{E}_{\mathrm{GW}}^{3}}\right)\>,$
(31)

where we used the leading order (Newtonian) orbital energy $E_{\rm
orb}=-\frac{1}{2}\eta m_{\rm tot}^{2}/a=-\frac{1}{2}\eta m_{\rm
tot}^{5/3}\omega_{\rm orb}^{2/3}=-\frac{1}{2}\mathcal{M}^{5/3}\omega_{\rm
orb}^{2/3}$ and Kepler’s law $a=m_{\rm tot}^{1/3}\omega^{-2/3}$, hence
$\omega_{\rm orb}dE_{\rm orb}/df_{\rm orb}=2\pi\omega_{\rm orb}dE_{\rm
orb}/d\omega_{\rm orb}=2\pi\frac{2}{3}E_{\rm
orb}=-\frac{2}{3}\pi\mathcal{M}^{5/3}\omega_{\rm orb}^{2/3}$. Thus:

$\langle\delta
M\rangle(f)\approx\frac{(2\pi)^{5/3}}{3}\int\limits_{f_{\text{in}}}^{f}\mathrm{d}f_{\mathrm{orb}}\>\mathcal{M}^{5/3}f_{\mathrm{orb}}^{2/3}\frac{\dot{E}_{Q}}{\dot{E}_{\mathrm{GW}}^{2}}\>,$
(32)

where $\mathcal{M}=m_{\rm tot}^{2/5}\mu^{3/5}$ is the chirp mass,
$\dot{E}_{\mathrm{GW}}$ is given by Eq. (18) and $\dot{E}_{Q}=\Delta
E_{Q}/T_{\mathrm{orb}}$, where $\Delta E_{Q}$ is given by Eqs. (21), (23), and
(24). Note that these expressions only account for the phase shift caused by
dynamical tides and do not include adiabatic tides. The effect of adiabatic
tide can be calculated separately, and has been done by Ref. Flanagan and
Hinderer (2008) for circular and by Ref. Bernaldez and Datta (2023) for
eccentric binaries.

Fig. 4 shows the orbital phase shift as a function of orbital frequency for
binary NSs with two $1.4~{}\mathrm{M}_{\odot}$ components. The figure shows
the results for several different $\Lambda-f_{0}$ pairs and two different
orbital configurations. The solid lines denote results obtained from Eq. (30),
where the input orbital parameters (i.e. $a$ and $e$ to be used for
calculating $\dot{E}_{Q}$ and $\dot{E}_{\mathrm{GW}}$) were extracted from our
numerical simulations in the adiabatic tide limit, i.e. with
$f_{0}=f_{\mathrm{ad}}\gg f_{\mathrm{orb}}$666Since
$f_{\mathrm{orb}}^{\mathrm{max}}\leq 800$ Hz in our configurations, we find
that the adiabatic limit is already approached for $f_{\mathrm{ad}}\geq 2.5$
kHz, while keeping the integration time reasonably low., while the dashed
lines denote results obtained with Eq. (32), for which the orbit-averaged
evolution was used without taking into account the effect of adiabatic tides.
There is a good agreement between the two formulae at moderately high orbital
frequencies, however, at the highest frequencies the orbit-averaged formula
overestimates the results calculated from our numerical simulations. This is
due to the effect of adiabatic tides, which influence the orbital evolution
even in the $f_{0}\gg f_{\mathrm{orb}}$ limit. It is clearly visible that
increasing the tidal deformability has a significant effect on the magnitude
of the phase shift. Changing the tidal deformability from $\Lambda=400$ to
$700$ increases the phase shift by a factor of $\sim 5$ near the merger of the
binary. This is due to the decrease of the f-mode frequency with increasing
$\Lambda$, which has a significant effect on the amount of energy transmitted
to tidal oscillation modes.

Figure 4: The mean phase shift $\langle\delta M\rangle$ caused by tidal
oscillations as a function of the orbital frequency for binary NSs with two
$1.4~{}\mathrm{M}_{\odot}$ components, calculated using Eq. (30) with the
input orbital parameters ($a$ and $e$, see Appendix B) extracted from our
numerical simulations run in the adiabatic tide limit, i.e. with
$f_{0}=f_{\mathrm{ad}}\gg f_{\mathrm{orb}}$ (solid lines) and the orbit-
averaged analytical calculation Eq. (32) (dashed lines). Different colors
correspond to different equations of state, and hence different values for
$\Lambda$ and $f_{0}$. For $\Lambda=(400,500,600,700)$ the corresponding
values for $f_{0}$ are $1.81$ kHz, $1.73$ kHz, $1.66$ kHz and $1.61$ kHz,
respectively. The filled circles denote the last apocenter passage before the
two NSs collide directly, which we defined as their distance dropping below
$28$ km. The black vertical lines denote the initial orbital frequencies for
the two different orbital evolutions.

An eccentric binary emits GWs in a wide range of frequencies. The GW spectrum
of such a binary exhibits a broad peak centered at approximately Wen (2003):

$f_{\mathrm{GW}}=n_{\mathrm{peak}}f_{\mathrm{orb}}\approx\frac{2(1+e)^{1.1954}}{(1-e^{2})^{3/2}}f_{\mathrm{orb}}\>.$
(33)

Note that this formula gives $2f_{\mathrm{orb}}$ for $e=0$, as is expected for
the circular case, while it tends to infinity as $e\to 1$. Fig. 5 shows the
dependence of the phase shift on the peak GW frequency $f_{\mathrm{GW}}$.
Although the initial orbital frequencies of the binary systems shown are below
$10$ Hz, the peak GW frequency is above $200$ Hz for both cases due to the
high initial eccentricity $e_{0}=0.9$. This initial peak GW frequency becomes
lower as we increase the initial semimajor axis, and so does the magnitude of
the phase shift.

Figure 5: Same as Fig. 4 but with the phase shift $\delta M$ shown as a
function of the peak GW frequency $f_{\mathrm{GW}}$ (Eq. 33). The inset shows
the high-frequency part of the evolution on a linear scale.

###  III.3 Comparison to numerical simulations

As already described in Sec. II, the Newtonian equations of motion were
derived from the Lagrangian in Eq. (1). We complemented these equations by
adding the dissipative term proportional to $\dot{\bm{Q}}$, which corresponds
to gravitational radiation generated by the f-mode oscillations. Additionally,
we have included the 2.5PN terms accounting for the orbital GW radiation back-
reaction. In the previous discussions we neglected the conservative 1PN and
2PN corrections to the orbital motion. We examine the effects of these post-
Newtonian terms at the end of this section in both the numerical simulations
and the analytical expressions derived below. The equations of motion are
given in Appendix B. The dynamical variables were transformed into a basis
where all the equations of motion are real. We solve these equations using a
fourth-order Runge-Kutta differential equation solver with adaptive timesteps.
The simulations are initialized at the apocenter with initial conditions
calculated from the Keplerian parameters, while the initial tidal deformations
are set to zero. The simulations are terminated if the separation becomes
smaller than $28$ km, which implies direct collision with a conservative NS
radius of $14$ km.

To compare the numerical simulation results with the analytical calculations a
few more additional steps are needed to extract the shift in the mean anomaly
for a given orbital frequency. Both $M$ and $f_{\mathrm{orb}}$ only exist in
the Keplerian picture and require special treatment to determine from the
numerical simulations, which also include post-Newtonian and tidal interaction
terms. We utilize the fact that the orbit is Keplerian to a good approximation
at the apocenter. Therefore, we calculate all the Keplerian parameters at each
apocenter passage using the energy and angular momentum of the binary system
(see Appendix B), from which the orbital frequency can be determined as well.
The mean anomaly increases by $2\pi$ at each apocenter passage777We find that
due to numerical inaccuracies, more accurate results can be obtained by
directly calculating the mean anomaly from the orbital parameters (see
Appendix B).. It was stressed below Eq. (26) that the analytical formula for
the phase shift is meaningful when the orbital evolutionary tracks are
compared at the same orbital frequency. However, the orbital frequency for the
evolution with and without dynamical tides slightly differ for a given
apocenter passage. Hence, in order to obtain the phase shift at a given
frequency, we interpolate the mean anomaly between adjacent apocenter
passages.

In addition to dynamical tides, the phase of the binary orbit is also
influenced by adiabatic tides, which can be obtained by taking the adiabatic
limit, $f_{0}\gg f_{\mathrm{orb}}$. For circular binaries the orbital phase
shift introduced by adiabatic tides was calculated by Flanagan & Hinderer to
be $\mathcal{O}(1)$ at $f_{\mathrm{orb}}\gtrsim 100$ Hz Flanagan and Hinderer
(2008), while it has been recently shown that the additional phase shift is
$\mathcal{O}(10^{-4}-10^{-2})$ in the adiabatic limit for moderate
eccentricities Bernaldez and Datta (2023). Since the analytic formulas, Eqs.
(30) and (32) only account for the extra phase shift caused by dynamical
tides, we extract this additional phase shift from our numerical simulations
as well. We calculate the excess phase shift caused by dynamical tides by
subtracting the phase shift from adiabatic tides

$\delta
M(f_{\mathrm{orb}})=M_{f_{0}}(f_{\mathrm{orb}})-M_{f_{\mathrm{ad}}}(f_{\mathrm{orb}})\>,$
(34)

where $M_{f_{0}}$ and $M_{f_{\mathrm{ad}}}$ denote the mean anomaly for f-mode
frequency $f_{0}$ and $f_{\mathrm{ad}}\gg f_{\mathrm{orb}}$, respectively. For
practical purposes it is sufficient to set $f_{\mathrm{ad}}$ to be a few times
$f_{0}$, due to the exponential suppression of tidal excitations with
increasing $f_{0}$ (see Fig. 1 and the discussion therein). We obtain $\delta
M$ as follows. For each set of initial parameters we run two separate
simulations, one with $f_{0}$ being appropriate for dynamical tides to be
relevant, and another initialised by setting $f_{0}=f_{\mathrm{ad}}\gg
f_{\mathrm{orb}}$. Then we subtract the mean anomalies evaluated at the same
frequencies. The expectation value of the phase shift evaluated in this way
should essentially be the same as the phase shift shown by the solid line in
Fig. 4, since there the orbital parameters $a$ and $e$ were already calculated
from the numerical simulation taking into account adiabatic tides as well. The
only region where we expect the numerical results to deviate from the analytic
calculations is where the energy transmitted to dynamical tides becomes
comparable to the energy radiated through GWs, and thus where our perturbative
assumption breaks down.

Fig. 6 shows the results obtained from this procedure for binary NSs with two
$1.4~{}\mathrm{M}_{\odot}$ components and three different configurations of
initial semi-major axis and tidal deformability. The initial eccentricity is
$e_{0}=0.9$ for all cases. The solid lines represent the semi-analytical
results obtained from Eq. (30), utilizing the numerically calculated orbital
evolution without dynamical tides, while the dashed lines represent the
numerical results obtained from the procedure described in the previous
paragraphs. For each combination of binary parameters, two results from
numerical simulations with slightly different parameters are shown. Due to the
random walk of phase shifts between consecutive orbits, the difference between
the two evolutions can become considerable, however, they both follow the
trend given by the semi-analytic formula. Even though the magnitude of the
phase shifts for the three different configurations differs by $2-3$ orders of
magnitude at low orbital frequencies, their difference shrinks to a single
order of magnitude by the merger. The phase shifts approach or even exceed 1
radian by the time the binaries merge for all three configurations.

Figure 6: The excess phase shift caused by dynamical tides $\delta M$ as a
function of the orbital frequency for binary NSs with two
$1.4~{}\mathrm{M}_{\odot}$ components obtained from the numerical simulation
(dashed lines) and calculated using the Eq. (30) (solid lines). For a given
color the two numerical results with different tones correspond to evolutions
with slightly different parameter configurations. $e_{0}=0.9$ for each case.

Note that a related measurable quantity is the phase shift measured at a
specific time, e.g. at each apocenter passage of the orbital evolution without
dynamical tides.888In this case the phase shift would be the difference
between the phases of the orbital evolutions with and without dynamical tides,
both measured at the same time defined by the apocenter passage in case of no
dynamical tides. Note that this quantity differs from the one shown in Figs.
4-7, since orbits with and without dynamical tides reach a given frequency at
different times.

Figure 7: The phase shift $\delta M$ as a function of the orbital frequency
for binary NSs with two $1.4~{}\mathrm{M}_{\odot}$ components, averaged for
$\sim 50$ orbital evolutions with slightly different parameters (error bars).
The colored regions show the areas covered by the results calculated from Eq.
(30), with the dashed lines representing their average. For the results with
$a_{0}=800$ km and $\Lambda=600$, the f-mode frequency was taken to be off the
quasi-universal relation, with $f_{0}=1.5$ kHz. $e_{0}=0.9$ for each case.

Given the stochastic nature of the evolution, particularly the sensitivity of
the impulsive tidal kicks at close approaches to the actual phase of the
stellar oscillations, we examine the expectation value and scatter of the
phase shift for an ensemble of simulations with slightly different initial
conditions. This is performed with care given the subtle dependence on the
variations of the initial conditions. Interestingly, for non-zero initial
deformations, varying the initial phase of the quadrupolar oscillation modes
does not change the time intervals between pericenter passages, which
ultimately determine the stochastic evolution of the tidal amplitudes for a
given $f_{0}$. The evolution of the mode amplitudes may change by modifying
either the f-mode frequency or the orbital evolution through $a_{0}$ or
$e_{0}$. Also note that variations in $e_{0}$ and $a_{0}$ are not independent,
changing $e_{0}$ to $e_{0}^{\prime}$ is equivalent to (i.e. results in the
same evolutionary path) modifying $a_{0}$ appropriately to some
$a_{0}^{\prime}$ while keeping $e_{0}$ fixed, i.e. along the evolutionary
paths shown in Figure 2. Hence, we are left with two independent parameters,
e.g. $f_{0}$ and $a_{0}$, that can be independently changed in order to get a
different evolution for the tidal amplitudes. The variation in these
parameters must also be kept sufficiently low in order to simulate
approximately the same orbital evolution, but variations chosen to be too
small may lead to negligible deviations in the evolution. These constraints
limit the range of orbital configurations which can be averaged over.

Fig. 7 shows the averaged values of phase shifts with error bars for different
parameter configurations, obtained by averaging the phase shifts from $\sim
50$ numerical simulations with slightly different parameter configurations.
The regions covered by the phase shifts calculated from the semi-analytical
formula Eq. (30) are indicated by the colored regions, with their average also
indicated by the dashed lines. For most parts of the results there is a good
agreement between the numerical averages and the analytical results, however,
there are a few exceptions. The semi-analytical results overestimate the
numerical average close to the merger for the simulations with $a_{0}=600$ km.
This can be attributed to the fact that Eq. (30) was obtained from a
perturbative approach, while the phase shifts at these stages can no longer be
considered small perturbations. Another discrepancy can be noticed at the low
frequency end of the parameter configuration with $a_{0}=800$ km and
$\Lambda=700$, where the numerical phase shifts are larger than the analytical
ones. This is due to numerical error which causes the phase difference to be
higher than $10^{-5}-10^{-3}$ depending on the orbital frequency even when
they should be lower according to the analytic formulas. To assess the extent
of these limitations we calculated the phase difference for a single
configuration with the same f-mode frequency, but with two different choices
for numerical precision. Since we compare the same evolution to itself, the
phase difference should be zero, however the numerical results differ. The
results from two different such simulations are shown by the dotted black
lines in Fig. 7. Even though this numerical artefact is well below the
observable region and while it is slowly suppressed by increasing the
numerical precision, it is not fully understood why it emerges. They might be
linked to the frequency interpolation of phases.

The semi-analytical phase shift estimates provide an accurate interpretation
of the numerical simulations also for tidal parameters off the quasi-universal
relation.

As we mentioned in the beginning of this subsection, a correct treatment of
the post-Newtonian interactions would require taking into account additional
corrections to equations of motion used in the numerical simulations and the
analytical formulas. For a simple test, we estimate the change in the phase
shift evolution caused solely by the 1PN and 2PN terms added to the orbital
evolution equations (see Eqs. 53 and 54. Fig. 8 shows the results for the
phase shift as a function of the orbital frequency when the 1PN and 2PN terms
are included in the orbital equations of motion. The numerical evaluation of
the phase shift is described in Appendix B. $f_{\mathrm{orb}}$ was calculated
using the PN corrections described in Ref. Blanchet (2014). There is a good
overall agreement between the analytic and numerical results, one major
difference being that the phase shifts for $a_{0}=800\,$km and $\Lambda=700$
are systematically increased compared to the analytic results. Since the
relative scatter of phase shifts is also reduced in this case, we think that
it might be caused by an additional, non-stochastic effect, although the exact
reason remains unclear. It should also be noted that with the inclusion of the
1PN and 2PN corrections, the orbital frequency at merger is reduced by about a
factor of two, resulting in a lower phase shift at merger, while the total
number of orbits is also reduced.

Figure 8: The phase shift $\delta M$ as a function of the orbital frequency
for binary NSs with two $1.4~{}\mathrm{M}_{\odot}$ components, averaged for
$\sim 50$ orbital evolutions with slightly different parameters (error bars).
The orbital evolutions were calculated using the 1PN and 2PN terms as well.
The colored regions show the areas covered by the results analytically
calculated from Eq. (30) without the 1PN and 2PN terms, with the dashed lines
representing their average. For the results with $a_{0}=800$ km and
$\Lambda=600$, the f-mode frequency was taken to be off the quasi-universal
relation, with $f_{0}=1.5$ kHz. $e_{0}=0.9$ for each case.

###  III.4 Comparison to the literature

Our results differ from the results of Vick and Lai Vick and Lai (2019) who
studied the GW phase shift in the time domain for binary NS systems. Their
Fig. 5 showed spikes in the GW phase shift around each pericenter passage,
which exceeded $\sim 10$ radians close to the merger. They also investigated
the differences in GW phases measured at the merger of the binary as a
function of the initial eccentricity in their Fig. 7. The latter figure showed
that the GW phase shift at merger has an oscillatory behaviour as a function
of $e_{0}$, with this phase shift also changing sign at high eccentricities.
We believe that these phase shift spikes and the oscillatory $e_{0}$
dependence are not directly measurable. Fig. 5 in Ref. Vick and Lai (2019) was
obtained by matching the GW phases at the start of their simulations, but in a
realistic scenario the GW signal of the template is matched to minimize the
SNR of the difference relative to the measured data. The corresponding
relevant quantity to predict the expected significance of the tidal
perturbation in a measurement is the difference SNR of the waveforms with and
without tidal effects. For practical purposes this approximately corresponds
to setting the time-shift to match the waveforms at the maximum strengths
(i.e. at pericenter passages close to the merger). This is expected to
diminish the spikes seen in Fig. 5 of that paper and the overall phase shift
is expected to be reduced. Larger phase differences then appear at lower
frequencies, where the GW signal is weaker and the detectors might be much
less sensitive. As for Fig. 7 in Ref. Vick and Lai (2019), since the GW phase
at merger is highly dependent on the assumptions of the merger, the
oscillatory behaviour shown in that figure is very sensitive to the uncertain
assumptions of the simulation and may not be a robust indicator of the
properties of the real physical system. In fact, the condition $r=2.5R$ used
in the simulations as the merger condition in that paper, with $R$ being the
NS radius, is only a rough estimate, similar to the one we used in our
simulations. Indeed, since dynamical tides remove energy from the system, the
frequency of the orbital evolution with dynamical tides at a given time can
only be larger than the frequency of the evolution without dynamical tides.
This means that the phase difference will always be positive999The definition
of the phase shift in that paper contains a minus sign, unlike in our
definition., as is the case in Fig. 5 of Ref. Vick and Lai (2019). Thus, the
negative values shown in Fig. 7 of Ref. Vick and Lai (2019) may be a
consequence of comparing the phases of the GW signals at different times based
on the arbitrary choice of the merger condition. A more adequate assessment of
the phase shift is shown in Fig. 6 of Ref. Vick and Lai (2019), which shows
the cumulative time advance of the orbit with dynamical tides at each
pericenter passage relative to the time of pericenter passages without such
tides, as a function of the number of consecutive pericenter passages. Here
the time advances monotonically increase with each subsequent orbit and the
corresponding phase advances of the mean anomaly at a given time can be
estimated to be $\mathcal{O}(1)$ near the merger, which is consistent with our
results.

It is also worth noting that several other studies consider a different
parameterization of the Lagrangian, utilizing the tidal overlap integral
$Q_{\xi}$ (see e.g. Refs. Lai (1994); Ho and Lai (1999); Vick and Lai (2019);
Yang (2019)). The connection between $\lambda$ and $Q_{\xi}$ is given by
$\lambda\approx 4\pi Q_{\xi}^{2}/(15\omega_{0}^{2})$ using the definition of
$Q_{\xi}$ as given in Ref. Yang (2019).

###  III.5 Detectability of the tidal oscillations

While the orbital phase shift can provide a general guide about the
observability of dynamical tides, a more rigorous test of detectability may be
obtained by calculating the SNR betwen the difference of the evolutionary
curves with and without the dynamical tides. We calculate the $h_{+}$ and
$h_{\times}$ components of the GW signal generated by the orbital motion using
the quadrupole formula, $h_{ij}=2\ddot{I}_{ij}/d$, where $d$ is the distance
of the binary and $I_{ij}$ is the trace-free quadrupole moment. We also
calculate the GW signal emitted by the stellar oscillations, $h^{Q}_{+}$ and
$h^{Q}_{\times}$. For GWs propagating perpendicularly to the orbital plane,
along the $z$ axis, the two polarizations in the transverse–traceless gauge
become:

$h_{+}=\frac{1}{d}\left(\ddot{I}_{xx}-\ddot{I}_{yy}\right)\>,\quad
h_{\times}=\frac{2}{d}\ddot{I}_{xy}\>,$ (35)

which gives:

$\displaystyle h_{+}=\frac{2\mu}{d}\left(\dot{r}^{2}\cos
2\varphi+r\ddot{r}\cos 2\varphi-4r\dot{r}\dot{\varphi}\sin 2\varphi\right.$
$\displaystyle-\left.2r^{2}\dot{\varphi}^{2}\cos
2\varphi-r^{2}\ddot{\varphi}\sin 2\varphi\right)\>,$ (36) $\displaystyle
h_{\times}=\frac{2\mu}{d}\left(\dot{r}^{2}\sin 2\varphi+r\ddot{r}\sin
2\varphi+4r\dot{r}\dot{\varphi}\cos 2\varphi\right.$
$\displaystyle-\left.2r^{2}\dot{\varphi}^{2}\sin
2\varphi+r^{2}\ddot{\varphi}\cos 2\varphi\right)\>.$ (37)

The same components for the GW signal induced by the tidal oscillations are

$h^{Q}_{+}=\frac{\sqrt{2}}{d}\ddot{Q}_{+}\>,\qquad
h^{Q}_{\times}=\frac{\sqrt{2}}{d}\ddot{Q}_{\times}\>,$ (38)

with the variables $Q_{+}$ and $Q_{\times}$ introduced in Appendix B.

Generally, the SNR can be calculated for any waveform using the Fourier-
spectrum of the measured strain, $\tilde{h}(f)$, assuming Wiener-filtering
(see e.g. Ref. Thrane  _et al._ (2014)):

$\mathrm{SNR}^{2}=4\int_{0}^{\infty}\frac{|\tilde{h}(f)|^{2}}{S_{h}(f)}\mathrm{d}f\>,$
(39)

where $S_{h}(f)$ is the one-sided spectral density of the detector
noise101010We used the updated Advanced LIGO sensitivity curve:
https://dcc.ligo.org/T1800044-v5/public and the complicated angular depdence
is incorporated in $\tilde{h}(f)$ For L-shaped detectors, like LIGO, Virgo and
KAGRA, the position and orientation averaged SNR is reduced by a factor of
$2/5$ compared to the SNR in the optimal (face-on) orientation and sky
position (see e.g. Ref. Dalal  _et al._ (2006)):

$\mathrm{SNR}_{\mathrm{ave}}^{\mathrm{LIGO}}=\frac{2}{5}\mathrm{SNR}_{\mathrm{opt}}^{\mathrm{LIGO}}\>,$
(40)

while for the Einstein Telescope this factor is reduced by $\sqrt{3}/2$ due to
the $60^{\circ}$ opening angle between the interferometer arms, but increased
by a factor of $\sqrt{3}$ due to essentially acting as 3 independent
detectors. Hence the above relation becomes Regimbau  _et al._ (2012):

$\mathrm{SNR}_{\mathrm{ave}}^{\mathrm{ET}}=\frac{3}{5}\mathrm{SNR}_{\mathrm{opt}}^{\mathrm{ET}}\>.$
(41)

Then the detection range of a binary system can be defined as the distance
where the averaged SNR is 8:

$d_{\rm
L,det}\,[\mathrm{Mpc}]=\frac{1}{8}\mathrm{SNR}_{\mathrm{ave}}(d=1~{}\mathrm{Mpc})\>.$
(42)

Calculating the detection range for directly measuring NS oscillations, we
obtain sub-megaparsec distances for a single aLIGO detector and a few Mpcs for
the Einstein Telescope, using design sensitivities, even for the highly
excited case of $a_{0}=600$ km, $\Lambda=700$.

In order to assess the SNR of the phase shift caused by dynamical tides let us
proceed similarly to Ref. Kocsis  _et al._ (2011) and calculate the SNR for
the subtracted spectrum:

$\delta\tilde{h}(f)=\tilde{h}_{\mathrm{dyn}}(f)-\tilde{h}_{\mathrm{ad}}(f)\>,$
(43)

where $\tilde{h}_{\mathrm{dyn}}(f)$ is the GW spectrum of the orbital
evolution with GW radiation reaction and dynamical tides, and
$\tilde{h}_{\mathrm{ad}}(f)$ is the same for the orbital evolution with GW
radiation reaction and only adiabatic tides. To mimic matched filtering we fix
the phase of the two signals at their highest peak (technically at the last
pericenter passage before the simulations are terminated at $r=28$ km, which
conservatively corresponds to a direct collision), conservatively truncate the
signals beyond this peak and calculate their Fourier-spectra.111111Note that
the difference of the two signals vanishes at the truncation boundary and we
have verified that the truncation does not introduce a significant spurious
powerlaw spectral tail. We do not account for the GW signal beyond this point
since there the assumed stellar oscillation and dissipation model may break,
the NSs may get tidally disrupted or be subject to other non-perturbative
effects. Since these limitations may not necessarily only apply strictly at
the last orbit, we examine the cumulative signal to noise ratio as a function
of the pericenter distance $r_{\rm p}$. Note that $r_{\rm p}$ decreases
monotonically during the time-evolution.

Figure 9: The cumulative signal-to-noise ratio (SNR) of dynamical tides for
sources at distance 40 Mpc initialized with $e_{0}=0.9$ and $a_{0}$ as
labelled until reaching pericenter distance $r_{\rm p}$ shown on the
horizontal axis, for different initial conditions and dimensionless tidal
deformabilities. Note that direct collision corresponds to
$r_{\mathrm{p,coll}}=28$ km. Solid lines represent the average of individual
SNRs, while the coloured patches represent the $68\%$ confidence intervals.
Circles correspond to pericenter distances with characteristic GW frequencies
(see Eq. 33) $500$ Hz (empty circle) and $1000$ Hz (filled circle). The
horizontal dashed lines show the SNR of the GW signals of orbital evolutions
without tidal interactions for the two different initial conditions,
respectively.

Fig. 9 shows the orientation and sky position averaged cumulative SNR for
measuring the tidal effect with a single aLIGO detector as a function of the
pericenter distance $r_{\mathrm{p}}$, calculated from the difference spectrum
$\delta\tilde{h}$. The distance of the binaries was chosen to be the same as
that of the source of GW170817, which was estimated to be $\sim 40$ Mpc Abbott
_et al._ (2019). The SNR without tidal interactions are also indicated by the
horizontal dashed lines for the same initial conditions. As expected from the
phase shifts shown in the previous figures, the SNR for detecting dynamical
tides can become significant. For the two configurations with $a_{0}=600$ km
($r_{\rm p0}=60$ km) it exceeds the SNR of the unperturbed orbit. At this
point the tidal effects obviously can no longer cosidered to be small
perturbations to the original GW-driven orbit plus adiabatic tides, as it is
assumed otherwise. The figure shows that the SNR accumulates to detectable
levels already at relatively large periapsis separations, where the post-
Newtonian treatment and the treatment of the tidal effects may be expected to
be valid. The phase shifts for the orbits with $a_{0}=800$ km ($r_{\rm p0}=80$
km), which experience weaker tidal interactions, also become significant at
small separations. Also note that these results only correspond to a single
aLIGO detector and considering the whole GW detector network would increase
these by a factor of $1.5-2$. This means that it might be feasible to observe
the effect of dynamical tides already with the current network of GW
detectors. The SNRs (or equivalently the detection ranges) of dynamical tides
increase by more than an order of magnitude when we calculate with the design
sensitivity of the Einstein Telescope. The detection ranges for the two
detectors and the different binary configurations at two different
characteristic GW frequencies are listed in Table 2.

$a_{0}$ [km], $\Lambda$ Detection range [Mpc] 800, 700 800, 600* 600, 500 600,
700 $d_{\rm L,det}^{\mathrm{LIGO}}$ ($\Lambda=0$) 65 65 46 46 $d_{\rm
L,det}^{\mathrm{ET}}$ ($\Lambda=0$) 817 817 629 629 $d_{\rm
L,det}^{\mathrm{LIGO}}$ ($f_{\mathrm{GW}}=500$ Hz) $3.1^{+0.4}_{-0.5}$
$4.8^{+1.0}_{-0.9}$ $25^{+13}_{-15}$ $43^{+12}_{-14}$ $d_{\rm
L,det}^{\mathrm{ET}}$ ($f_{\mathrm{GW}}=500$ Hz) $50.9^{+8.5}_{-8.0}$
$82^{+13}_{-17}$ $390^{+260}_{-280}$ $760^{+220}_{-220}$ $d_{\rm
L,det}^{\mathrm{LIGO}}$ ($f_{\mathrm{GW}}=1$ kHz) $17.2^{+1.9}_{-3.3}$
$26.9^{+5.6}_{-5.3}$ $46^{+14}_{-17}$ $77^{+10}_{-12}$ $d_{\rm
L,det}^{\mathrm{ET}}$ ($f_{\mathrm{GW}}=1$ kHz) $280^{+34}_{-49}$
$440^{+77}_{-86}$ $770^{+280}_{-280}$ $1270^{+180}_{-180}$

Table 2: Detection ranges for different binary configurations at two different
GW frequencies, calculated by using the design sensitivity curves of a single
aLIGO detector as well as the Einstein Telescope. The detection ranges for
$\Lambda=0$ were calculated from the orbit-averaged results O’Leary  _et al._
(2009).

It is also important to point out that due to the strong dependence of $\delta
M$ on the f-mode frequency, the phase shift and therefore the SNR of a binary
system with $\Lambda=600$ greatly increases if we set the f-mode frequency to
$1.5$ kHz instead of $\sim 1.66$ kHz predicted by the quasi-universal relation
between $\Lambda$ and $f_{0}$. They even increase above the values obtained
for the binary with $\Lambda=700$ and $f_{0}\approx 1.61$ kHz. Due to this
strong dependence on the f-mode frequency the measurement of $\delta M$ could
be used to test the quasi-universal relation between $\Lambda$ and $f_{0}$.

Figure 10: Orientation and sky location averaged signal-to-noise ratio of the
subtracted spectrum $\delta h$ (see Eq. 43), calculated from GW strains
truncated at a final periapsis separation of $r_{\mathrm{p}}=40$ km, as a
function of initial orbital parameters $r_{\mathrm{p0}}$ and $e_{0}$. The two
panels correspond to $\Lambda=500$ (left) and $\Lambda=700$ (right) with
corresponding f-mode frequencies of $1.73$ kHz and $1.61$ kHz, respectively.
Both panels show values averaged over several, slightly different initial
conditions in order to remove stochasticity. The source distance is set to
$40$ Mpc, corresponding to the estimated distance of GW170817. The two contour
lines correspond to SNRs of $8$ (solid line) and $2/5\cdot 8$ (dashed line),
indicating the detectability threshold of a single aLIGO detector with random
orientation, and one with optimal orientation, respectively.

Figure 10 shows the orientation and sky location averaged SNR of the
subtracted spectrum $\delta h$ (Eq. 43), calculated from GW strains truncated
at $r_{\mathrm{p}}=40$ km, as a function of initial orbital parameters
$r_{\mathrm{p0}}$ and $e_{0}$. The distance of the binary is set to $40$ Mpc.
Fig. 10 shows that the SNR is greatly increased in the highly eccentric case,
compared to the circular or moderately eccentric case (the SNR is increased by
a factor of $\sim 2$ between $e_{0}=0.7$ and $e_{0}=0.98$). While increasing
the initial pericenter distance quickly suppresses the effect of dynamical
tides, the SNR can still reach the detectability threshold (SNR$\sim 8$) for
binaries with $r_{\mathrm{p0}}\lesssim 68$ km ($r_{\mathrm{p0}}\lesssim 76$
km) and a tidal deformability parameter of $\Lambda=500$ ($700$), for a single
optimally oriented aLIGO detector.

However we caution that the assumed detectability limit of SNR$=8$ is usually
assumed for GW sources with available matched filtering templates, but it is
unclear whether a matched filtering search is feasible for NS binaries with
dynamical tides. The detection and measurement of NSs with dynamical tides may
be more difficult.

###  III.6 Implications for the equation of state

While dynamical tides affect the evolution of circular binaries as well, they
only become significant very close to the merger Hinderer  _et al._ (2016);
Schmidt and Hinderer (2019). Therefore, the tidal phase shift is mostly
determined by adiabatic tides, which only depend on the mass-weighted tidal
deformability:

$\tilde{\Lambda}=\frac{16}{13}\frac{(m_{1}+12m_{2})m_{1}^{4}\Lambda_{1}+(m_{2}+12m_{1})m_{2}^{4}\Lambda_{2}}{(m_{1}+m_{2})^{5}}\>.$
(44)

Hence, the measurement of adiabatic tides only provide information about this
specific combination of binary parameters, while $\Lambda_{1}$ and
$\Lambda_{2}$ could only be constrained by imposing additional constraints on
the EoS. On the other hand, dynamical tides also carry information about the
f-mode frequency of tidal oscillations. This could resolve the degeneracy
between the binary parameters and could be used together with quasi-universal
relations to determine $\Lambda_{1}$ and $\Lambda_{2}$ independently, thus
also providing information about the radii, $R_{1}$ and $R_{2}$ of the two
NSs.

An interesting open question about the EoS of dense strongly interacting
matter is whether it contains a strong first-order phase transition at
densities available for NSs Annala  _et al._ (2020); Takatsy  _et al._ (2023).
In case it does, a new stable branch of compact stars could emerge, containing
some form of exotic matter (e.g. quark matter) inside their cores. The two
branches are then separated by an unstable branch, characterized by decreasing
NS mass with increasing central energy density. Thus, two compact stars having
the same mass but positioned on different stable branches, therefore having
different radii, can exist, referred to as _twin stars_ Glendenning and
Kettner (2000); Alford and Han (2016); Christian and Schaffner-Bielich (2022).

Figure 11: The phase shift at different peak GW frequencies for twin stars
(see text) with masses of $1.4~{}\mathrm{M}_{\odot}$ but different radii,
resulting in different tidal deformabilities and f-mode frequencies. The phase
shifts were obtained from Eq. (30), where the input orbital parameters were
extracted from our numerical simulations run in the adiabatic tide limit, i.e.
with $f_{0}=f_{\mathrm{ad}}\gg f_{\mathrm{orb}}$. The mass-weighted tidal
parameter is kept constant ($\tilde{\Lambda}=500$), while different colors
correspond to various differences between the radii of the twin stars. The two
different markers correspond to different initial orbital configurations,
where $e_{0}=0.9$ for each case. Figure 12: The phase shift at different peak
GW frequencies for binary NSs with a $1.4\mathrm{M}_{\odot}$ and a
$2\mathrm{M}_{\odot}$ component. The phase shifts were obtained from Eq. (30),
where the input orbital parameters were extracted from our numerical
simulations run in the adiabatic tide limit, i.e. with
$f_{0}=f_{\mathrm{ad}}\gg f_{\mathrm{orb}}$. The mass-weighted tidal parameter
is kept constant ($\tilde{\Lambda}\approx 194$, calculated from the equal-
radius case with $R_{\mathrm{eq}}=12.5$ km), while different colors correspond
to different slopes for the $M-R$ relation resulting in various differences
between the radii of the two NSs ($R(2.0~{}\mathrm{M}_{\odot})\geq
R(1.4~{}\mathrm{M}_{\odot})$). The two different markers correspond to
different initial orbital configurations, where $e_{0}=0.9$ for each case.

With only the mass-weighted tidal deformability, Eq. (44), we cannot
distinguish between a binary system consisting of twin stars and one comprised
of identical NSs from a single stable branch. On the other hand, measuring the
effect of dynamical tides might enable us to break the degeneracy and to
differentiate between these two possibilities.

Fig. 11 shows the phase shift caused by dynamical tides as a function of
pericenter separation for different $a_{0}$ and binaries with the same
component mass $m_{1}=m_{2}=1.4\mathrm{M}_{\odot}$ but different NS radii with
$\Delta R=R_{1}-R_{2}$ as shown in the legend. Lines with different colors
correspond to binary NSs having the same mass-weighted tidal deformability but
different values for the individual tidal deformabilities, and hence different
radii. The mass of both NSs was set to $1.4~{}\mathrm{M}_{\odot}$, while we
also utilized the quasi-universal relation of Ref. Godzieba  _et al._ (2021)
between the compactness $C$ and $\Lambda$, in order to determine the
individual radii of NSs given a fixed value for $\tilde{\Lambda}$. The f-mode
frequencies were determined from the quasi-universal relation of Ref. Pradhan
_et al._ (2022) between $\Lambda(1.4~{}\mathrm{M}_{\odot})$ and
$f_{0}(1.4~{}\mathrm{M}_{\odot})$. The value of the mass-weighted tidal
deformability was set to $\tilde{\Lambda}=500$, which corresponds to a radius
of $R_{\mathrm{eq}}\approx 12.2$ km for two NSs with identical radii.
Increasing the value of $\Delta R=R_{1}-R_{2}$ increases the effect of
dynamical tides. For larger values of $\Delta R$ the respective phase shift
contribution becomes dominated by the NS with the larger radius, hence the one
with the larger tidal deformability and lower f-mode frequency. A radius
difference of $2.5$ km increases the phase shift by a factor of $1.5-2$ for
$a_{0}=600$ km and $e_{0}=0.9$, and by a factor of $2.5-3$ for a binary with
$a_{0}=800$ km and $e_{0}=0.9$. Smaller radius differences are more difficult
to identify using the GW phase shift measurement; a radius difference of $1$
km results in a $15-20\%$ and a $25-45\%$ increase in the phase shift relative
to that of a binary with identical stellar radii, respectively. Note that this
difference might still be possible to detect, provided a signal with
sufficiently large SNR.

Similarly, the slope of the mass–radius curve could also be measured through
observations of non-equal mass NS binaries. Similar to Fig. 11, Fig. 12 shows
the phase shift for a binary with a $1.4~{}\mathrm{M}_{\odot}$ and a
$2~{}\mathrm{M}_{\odot}$ component. Here, the f-mode frequency and the tidal
deformability of the $2~{}\mathrm{M}_{\odot}$ NS was determined from the
relations in Refs. Pradhan  _et al._ (2023); Godzieba  _et al._ (2021) between
the compactness and $f_{0}$, and the compactness and $\Lambda$, respectively.
The radius of the two binary components were chosen so that
$R_{2.0\mathrm{M}_{\odot}}\geq R_{1.4\mathrm{M}_{\odot}}$, consistently with
recent EoS constraints Miller  _et al._ (2021); Takatsy  _et al._ (2023);
Koehn  _et al._ (2024). In the equal-radius case the radii were set to be
$R_{\mathrm{eq}}=12.5$ km, while in all other cases the radii were determined
in a way to keep the mass-weighted tidal deformability fixed
($\tilde{\Lambda}\approx 194$). An interesting trend can be observed in Fig.
12, since the phase shift first decreases with increasing $\Delta R$, then it
starts to increase, resulting in the radius differences of $1$ km and $2.5$ km
having almost the same effect. This is due to the increase of
$f_{0}(1.4~{}\mathrm{M}_{\odot})$ first dominating and later the decrease of
$f_{0}(2~{}\mathrm{M}_{\odot})$ taking over. The largest difference in the
phase shift is $25-30\%$ and $30-40\%$ for the two binary configurations with
$a_{0}=600$ km and $800$ km, respectively.

##  IV Conclusion

While dynamical tides only become relevant during the last couple of orbits
for circular inspirals, orbital eccentricity can increase these effects
considerably due to the fact that the periapsis distance may be much smaller
for a much larger number of periods for highly eccentric orbits. In the
eccentric case the spectrum of the driving tidal forces extend to higher
orbital harmonics. In this study, we have utilized a post-Newtonian model with
realistic NS parameters to investigate the effect of dynamical tides in
eccentric NS binary systems.

By calculating the energy transferred to tidal oscillations relative to the
amount of energy radiated through GWs we have found a strong, exponential
dependence on the dimensionless pericenter distance, with dynamical tides
being strongly suppressed for binaries formed with large initial pericenter
distances (Figs. 2–3). We have shown that the energy transferred to tidal
oscillations by dynamical tides are below $10^{-3}$ relative to the energy
loss due to GWs for NS-BH binaries in most cases (Fig. 3). For NS-NS binaries
this ratio increases beyond $10\%$ for periapsis distance of $r_{\rm
p}=28~{}\mathrm{km}-40~{}\mathrm{km}$ depending on eccentricity.

While GWs generated by tidal oscillations driven by dynamical tides typically
cannot be detected directly, their effect on the orbital evolution may be
observed through a phase shift in the GW signal of the inspiral. Since the
decay timescale of the tidal oscillations is much larger than the orbital
period, the amplitude of the oscillations varies stochastically from one orbit
to another, causing stochastic variations in the orbital phase shift as well.
We have demonstrated that the expectation value of these phase shifts are
accurately modelled analytically for the Newtonian model complemented with the
leading-order radiation terms (Fig. 4). We have also shown that the phase
shifts calculated from the analytical formula agree well with numerical
calculations using our pseudo-Newtonian model (Fig. 7). We have also found
that including the 1PN and 2PN terms in the equations of motion does not
change the frequency-dependent evolution of the phase shift significantly,
although the orbital frequency before merger is reduced, and therefore so is
the phase shift before merger as well (Fig. 8).

We have shown that the SNR associated with the GW phase shift caused by the
dynamical tides can be significant for NSs with an initial eccentricity of 0.9
with periapsis distances of 60-80km at similar distances to the detected
source of GW170817 (40 Mpc) even with a single aLIGO detector (Fig. 9), and
the detection range of the dynamical tides for ET lies between a few hundred
Mpc to $\sim 1$ Gpc for binary NSs with initial pericenter distances
$r_{\mathrm{p0}}\lesssim 80$ km and tidal deformabilities in the reasonable
range $\Lambda\sim 400-700$ (see Table 2). Below this range of tidal
deformability, the f-mode frequency is larger which suppresses dynamical tides
making these phase shifts more difficult to observe. These eccentric sources
form in the LIGO frequency band and have a signal to noise ratio well above 8
for aLIGO and may be expected to be detactable if neglecting the perturbations
of dynamical tides Gondán and Kocsis (2019, 2021). However dynamical tides
imprint stochastic variations on the orbital phase, particularly at periapsis
passages, which may make it more difficult to detect these sources with
matched filtering.

Currently it is not well understood how likely it may be to detect highly
eccentric neutron star binaries in the Universe. These systems may form during
GW captures in dense, high velocity dispersion environments O’Leary  _et al._
(2009); Gondán  _et al._ (2018a); Gondán and Kocsis (2021) or GW captures
during binary-single or binary-binary interactions Samsing  _et al._ (2014);
Samsing (2018); Zevin  _et al._ (2019) especially in flattened configurations
such as in active galactic nuclei (AGN) accretion disks Samsing  _et al._
(2022). The initial periapsis distribution of such highly eccentric binaries
with $e_{0}\sim 0.99$ peaks at
$r_{p0}=100\mathrm{km}(m/1.4\mathrm{M}_{\odot})(v/500\mathrm{km/s})^{-4/7}$,
where for GW captures $v$ is the escape velocity from the system and for
binary-single interactions it is the initial orbital velocity of the binary
Gondán and Kocsis (2021); Kocsis (2022). In comparison, we found that the
signal-to-noise ratio of detecting dynamical tides in the GW signal is
significant below an initial periapsis of $r_{p0}\lesssim 70$km (Fig. 10).
Thus the existence of highly eccentric NS binaries with strong dynamical tides
hinges on the existence of dense clusters abundant in neutron stars with high
escape velocities or systems with very tight binaries undergoing scattering
interactions. The rate of the former is theoretically expected to be low
(Tsang, 2013); the rates in AGN may possibly be more substantial (Tagawa  _et
al._ , 2021; Samsing  _et al._ , 2022).

NS seismology, specifically the measurement of the f-mode frequency of NSs in
this case, can be immensely useful in characterising the EoS of dense strongly
interacting matter, since it can serve as independent and therefore
complementary information to measurements of masses, radii and tidal
deformability, all characterising the structure of NSs in different ways.
While these properties can be connected through quasi-universal relations,
deviations from these relations are possible, especially for equations of
state with a strong first-order phase transition. We have shown that
deviations from these quasi-universal relations can cause significant
differences in the GW phase shift due to the strong dependence of these
effects on the f-mode frequency of the NSs, which in turn allows us to measure
these deviations.

Another interesting possibility provided by f-mode frequency measurements is
the ability of telling equal-radius and twin-star binaries apart (i.e.
binaries consisting of NSs with the same mass but different radii). We have
shown that while there is a complete degeneracy between equal-mass binary NSs
with the same or different component radii with respect to the the mass-
weighted tidal deformability, characterizing the magnitude of adiabatic tides,
this degeneracy is broken by dynamical tides and the two branches can be
distinguished by measuring the corresponding phase shift (Fig. 11). We have
also shown that dynamical tides might also allow us to measure the slope of
the mass–radius relation with unequal-mass binaries, although this effect was
found to be much smaller (Fig. 12).

While we do not expect major deviations from our leading-order calculations
with respect to the magnitude of the phase shift caused by dynamical tides, an
improvement on our assumptions would be necessary in order to be able to build
accurate waveform models. First of all, the effect of higher-order post-
Newtonian terms modify the tidal excitation spectrum near the merger, which
have been neglected in our analytical model. In addition, while back-reaction
forces from GW radiation are expected to be accurate for the orbital
evolution, as well as for the freely-oscillating part of the tidal evolution,
cross-terms can become important near the pericenter passages, hence these
interactions should also be revisited. We note however, that in our
calculations the quadrupole moments $Q_{m}$ always remained smaller than $\sim
10\%$ of the unperturbed moment of inertia of the NSs, even for the strongest
tides, while they were below $\sim 1\%$ a few orbits before the merger.
Therefore we consider our results to be reliable from the point of view of
parametric instabilities or other non-perturbative effects that can lead to
the saturation of tides. While the spin frequency of inspiralling NSs is
expected to be low and usually aligned with the orbital angular momentum, due
to the formation and the reduced inspiral timescale of eccentric binaries, the
existence of binary NSs with considerable anti-aligned spins might be
possible. In this case, the spin of the NS can effectively reduce the
resonance frequency of the f-mode oscillations, leading to even stronger
dynamical tides. This scenario should also be investigated in a follow-up
study.

##  Appendix A The gravitational potential

Let us define the external and induced Newtonian multipolar moments by first
expanding the external gravitational potential (which is non-singular in the
origin) in a Taylor-series:

$\displaystyle\Phi_{\rm ext}({\bf r})$ $\displaystyle=\Phi_{\rm
ext}(0)+{r}_{i}\partial_{i}\Phi_{\rm ext}\bigr{\rvert}_{0}$
$\displaystyle+\frac{1}{2}{r}_{i}{r}_{j}\partial_{i}\partial_{j}\Phi_{\rm
ext}\bigr{\rvert}_{0}+\ldots\>.$ (45)

Since the second term simply generates an overall translation, it can be
transformed out. Then, let us write equivalently:

$\displaystyle\Phi_{\rm ext}({\bf r})=\Phi_{\rm
ext}(0)+\frac{1}{2}{r}_{i}{r}_{j}\mathcal{E}_{ij}+\ldots\>.$ (46)

Similarly, one can expand the potential generated by the massive object under
the influence of the external potential (which can be singular in the origin
but vanishes at infinity) in a multipole-series:

$\Phi_{\rm ind}({\bf
r})=-\left(\frac{m}{r}+\frac{3}{2}Q_{ij}\frac{{r}_{i}{r}_{j}}{r^{5}}+\ldots\right)\quad,$
(47)

where the dipole moment has been removed by choosing the center-of-mass frame
as the reference frame. Using this convention for the quadrupole moment and
Eq. (46) the interaction energy becomes

$U_{\mathrm{int}}=\frac{1}{2}Q_{ij}\mathcal{E}_{ij}\>,$ (48)

where we have used the trace-free property of $\mathcal{E}_{ij}$.

##  Appendix B Equations of motion

For our numerical simulations we find it convenient to introduce the variables
$Q_{+}$ and $Q_{\times}$ defined by

$Q_{+}=\frac{Q_{-2}+Q_{2}}{\sqrt{2}}\>,\quad
Q_{\times}=\frac{Q_{-2}-Q_{2}}{i\sqrt{2}}\>.$ (49)

This way our equations of motion will become purely real. Then, considering
only the tidal deformation of object $1$ (similar equations hold for object
$2$), the complete equations of motion including 1PN, 2PN and the 2.5PN
leading order radiation terms for the orbital evolution are

$\displaystyle{\ddot{Q}}_{+}+2\gamma_{0}{\dot{Q}}_{+}+\omega_{0}^{2}Q_{+}$
$\displaystyle=$
$\displaystyle{3\over\sqrt{2}}\lambda_{1}\omega_{0}^{2}{m_{2}\over r^{3}}\cos
2\varphi\ ,$ (50)
$\displaystyle{\ddot{Q}}_{\times}+2\gamma_{0}{\dot{Q}}_{\times}+\omega_{0}^{2}Q_{\times}$
$\displaystyle=$
$\displaystyle{3\over\sqrt{2}}\lambda_{1}\omega_{0}^{2}{m_{2}\over r^{3}}\sin
2\varphi\ ,$ (51)
$\displaystyle{\ddot{Q}}_{0}+2\gamma_{0}{\dot{Q}}_{0}+\omega_{0}^{2}Q_{0}$
$\displaystyle=$ $\displaystyle\sqrt{3\over
2}\lambda_{1}\omega_{0}^{2}{m_{2}\over r^{3}}\ ,$ (52)

for the tidal variables of object $1$ and

$\displaystyle{\ddot{r}}$ $\displaystyle=$ $\displaystyle
r{\dot{\varphi}}^{2}-\frac{9}{2\sqrt{2}}\frac{m_{2}}{\mu r^{4}}\left(Q_{+}\cos
2\varphi+Q_{\times}\sin 2\varphi+\frac{Q_{0}}{\sqrt{3}}\right)$ (53)
$\displaystyle-$ $\displaystyle\frac{m_{\mathrm{tot}}}{r^{2}}\left(1+A_{\rm
PN}+A_{5/2}+B_{\rm PN}\dot{r}+B_{5/2}\dot{r}\right)\ ,$
$\displaystyle{\ddot{\varphi}}$ $\displaystyle=$
$\displaystyle-2\frac{\dot{\varphi}\dot{r}}{r}-{3\over\sqrt{2}}\frac{m_{2}}{\mu
r^{5}}\left(Q_{+}\sin 2\varphi-Q_{\times}\cos 2\varphi\right)$ (54)
$\displaystyle-$
$\displaystyle\frac{m_{\mathrm{tot}}}{r^{2}}\dot{\varphi}\left(B_{\rm
PN}+B_{5/2}\right)\ ,$

for the orbital parameters, with the coefficients Lincoln and Will (1990)

$\displaystyle A_{5/2}$
$\displaystyle=-\frac{8\mu}{5r}\dot{r}\left(18v^{2}+\frac{2m_{\mathrm{tot}}}{3r}-25\dot{r}^{2}\right),$
(55) $\displaystyle B_{5/2}$ $\displaystyle=\
\frac{8\mu}{5r}\left(6v^{2}-\frac{2m_{\mathrm{tot}}}{r}-15\dot{r}^{2}\right),$
(56) $\displaystyle A_{\rm PN}$ $\displaystyle=\
(1+3\eta)v^{2}-2(2+\eta)\frac{m_{\mathrm{tot}}}{r}-\frac{3}{2}\eta\dot{r}^{2}$
$\displaystyle+\frac{3}{4}(12+29\eta)\left(\frac{m_{\mathrm{tot}}}{r}\right)^{2}+\eta(3-4\eta)v^{4}$
$\displaystyle+\frac{15}{8}\eta(1-3\eta)\dot{r}^{4}-\frac{3}{2}\eta(3-4\eta)v^{2}\dot{r}^{2}$
$\displaystyle-\frac{1}{2}\eta(13-4\eta)\frac{m_{\mathrm{tot}}}{r}v^{2}$
$\displaystyle-(2+25\eta+2\eta^{2})\frac{m_{\mathrm{tot}}}{r}\dot{r}^{2},$
(57) $\displaystyle B_{\rm PN}$
$\displaystyle=-2(2-\eta)\dot{r}-\frac{1}{2}\dot{r}\left[\vphantom{\frac{m_{\mathrm{tot}}}{r}}\eta(15+4\eta)v^{2}\right.$
$\displaystyle\left.-(4+41\eta+8\eta^{2})\frac{m_{\mathrm{tot}}}{r}-3\eta(3+2\eta)\dot{r}^{2}\right],$
(58)

where $v^{2}=\dot{r}^{2}+r^{2}\dot{\varphi}^{2}$,
$m_{\mathrm{tot}}=m_{1}+m_{2}$ and $\mu=m_{1}m_{2}/m_{\mathrm{tot}}$ are the
total and reduced masses, and $\eta=\mu/m_{\mathrm{tot}}$ is the symmetric
mass ratio. Note that the 1PN and 2PN terms are neglected (i.e.
$A_{\mathrm{PN}}=B_{\mathrm{PN}}=0$) during most of our analysis, and only
taken into account at the end of Sec. III.3. We solve these equations using a
fourth-order Runge-Kutta differential equation solver with adaptive timesteps.

Fig. 13 shows the typical time evolution of the quadrupole moment. The
stochastic motion of the oscillation amplitude from orbit to orbit can be
clearly observed. Fig. 14 shows the analytic spectrum calculated from Eq. 13,
at different stages of the binary evolution.

Figure 13: The time evolution of the quadrupole moment Re($Q_{-2}$) for a
binary NS with initial orbital parameters $a_{0}=600$ km and $e_{0}$, and
tidal parameters $\Lambda=400$ and $f_{0}=1.8$ kHz. Figure 14: The analytic
spectrum of the excited quadrupole oscillation $Q_{-2}$ at different stages of
the binary evolution. The initial orbital parameters are $a_{0}=600$ km and
$e_{0}$, while tidal parameters were chosen as $\Lambda=400$ and $f_{0}=1.8$
kHz. The datapoints in the discrete spectra are connected to guide the eye.
The inset shows the leading-order evolution of the binary parameters without
tidal interactions, calculated from Eq. 14.

When $A_{\mathrm{PN}}=B_{\mathrm{PN}}=0$, the Keplerian orbital parameters are
calculated at the apocenter from the specific orbital energy and angular
momentum:

$\displaystyle
E_{\mathrm{orb}}=\frac{v_{\mathrm{a}}^{2}}{2}-\frac{m_{\mathrm{tot}}}{r_{\mathrm{a}}}\>$
$\displaystyle,\quad
L_{\mathrm{orb}}=r_{\mathrm{a}}^{2}\dot{\varphi}_{\mathrm{a}}\>,$
$\displaystyle
e=\sqrt{1+\frac{2E_{\mathrm{orb}}L_{\mathrm{orb}}^{2}}{m_{\mathrm{tot}}^{2}}}\>$
$\displaystyle,\quad a=\frac{m_{\mathrm{tot}}}{2E_{\mathrm{orb}}}\>,$ (59)

where the “a” subscript refers to the apocenter. During the calculation of the
mean anomaly $M$, for small $\delta M$ differences, due to numerical
inaccuracies in determining the apocenter, we find that more accurate results
can be obtained by directly calculating $M$ from the true anomaly $\varphi$,
using

$\displaystyle M=$
$\displaystyle\mathrm{atan2}\left(-\sqrt{1-e^{2}}\sin\varphi,-e-\cos\varphi\right)+$
$\displaystyle+\pi-e\frac{\sqrt{1-e^{2}}\sin\varphi}{1+e\cos\varphi}\>.$ (60)

Including the 1PN and 2PN terms introduces additional effects, such as non-
negligible apsidal precession, which modifies Eq. (60). We circumvent the
complications by directly calculating the direction of the periastron using
the Runge-Lenz vector, defined as

$\mathbf{A}=\mathbf{p}\times\mathbf{L}-\mu^{2}m_{\mathrm{tot}}\frac{\mathbf{r}}{r}\>,$
(61)

and adjusting for its direction when calculating the true anomaly and the mean
anomaly afterwards. The components of the Runge-Lenz vector using polar
coordinates are

$\displaystyle\mathbf{A}$
$\displaystyle=\mu^{2}r^{2}\dot{\varphi}\begin{pmatrix}\dot{r}\sin\varphi+r\dot{\varphi}\cos\varphi\\\
-\dot{r}\cos\varphi+r\dot{\varphi}\sin\varphi\end{pmatrix}$
$\displaystyle-\mu^{2}m_{\mathrm{tot}}\begin{pmatrix}\cos\varphi\\\
\sin\varphi\end{pmatrix}\>.$ (62)

The 2-body motion with 1PN and 2PN interactions included can be projected into
a quasi-Keplerian picture with modified orbital parameters, as well as orbital
frequency Blanchet (2014). In our analytic calculations we take into account
the modification of $\omega_{\mathrm{orb}}$ due to PN interactions to the 2PN
order as described in Ref. Blanchet (2014) (see Eq. (347a) in Section 10.2).

###### Acknowledgements.

J. T. and P. K. acknowledge support by the National Research, Development and
Innovation (NRDI) fund of Hungary, financed under the K 21 funding scheme,
Project No. K 138277. This work was supported by the UK Science and Technology
Facilities Council Grant Number ST/W000903/1 (to B.K.).

## References

  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2023) R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13, 041039 (2023), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017b), arXiv:1710.05833 [astro-ph.HE] .
  * Abbott  _et al._ (2017c) B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Astrophys. J. Lett. 848, L13 (2017c), arXiv:1710.05834 [astro-ph.HE] .
  * Abbott  _et al._ (2017d) B. P. Abbott _et al._ (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), Nature 551, 85 (2017d), arXiv:1710.05835 [astro-ph.CO] .
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Kasen  _et al._ (2017) D. Kasen, B. Metzger, J. Barnes, E. Quataert, and E. Ramirez-Ruiz, Nature 551, 80 (2017), arXiv:1710.05463 [astro-ph.HE] .
  * Annala  _et al._ (2018) E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018), arXiv:1711.02644 [astro-ph.HE] .
  * Most  _et al._ (2018) E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120, 261103 (2018), arXiv:1803.00549 [gr-qc] .
  * De  _et al._ (2018) S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Phys. Rev. Lett. 121, 091102 (2018), [Erratum: Phys.Rev.Lett. 121, 259902 (2018)], arXiv:1804.08583 [astro-ph.HE] .
  * Tews  _et al._ (2018a) I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98, 045804 (2018a), arXiv:1804.02783 [nucl-th] .
  * Akmal  _et al._ (1998) A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998), arXiv:nucl-th/9804027 .
  * Wiringa  _et al._ (1988) R. B. Wiringa, V. Fiks, and A. Fabrocini, Phys. Rev. C 38, 1010 (1988).
  * Huovinen and Petreczky (2010) P. Huovinen and P. Petreczky, Nucl. Phys. A 837, 26 (2010).
  * Bazavov  _et al._ (2012) A. Bazavov _et al._ (HotQCD), Phys. Rev. D 86, 034509 (2012).
  * Lynn  _et al._ (2016) J. E. Lynn, I. Tews, J. Carlson, S. Gandolfi, A. Gezerlis, K. E. Schmidt, and A. Schwenk, Phys. Rev. Lett. 116, 062501 (2016).
  * Tews  _et al._ (2018b) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018b).
  * Tews  _et al._ (2013) I. Tews, T. Krüger, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 110, 032504 (2013), arXiv:1206.0025 [nucl-th] .
  * Tews  _et al._ (2019) I. Tews, J. Margueron, and S. Reddy, Eur. Phys. J. A 55, 97 (2019).
  * Kurkela  _et al._ (2010) A. Kurkela, P. Romatschke, and A. Vuorinen, Phys. Rev. D 81, 105021 (2010), arXiv:0912.1856 [hep-ph] .
  * Mogliacci  _et al._ (2013) S. Mogliacci, J. O. Andersen, M. Strickland, N. Su, and A. Vuorinen, JHEP 12, 055 (2013).
  * Dietrich  _et al._ (2020) T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020), arXiv:2002.11355 [astro-ph.HE] .
  * Huth  _et al._ (2022) S. Huth _et al._ , Nature 606, 276 (2022), arXiv:2107.06229 [nucl-th] .
  * Marczenko  _et al._ (2023) M. Marczenko, L. McLerran, K. Redlich, and C. Sasaki, Phys. Rev. C 107, 025802 (2023), arXiv:2207.13059 [nucl-th] .
  * Takatsy  _et al._ (2023) J. Takatsy, P. Kovacs, G. Wolf, and J. Schaffner-Bielich, Phys. Rev. D 108, 043002 (2023), arXiv:2303.00013 [astro-ph.HE] .
  * Flanagan and Hinderer (2008) E. E. Flanagan and T. Hinderer, Phys. Rev. D 77, 021502 (2008), arXiv:0709.1915 [astro-ph] .
  * Love (1909) A. E. H. Love, Proceedings of the Royal Society of London Series A 82, 73 (1909).
  * Damour and Nagar (2009) T. Damour and A. Nagar, Phys. Rev. D 80, 084035 (2009), arXiv:0906.0096 [gr-qc] .
  * Hinderer  _et al._ (2016) T. Hinderer _et al._ , Phys. Rev. Lett. 116, 181101 (2016), arXiv:1602.00599 [gr-qc] .
  * Steinhoff  _et al._ (2016) J. Steinhoff, T. Hinderer, A. Buonanno, and A. Taracchini, Phys. Rev. D 94, 104028 (2016), arXiv:1608.01907 [gr-qc] .
  * Schmidt and Hinderer (2019) P. Schmidt and T. Hinderer, Phys. Rev. D 100, 021501 (2019), arXiv:1905.00818 [gr-qc] .
  * Andersson and Pnigouras (2021) N. Andersson and P. Pnigouras, Mon. Not. Roy. Astron. Soc. 503, 533 (2021), arXiv:1905.00012 [gr-qc] .
  * Shibata (1994) M. Shibata, Progress of Theoretical Physics 91, 871 (1994).
  * Lai (1994) D. Lai, Mon. Not. Roy. Astron. Soc. 270, 611 (1994), arXiv:astro-ph/9404062 .
  * Kokkotas and Schaefer (1995) K. D. Kokkotas and G. Schaefer, Mon. Not. Roy. Astron. Soc. 275, 301 (1995), arXiv:gr-qc/9502034 .
  * Ajith  _et al._ (2008) P. Ajith _et al._ , Phys. Rev. D 77, 104017 (2008), [Erratum: Phys.Rev.D 79, 129901 (2009)], arXiv:0710.2335 [gr-qc] .
  * Khan  _et al._ (2016) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Field  _et al._ (2014) S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, and M. Tiglio, Phys. Rev. X 4, 031006 (2014), arXiv:1308.3565 [gr-qc] .
  * Bohé  _et al._ (2017) A. Bohé _et al._ , Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Yang  _et al._ (2018) H. Yang, W. E. East, V. Paschalidis, F. Pretorius, and R. F. P. Mendes, Phys. Rev. D 98, 044007 (2018), arXiv:1806.00158 [gr-qc] .
  * Yang (2019) H. Yang, Phys. Rev. D 100, 064023 (2019), arXiv:1904.11089 [gr-qc] .
  * Vick and Lai (2019) M. Vick and D. Lai, Phys. Rev. D 100, 063001 (2019), arXiv:1906.08780 [astro-ph.HE] .
  * Peters (1964) P. C. Peters, Phys. Rev. 136, B1224 (1964).
  * O’Leary  _et al._ (2009) R. M. O’Leary, B. Kocsis, and A. Loeb, Mon. Not. Roy. Astron. Soc. 395, 2127 (2009), arXiv:0807.2638 [astro-ph] .
  * Kocsis  _et al._ (2011) B. Kocsis, N. Yunes, and A. Loeb, Phys. Rev. D 84, 024032 (2011), arXiv:1104.2322 [astro-ph.GA] .
  * Gondán  _et al._ (2018a) L. Gondán, B. Kocsis, P. Raffai, and Z. Frei, Astrophys. J. 860, 5 (2018a), arXiv:1711.09989 [astro-ph.HE] .
  * Takátsy  _et al._ (2019) J. Takátsy, B. Bécsy, and P. Raffai, Mon. Not. Roy. Astron. Soc. 486, 570 (2019), arXiv:1812.04012 [astro-ph.HE] .
  * Rasskazov and Kocsis (2019) A. Rasskazov and B. Kocsis, Astrophys. J. 881, 20 (2019), arXiv:1902.03242 [astro-ph.HE] .
  * Gondán and Kocsis (2021) L. Gondán and B. Kocsis, Mon. Not. Roy. Astron. Soc. 506, 1665 (2021), arXiv:2011.02507 [astro-ph.HE] .
  * Zevin  _et al._ (2021) M. Zevin, I. M. Romero-Shaw, K. Kremer, E. Thrane, and P. D. Lasky, Astrophys. J. Lett. 921, L43 (2021), arXiv:2106.09042 [astro-ph.HE] .
  * Samsing  _et al._ (2014) J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, Astrophys. J. 784, 71 (2014), arXiv:1308.2964 [astro-ph.HE] .
  * Samsing (2018) J. Samsing, Phys. Rev. D 97, 103014 (2018), arXiv:1711.07452 [astro-ph.HE] .
  * Zevin  _et al._ (2019) M. Zevin, J. Samsing, C. Rodriguez, C.-J. Haster, and E. Ramirez-Ruiz, Astrophys. J. 871, 91 (2019), arXiv:1810.00901 [astro-ph.HE] .
  * Arca-Sedda  _et al._ (2021) M. Arca-Sedda, G. Li, and B. Kocsis, Astron. Astrophys. 650, A189 (2021), arXiv:1805.06458 [astro-ph.HE] .
  * Wen (2003) L. Wen, Astrophys. J. 598, 419 (2003), arXiv:astro-ph/0211492 .
  * Antonini and Perets (2012) F. Antonini and H. B. Perets, Astrophys. J. 757, 27 (2012), arXiv:1203.2938 [astro-ph.GA] .
  * Antognini  _et al._ (2014) J. M. Antognini, B. J. Shappee, T. A. Thompson, and P. Amaro-Seoane, Mon. Not. Roy. Astron. Soc. 439, 1079 (2014), arXiv:1308.5682 [astro-ph.HE] .
  * Silsbee and Tremaine (2017) K. Silsbee and S. Tremaine, Astrophys. J. 836, 39 (2017), arXiv:1608.07642 [astro-ph.HE] .
  * Hoang  _et al._ (2018) B.-M. Hoang, S. Naoz, B. Kocsis, F. A. Rasio, and F. Dosopoulou, Astrophys. J. 856, 140 (2018), arXiv:1706.09896 [astro-ph.HE] .
  * Petrovich and Antonini (2017) C. Petrovich and F. Antonini, Astrophys. J. 846, 146 (2017), arXiv:1705.05848 [astro-ph.HE] .
  * Rodriguez and Antonini (2018) C. L. Rodriguez and F. Antonini, Astrophys. J. 863, 7 (2018), arXiv:1805.08212 [astro-ph.HE] .
  * Fragione and Kocsis (2019) G. Fragione and B. Kocsis, Mon. Not. Roy. Astron. Soc. 486, 4781 (2019), arXiv:1903.03112 [astro-ph.GA] .
  * Fragione  _et al._ (2019) G. Fragione, E. Grishin, N. W. C. Leigh, H. B. Perets, and R. Perna, Mon. Not. Roy. Astron. Soc. 488, 47 (2019), arXiv:1811.10627 [astro-ph.GA] .
  * Tagawa  _et al._ (2021) H. Tagawa, B. Kocsis, Z. Haiman, I. Bartos, K. Omukai, and J. Samsing, Astrophys. J. Lett. 907, L20 (2021), arXiv:2010.10526 [astro-ph.HE] .
  * Samsing  _et al._ (2022) J. Samsing, I. Bartos, D. J. D’Orazio, Z. Haiman, B. Kocsis, N. W. C. Leigh, B. Liu, M. E. Pessah, and H. Tagawa, Nature 603, 237 (2022), arXiv:2010.09765 [astro-ph.HE] .
  * Fabj and Samsing (2024) G. Fabj and J. Samsing, arXiv e-prints , arXiv:2402.16948 (2024), arXiv:2402.16948 [astro-ph.GA] .
  * Li  _et al._ (2023) J. Li, A. M. Dempsey, H. Li, D. Lai, and S. Li, Astrophys. J. Lett. 944, L42 (2023), arXiv:2211.10357 [astro-ph.HE] .
  * Rowan  _et al._ (2023) C. Rowan, T. Boekholt, B. Kocsis, and Z. Haiman, Mon. Not. Roy. Astron. Soc. 524, 2770 (2023), arXiv:2212.06133 [astro-ph.GA] .
  * Whitehead  _et al._ (2023) H. Whitehead, C. Rowan, T. Boekholt, and B. Kocsis, (2023), arXiv:2309.11561 [astro-ph.GA] .
  * Messick  _et al._ (2017) C. Messick _et al._ , Phys. Rev. D 95, 042001 (2017), arXiv:1604.04324 [astro-ph.IM] .
  * Davies  _et al._ (2020) G. S. Davies, T. Dent, M. Tápai, I. Harry, C. McIsaac, and A. H. Nitz, Phys. Rev. D 102, 022004 (2020), arXiv:2002.08291 [astro-ph.HE] .
  * Aubin  _et al._ (2021) F. Aubin _et al._ , Class. Quant. Grav. 38, 095004 (2021), arXiv:2012.11512 [gr-qc] .
  * Tiwari  _et al._ (2016) V. Tiwari _et al._ , Phys. Rev. D 93, 043007 (2016), arXiv:1511.09240 [gr-qc] .
  * Ramos-Buades  _et al._ (2020) A. Ramos-Buades, S. Tiwari, M. Haney, and S. Husa, Phys. Rev. D 102, 043005 (2020), arXiv:2005.14016 [gr-qc] .
  * Cornish  _et al._ (2021) N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse, Phys. Rev. D 103, 044006 (2021), arXiv:2011.09494 [gr-qc] .
  * Dálya  _et al._ (2021) G. Dálya, P. Raffai, and B. Bécsy, Class. Quant. Grav. 38, 065002 (2021), arXiv:2006.06256 [astro-ph.HE] .
  * Klimenko  _et al._ (2016) S. Klimenko _et al._ , Phys. Rev. D 93, 042004 (2016), arXiv:1511.05999 [gr-qc] .
  * Gondán  _et al._ (2018b) L. Gondán, B. Kocsis, P. Raffai, and Z. Frei, Astrophys. J. 855, 34 (2018b), arXiv:1705.10781 [astro-ph.HE] .
  * Gondán and Kocsis (2019) L. Gondán and B. Kocsis, Astrophys. J. 871, 178 (2019), arXiv:1809.00672 [astro-ph.HE] .
  * Romero-Shaw  _et al._ (2023) I. M. Romero-Shaw, D. Gerosa, and N. Loutrel, Mon. Not. Roy. Astron. Soc. 519, 5352 (2023), arXiv:2211.07528 [astro-ph.HE] .
  * Gamba  _et al._ (2023) R. Gamba, M. Breschi, G. Carullo, S. Albanesi, P. Rettegno, S. Bernuzzi, and A. Nagar, Nature Astron. 7, 11 (2023), arXiv:2106.05575 [gr-qc] .
  * Gayathri  _et al._ (2022) V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczepanczyk, I. Bartos, M. Campanelli, S. Klimenko, C. O. Lousto, and R. O’Shaughnessy, Nature Astron. 6, 344 (2022), arXiv:2009.05461 [astro-ph.HE] .
  * Romero-Shaw  _et al._ (2022) I. M. Romero-Shaw, P. D. Lasky, and E. Thrane, Astrophys. J. 940, 171 (2022), arXiv:2206.14695 [astro-ph.HE] .
  * Gupte  _et al._ (2024) N. Gupte _et al._ , (2024), arXiv:2404.14286 [gr-qc] .
  * Iglesias  _et al._ (2022) H. L. Iglesias _et al._ , (2022), arXiv:2208.01766 [gr-qc] .
  * Rodriguez  _et al._ (2018) C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, and F. A. Rasio, Phys. Rev. Lett. 120, 151101 (2018), arXiv:1712.04937 [astro-ph.HE] .
  * Samsing  _et al._ (2018a) J. Samsing, D. J. D’Orazio, A. Askar, and M. Giersz, (2018a), arXiv:1802.08654 [astro-ph.HE] .
  * Samsing  _et al._ (2018b) J. Samsing, A. Askar, and M. Giersz, Astrophys. J. 855, 124 (2018b), arXiv:1712.06186 [astro-ph.HE] .
  * Antonini and Gieles (2020) F. Antonini and M. Gieles, Mon. Not. Roy. Astron. Soc. 492, 2936 (2020), arXiv:1906.11855 [astro-ph.HE] .
  * Kremer  _et al._ (2020) K. Kremer, C. S. Ye, N. Z. Rui, N. C. Weatherford, S. Chatterjee, G. Fragione, C. L. Rodriguez, M. Spera, and F. A. Rasio, Astrophys. J. Suppl. 247, 48 (2020), arXiv:1911.00018 [astro-ph.HE] .
  * Parisi and Sturani (2018) A. Parisi and R. Sturani, Phys. Rev. D 97, 043015 (2018), arXiv:1705.04751 [gr-qc] .
  * Wang and Lai (2020) J.-S. Wang and D. Lai, Phys. Rev. D 102, 083005 (2020), arXiv:2009.08300 [astro-ph.HE] .
  * Bernaldez and Datta (2023) J. P. Bernaldez and S. Datta, Phys. Rev. D 108, 124014 (2023), arXiv:2303.01398 [gr-qc] .
  * Rathore  _et al._ (2005) Y. Rathore, R. D. Blandford, and A. E. Broderick, Mon. Not. Roy. Astron. Soc. 357, 834 (2005), arXiv:astro-ph/0406102 .
  * Fuller and Lai (2011) J. Fuller and D. Lai, Mon. Not. Roy. Astron. Soc. 412, 1331 (2011), arXiv:1009.3316 [astro-ph.HE] .
  * Burkart  _et al._ (2014) J. Burkart, E. Quataert, and P. Arras, Mon. Not. Roy. Astron. Soc. 443, 2957 (2014), arXiv:1312.4966 [astro-ph.SR] .
  * Chan  _et al._ (2014) T. K. Chan, Y. H. Sham, P. T. Leung, and L. M. Lin, Phys. Rev. D 90, 124023 (2014), arXiv:1408.3789 [gr-qc] .
  * Godzieba  _et al._ (2021) D. A. Godzieba, R. Gamba, D. Radice, and S. Bernuzzi, Phys. Rev. D 103, 063036 (2021), arXiv:2012.12151 [astro-ph.HE] .
  * Pradhan  _et al._ (2022) B. K. Pradhan, D. Chatterjee, M. Lanoye, and P. Jaikumar, Phys. Rev. C 106, 015805 (2022), arXiv:2203.03141 [astro-ph.HE] .
  * Pradhan  _et al._ (2023) B. K. Pradhan, D. Chatterjee, and D. E. Alvarez-Castillo, (2023), arXiv:2309.08775 [nucl-th] .
  * Alford and Harris (2019) M. G. Alford and S. P. Harris, Phys. Rev. C 100, 035803 (2019), arXiv:1907.03795 [nucl-th] .
  * Alford  _et al._ (2020) M. Alford, A. Harutyunyan, and A. Sedrakian, Particles 3, 500 (2020), arXiv:2006.07975 [nucl-th] .
  * Most  _et al._ (2022) E. R. Most, A. Haber, S. P. Harris, Z. Zhang, M. G. Alford, and J. Noronha, (2022), arXiv:2207.00442 [astro-ph.HE] .
  * Cutler  _et al._ (1990) C. Cutler, L. Lindblom, and R. J. Splinter, Astrophys. J. 363, 603 (1990).
  * Kochanek (1992) C. S. Kochanek, Astrophys. J. 398, 234 (1992).
  * Bildsten and Cutler (1992) L. Bildsten and C. Cutler, Astrophys. J. 400, 175 (1992).
  * Ho and Lai (1999) W. C. G. Ho and D. Lai, Mon. Not. Roy. Astron. Soc. 308, 153 (1999), arXiv:astro-ph/9812116 .
  * Foucart  _et al._ (2019) F. Foucart _et al._ , Phys. Rev. D 99, 044008 (2019), arXiv:1812.06988 [gr-qc] .
  * Ma  _et al._ (2020) S. Ma, H. Yu, and Y. Chen, Phys. Rev. D 101, 123020 (2020), arXiv:2003.02373 [gr-qc] .
  * Steinhoff  _et al._ (2021) J. Steinhoff, T. Hinderer, T. Dietrich, and F. Foucart, Phys. Rev. Res. 3, 033129 (2021), arXiv:2103.06100 [gr-qc] .
  * Lai and Wu (2006) D. Lai and Y. Wu, Phys. Rev. D 74, 024007 (2006), arXiv:astro-ph/0604163 .
  * Abbott  _et al._ (2020) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020), arXiv:2001.01761 [astro-ph.HE] .
  * Lyne  _et al._ (2004) A. G. Lyne _et al._ , Science 303, 1153 (2004), arXiv:astro-ph/0401086 .
  * Burkart  _et al._ (2012) J. Burkart, E. Quataert, P. Arras, and N. N. Weinberg, Mon. Not. Roy. Astron. Soc. 421, 983 (2012), arXiv:1108.3822 [astro-ph.SR] .
  * Burkart  _et al._ (2013) J. Burkart, E. Quataert, P. Arras, and N. N. Weinberg, Mon. Not. Roy. Astron. Soc. 433, 332 (2013), arXiv:1211.1393 [astro-ph.SR] .
  * Farr  _et al._ (2017) W. M. Farr, S. Stevenson, M. Coleman Miller, I. Mandel, B. Farr, and A. Vecchio, Nature 548, 426 (2017), arXiv:1706.01385 [astro-ph.HE] .
  * Cutler  _et al._ (1994) C. Cutler, D. Kennefick, and E. Poisson, Phys. Rev. D 50, 3816 (1994).
  * Lincoln and Will (1990) C. W. Lincoln and C. M. Will, Phys. Rev. D 42, 1123 (1990).
  * Blanchet (2014) L. Blanchet, Living Rev. Rel. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Vines and Flanagan (2013) J. E. Vines and E. E. Flanagan, Phys. Rev. D 88, 024046 (2013), arXiv:1009.4919 [gr-qc] .
  * Vines  _et al._ (2011) J. Vines, E. E. Flanagan, and T. Hinderer, Phys. Rev. D 83, 084051 (2011), arXiv:1101.1673 [gr-qc] .
  * Henry  _et al._ (2020) Q. Henry, G. Faye, and L. Blanchet, Phys. Rev. D 101, 064047 (2020), arXiv:1912.01920 [gr-qc] .
  * Yagi and Yunes (2013) K. Yagi and N. Yunes, Phys. Rev. D 88, 023009 (2013), arXiv:1303.1528 [gr-qc] .
  * Yagi and Yunes (2017) K. Yagi and N. Yunes, Phys. Rept. 681, 1 (2017), arXiv:1608.02582 [gr-qc] .
  * Maselli  _et al._ (2013) A. Maselli, V. Cardoso, V. Ferrari, L. Gualtieri, and P. Pani, Phys. Rev. D 88, 023007 (2013), arXiv:1304.2052 [gr-qc] .
  * Urbanec  _et al._ (2013) M. Urbanec, J. C. Miller, and Z. Stuchlik, Mon. Not. Roy. Astron. Soc. 433, 1903 (2013), arXiv:1301.5925 [astro-ph.SR] .
  * Press and Teukolsky (1977) W. H. Press and S. A. Teukolsky, Astrophys. J. 213, 183 (1977).
  * Samsing  _et al._ (2017) J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, Astrophys. J. 846, 36 (2017), arXiv:1609.09114 [astro-ph.HE] .
  * Chirenti  _et al._ (2017) C. Chirenti, R. Gold, and M. C. Miller, Astrophys. J. 837, 67 (2017), arXiv:1612.07097 [astro-ph.HE] .
  * Xu and Lai (2017) W. Xu and D. Lai, Phys. Rev. D 96, 083005 (2017), arXiv:1708.01839 [astro-ph.HE] .
  * Andersson and Ho (2018) N. Andersson and W. C. G. Ho, Phys. Rev. D 97, 023016 (2018), arXiv:1710.05950 [astro-ph.HE] .
  * Kuan  _et al._ (2021) H.-J. Kuan, A. G. Suvorov, and K. D. Kokkotas, Mon. Not. Roy. Astron. Soc. 506, 2985 (2021), arXiv:2106.16123 [gr-qc] .
  * Thrane  _et al._ (2014) E. Thrane, N. Christensen, R. M. S. Schofield, and A. Effler, Phys. Rev. D 90, 023013 (2014), arXiv:1406.2367 [astro-ph.IM] .
  * Dalal  _et al._ (2006) N. Dalal, D. E. Holz, S. A. Hughes, and B. Jain, Phys. Rev. D 74, 063006 (2006), arXiv:astro-ph/0601275 .
  * Regimbau  _et al._ (2012) T. Regimbau _et al._ , Phys. Rev. D 86, 122001 (2012), arXiv:1201.3563 [gr-qc] .