Inspiralling eccentric binary neutron stars: orbital motion and tidal resonance

TL;DR

This paper investigates the orbital dynamics and tidal resonances in eccentric binary neutron stars, revealing a significant phase shift in gravitational waves that could be detected by future detectors.

Contribution

It introduces a model for eccentric binary neutron star orbits including tidal effects and identifies a novel orbital resonance affecting gravitational wave signals.

Findings

Tidal interactions cause orbital precession beyond Post-Newtonian effects.

A new orbital resonance at one third of the f-mode frequency is identified.

Resonance can induce a detectable phase shift in gravitational waveforms.

Abstract

We study the orbital evolution of eccentric binary neutron stars. The orbit is described as a Quasi-Keplarian orbit with perturbations due to tidal couplings. We find that the tidal interaction between stars contributes to orbital precession, in addition to the Post-Newtonian procession. The coupling between the angular and radial motion of the binary also excites a series of harmonics in the stars' oscillation. In the small eccentricity limit, this coupling mainly gives rise to an additional orbital resonance, with the orbital frequency being one third of the f-mode frequency. For a binary with initial eccentricity $\sim 0.2$ at $50$Hz orbital frequency, the presence of this tidal resonance introduces $\sim 0.5$ phase shift in the gravitational waveform till merger, subject to uncertainties in neutron star equation of state and the distribution of binary component masses. Such phase shift in the late-inspiral stage is likely detectable with third-generation gravitational-wave detectors.

Figures

Equations

$$\mathcal{H}$$

$$={\mathcal{H}}_{\mathrm{orb}}+{\mathcal{H}}_{\mathrm{mode}}+{\mathcal{H}}_{\mathrm{int}},$$

$$Q_{ij}=\sum _n{\mathcal{O}}_{ij}^n{Q}_n,$$

$${\mathcal{O}}_{ij}^n=\int \rho d^{3}x({\xi }_i^n{x}_j+{\xi }_j^n{x}_i)$$

$${\xi }^{(n)}=\left({\xi }_R^{(n)}(r){\mathbf{e}}_{\mathbf{r}}+{\xi }_S^{(n)}(r)r\nabla \right)Y_{lm}(\theta ,\varphi ).$$

$$\xi =\sum _n{Q}_n{\xi }^n,$$

$${\mathcal{O}}_{ab}^n=Q_{\xi }^n\int d\Omega ({\mathbf{e}}_a\cdot {\mathbf{e}}_r)({\mathbf{e}}_b\cdot {\mathbf{e}}_r).$$

$$Q_{\xi }^n=2{\int }_0^{R_{\ast }}dr{r}^{\ell +1}\rho [{\xi }_R^{(n)}+(\ell +1){\xi }_S^{(n)}].$$

$$\mathbf{x}$$

$$\mathbf{p}$$

$$E$$

$$J$$

$$P_{\mathrm{orb}}^E+P_{\mathrm{mode}}^E+P_{\times }^E$$

$$=\dot{E}+{\left[{ϵ}^2\frac{\partial {\mathcal{H}}_{\mathrm{orb}}}{\partial \mathbf{x}}\delta x\right]}_{,t}+{ϵ}^2{\left[\frac{\partial {\mathcal{H}}_{\mathrm{orb}}}{\partial \mathbf{p}}\delta p\right]}_{,t}$$

$$+{ϵ}^2{\dot{E}}_{\mathrm{mode}}+{ϵ}^2{\dot{E}}_{\mathrm{int}}$$

$$=\dot{E}+\frac{\partial {\mathcal{H}}_{{ϵ}^2}}{\partial E}\dot{E}+\frac{\partial {\mathcal{H}}_{{ϵ}^2}}{\partial J}\dot{J},$$

$$P_{\mathrm{orb}}^J+P_{\mathrm{mode}}^J+P_{\times }^J$$

$$=\dot{J}+{\left[{ϵ}^2\frac{\partial {\mathcal{J}}_{\mathrm{orb}}}{\partial \mathbf{x}}\delta x\right]}_{,t}+{ϵ}^2{\left[\frac{\partial {\mathcal{J}}_{\mathrm{orb}}}{\partial \mathbf{p}}\delta p\right]}_{,t}$$

$$+{ϵ}^2{\dot{J}}_{\mathrm{mode}}$$

$$=\dot{J}+\frac{\partial {\mathcal{J}}_{{ϵ}^2}}{\partial E}\dot{E}+\frac{\partial {\mathcal{J}}_{{ϵ}^2}}{\partial J}\dot{J},$$

$${\ddot{Q}}^n+\gamma {\dot{Q}}^n+{\omega }_n^2{Q}^n=-ϵ{\mathcal{E}}_Q({\mathbf{x}}_{\mathrm{orb}},{\mathbf{p}}_{\mathrm{orb}}),$$

$$\dot{\delta \mathbf{p}}=-\frac{{\partial }^2{\mathcal{H}}_{\mathrm{orb}}}{\partial \mathbf{x}\partial \mathbf{x}}\delta \mathbf{x}-\frac{{\partial }^2{\mathcal{H}}_{\mathrm{orb}}}{\partial \mathbf{x}\partial \mathbf{p}}\delta \mathbf{p}-\frac{\partial {\mathcal{H}}_{\mathrm{int}}}{\partial \mathbf{x}},$$

$$\dot{\delta \mathbf{x}}=\frac{{\partial }^2{\mathcal{H}}_{\mathrm{orb}}}{\partial \mathbf{p}\partial \mathbf{x}}\delta \mathbf{x}+\frac{{\partial }^2{\mathcal{H}}_{\mathrm{orb}}}{\partial \mathbf{p}\partial \mathbf{p}}\delta \mathbf{p}.$$

$$Q^n=$$

$$+Q_{\mathrm{init}}^n{e}^{-i{\omega }_{n}t-\gamma t}=Q_{\mathrm{driven}}^n+Q_{\mathrm{free}}^n,$$

$$Q_{\mathrm{driven}}^n\approx Q_{\mathrm{ad}}^n=-ϵ\frac{{\mathcal{E}}_Q(t)}{{\omega }_n^2-2i\gamma \omega -{\omega }^2}.$$

$$P^E$$

$$P_i^J$$

$$P_{\mathrm{mode}}^E$$

$$=-\frac{8\pi }{75}{Q}_{\xi }^2⟨\sum _m|{\stackrel{...}{Q}}_m{|}^2⟩$$

$$=-\frac{4\pi M_{\ast }^2{Q}_{\xi }^4{n}^6}{75q^6}\sum _m{W}_{2m}^2\sum _k\frac{k^6(|c_{3,k}^m{|}^2+|s_{3,k}^m{|}^2)}{[(kn{)}^2-{\omega }_f^2{]}^2},$$

$$P_{\mathrm{mode}}^J$$

$$=-\frac{i16\pi }{75}{Q}_{\xi }^2⟨{\ddot{Q}}_2{\stackrel{...}{Q}}_{-2}-{\ddot{Q}}_{-2}{\stackrel{...}{Q}}_2⟩$$

$$=-\frac{32\pi }{75}\frac{W_{22}^2{M}_{\ast }^2}{q^6}{Q}_{\xi }^4{n}^5\sum _k\frac{k^5{s}_{3,k}^2{c}_{3,k}^{-2}}{[(kn{)}^2-{\omega }_f^2{]}^2},$$

$$P_{\times }^E=$$

$$=\frac{4\pi }{25}\mu \frac{M^{3/2}}{q^{5/2}}\frac{M_{\ast }{Q}_{\xi }^2}{q^3}{n}^3\sum _k{k}^3\frac{(4s_{3,k}^2-0.5e{s}_{2,k}^3+0.5e{s}_{2,k}^1)c_{3,k}^2+(4c_{3,k}^2-0.5e{c}_{2,k}^3+0.5e{c}_{2,k}^1)s_{3,k}^2}{(kn{)}^2-{\omega }_f^2}$$

$$-\frac{8\pi }{75}\mu \frac{M^{3/2}}{q^{5/2}}e\frac{M_{\ast }{Q}_{\xi }^2}{q^3}{n}^3\sum _k{k}^3\frac{s_{2,k}^1{c}_{3,k}^2}{(kn{)}^2-{\omega }_f^2},$$

$$P_{\times }^J=$$

$$=$$

$$+\frac{4\pi }{25}\frac{\mu M^{3/2}}{q^{5/2}}\frac{M_{\ast }{Q}_{\xi }^2}{q^3}{n}^2\sum _k{k}^2\frac{(4c_{3,k}^2-0.5e{c}_{2,k}^3+0.5e{c}_{2,k}^1)c_{3,k}^2+(4s_{3,k}^2-0.5e{s}_{2,k}^3+0.5e{s}_{2,k}^1)s_{3,k}^2}{(kn{)}^2-{\omega }_f^2}.$$

$$Q_{ij}={\mathcal{T}}_{ij\alpha \beta }{\mathcal{E}}^{\alpha \beta }.$$

Full text

I Introduction

Neutron star mergers produce copious gravitational wave (GW) and electromagnetic (EM) radiation, that encode critical information about state of matter under the extreme physical conditions (gravity, pressure, temperature, etc.). In particular, the detection of GW170817 Abbott et al. (2017a) and subsequent observations from gamma-ray band to radio-wave frequencies has ushered in a new era of multi- messenger astronomy Abbott et al. (2017b, c, d); Coulter et al. (2017); Troja et al. (2017); Alexander et al. (2017); Kochanek and Piran (1993); Li and Paczyński (1998). As the sensitivity of current-generation GW detectors degrades significantly above $1$kHz Miao et al. (2018); Martynov et al. (2019), directly probing the post-merger GW signal from binary neutron stars (BNSs) is less likely in a few years. The inspiral part of the GW signal, however, not only provides information about the binary orbital parameters, such as component masses, eccentricity, distance, etc., but also the matter properties through the measurement of tidal deformability and possible resolution of the merger frequency. For example, the tidal love number constraint from GW170817 has greatly refined the possible parameter range of NS EOS.

Most of the compact binary mergers are believed to be (nearly) circular 111A recent publication has argued that Post-Newtonian radiation reaction excites small eccentricities $\mathcal{O}(0.01)$ in the late inspiral stage Loutrel et al. (2018), which appears to be interesting to be the discussion in this work.. However, recent studies have shown the possibility of forming eccentric compact binaries in the LIGO (Laser Interferometric Gravitational Wave Observatory) band by resonant and hierarchical triple and quadruple systems in globular clusters Wen (2003); Seto (2013); Naoz (2016); Antonini et al. (2017); Rodriguez and Antonini (2018); Samsing and Ramirez-Ruiz (2017); Samsing (2018); Liu and Lai (2018); Hoang et al. (2018), in addition to dynamically captured binaries East et al. (2013). These events may be rare comparing to circular inspirals, but they contain important information about their formation channel(s), environment and distribution, through the spin and eccentricity measurement. For eccentric binary neutron stars, we shall show that the coupling between radial and angular orbital motions give rises to new tidal resonances, that can be possibly detected by third-generation GW detectors such as Einstein Telescope Punturo et al. (2010) and Cosmic Explorer Abbott et al. (2017e).

The motion of inspiring eccentric binary black holes (BBHs), within the Post-Newtonian framework, can be described by the Quasi-Keplerian (QK) orbits Blanchet (2014). At the zeroth order, the QK orbit coincides with the Newtonian elliptical orbit. With higher order Post-Newtonian effect included, the radial motion generically oscillates with a different period from the angular motion, which is known as the general relativistic precession. In addition, radial and angular motions also receive periodic corrections that can be expanded in a Post-Newtonian manner. Alternatively, it is possible to utilize the Effective-One-Body (EOB) framework, which resembles the Post-Newtonian expansion in the EOB spacetime Hinderer and Babak (2017). Based on these frameworks, different waveform models have been developed to characterize the GW signature of BBHs, which have achieved in decent accuracy for low-medium eccentricity binaries Huerta et al. (2018, 2017a); Cao and Han (2017); Hinder et al. (2018), with promising potential to realize fast frequency-domain waveform model with arbitrary eccentricities Huerta et al. (2017b); Moore et al. (2018).

For eccentric BNSs, the tidal coupling between stars introduces extra complication in the equation of motion. For example, in Yang et al. (2018a) we have studied highly eccentric BNS in the inspiral stage, where the f-modes of NSs are excited during the close pericenter passages, and subsequently play important roles in the orbital evolution. With sufficiently accurate orbit model, it might be possible to combine f-mode emission from different pericenter encounter cycles to boost the signal-to-noise ratio of detection, e.g., with the coherent stacking method Yang et al. (2017, 2018b); Berti et al. (2018). While the formation of such systems may require rather restrictive initial parameters, it is still important to characterize the tidal effect in medium-low eccentricity BNSs, in addition to the understanding of equilibrium and dynamic tide in the circular orbit limit Steinhoff et al. (2016); Flanagan and Hinderer (2008).

In this work, we adopt the QK orbit as the unperturbed solution without the tidal effect, and compute the NS oscillation in response to the orbital motion and tidal couplings. For eccentric orbits, the tidal budge induced on a star generally does not point to the companion star in the binary. As a result, the binary orbital angular momentum continuously exchanges with the NS mode angular momentum within orbital timescales. These NS oscillations also back-react on the binary orbit, giving rise to perturbations of the QK orbit that affect the GW radiation. In the low eccentricity limit, we explicitly evaluate these perturbation terms and determine their influence on the orbital phases in the time domain. Although the discussion of the NS mode and the orbital dynamics has been specified to BNSs, the result is also applicable to BH-NS binary Yang et al. (2018c), with the oscillation of the BH set to be zero. The construction of a frequency-domain waveform model will be left to further work.

This paper is organized as follows. In Sec. II we present the basic formalism to describe the motion of an eccentric BNS system under the influence of tidal couplings. In Sec. III we apply this formalism to the Newtonian orbits, characterize the importance of tidal effects for low-eccentricity systems and discuss the tidal resonance effect. In Sec. IV we repeat the exercise for Post-Newtonian orbits. We conclude in Sec. V. Natural units with $c=1,G=1$ are used through out the analysis.

II Basic Formalism

The Hamiltonian of a BNS system, including the leading-order tidal excitations of the stars, can be expressed as Steinhoff et al. (2016); Flanagan and Hinderer (2008)

[TABLE]

where the tidal field is labeled with $\epsilon$ for book keeping purpose, and $\omega_{n}$ is the eigenfrequency for mode $n$. In stellar perturbation theory Unno et al. (1979) the modes of a three-dimensional star are often indexed by the angular nodal number $\ell$, azimuthal nodal number $m$ and radial nodal number $n_{r}$. Here $n$ represents the collection of these three indices: $n=\{\ell mn_{r}\}$. In particular, the gravitational response of f-mode generally dominates over other modes (e.g., p-modes and g-modes) of the NS, so that in this work we primarily focus on the f-mode excitations. The electric part of the tidal tensor, $\mathcal{E}_{ij}$, can be evaluated based on the relative displacement of the binary Poisson and Vlasov (2010); Yang et al. (2018a), and the local spacetime of the target star is influenced by the tidal environment generated by its companion (also see the application in extreme mass-ratio inspiral systems in Yang and Casals (2017); Bonga et al. (2019)). The orbital Hamiltonian $\mathcal{H}_{\rm orb}$, expanded in the Post-Newtonian format, can be found in Memmesheimer et al. (2004) up to the $3$rd Post-Newtonian order. Because of the centre-of-mass conservation, the orbital motion is fully characterized by the relative distance and orbital angle, which orginally motivated the development of the EOB formalism Damour (2014). The modal quadrupole moment $Q^{n}_{ij}$ induced by the tidal field, and the modal excitation displacement $Q^{n}$, are related to each other through

[TABLE]

with the overlap tensor defined as

[TABLE]

where $\rho$ is the local density within the star and ${\bf\xi}^{n}$ is the displacement vector associated with the eigenmode wavefunction, that can be expressed in vector spherical harmonics:

[TABLE]

The determination of $\xi^{(n)}_{R}$ and $\xi^{(n)}_{S}$ is discussed in Appendix A of Yang et al. (2018a). Notice that Eq. (II.2) can be derived by combining Eq. (II.3) the expansion

[TABLE]

and it can be further represented as

[TABLE]

with

[TABLE]

At this point, all variables are still real. But later on when we talk about specific modes, it is often convenient to use complex wave functions for decomposition. In those cases, $Q_{n}$ can be complex-valued, although the total displacement $\xi$ remains real.

II.1 Orbital description

Under the influence of tidal field, the binary orbit is no longer eccentric, even in the Newtonian limit. However, as the tidal interaction is generally weaker than the direct point-mass gravitational attraction, we can expand the orbit as power laws in the tidal-coupling coefficient $\epsilon$

[TABLE]

where ${\bf x},{\bf p}$ follow an eccentric orbit in the Newtonian description or QK orbit in the Post-Newtonian description, and the orbital evolution is determined by the conserved energy and angular momentum $E,J$.

[TABLE]

where $\mathcal{J}_{\rm orb}$ can also be found in Memmesheimer et al. (2004) up to the 3rd Post-Newtonian order. The perturbation starts at the second order in $\epsilon$ as $Q$ starts at linear order in $\epsilon$, so that the back-reaction of $Q$ starts at the second order. As the binary orbit decays due to GW radiation generated by orbital motion and mode oscillation, we have

[TABLE]

where different components of the energy and angular momentum flux $P^{E,J}$ are discussed in Sec. II.2. To evolve the orbit, we need to relate them to the secular change of $E$ and $J$. This means that we need to solve the equation of motion for $\delta{\bf x},\delta{\bf p},Q$ as a function of $E,J$ and time $t$. Based on the Hamiltonian, the equations of motion are

[TABLE]

where $\gamma$ is an infinitesimal positive damping rate due to dissipations in the star that have not been accounted for in the Hamiltonian formalism, and $\mathcal{E}_{Q}:=\mathcal{E}_{ij}\mathcal{O}^{*}_{ij}/2$. Notice that we are taking complex conjugation for $\mathcal{O}_{ij}$ here because the complex wave function is being used. The solution of $Q$ is

[TABLE]

where $Q^{n}_{\rm init}$ corresponds to the initial value of $Q^{n}$.

Like a normal weakly-damped harmonic oscillator, $Q^{n}$ generically contains a “driven” part that is proportional to $\epsilon$ and a “free” that satisfies the homogeneous equation of motion (the left hand side of the first line of Eq. (II.1)). Under a periodic driving force with slowly varying frequency $\omega/(2\pi)$ ($\dot{\omega}/\omega\ll\omega$), the driven part $Q_{\rm driven}$ can often be approximated by its adiabatic value

[TABLE]

On the other hand, if the initial oscillation is zero or the initial time of integration is far in the past, such “free” part may be neglected. However, as we shall see later, when the system evolves across mode resonances, certain level of free mode oscillation will be excited as well.

II.2 The radiative terms

The radiation terms that appear in Eq. (II.1) can be evaluated from the quadruple formula:

[TABLE]

Because we have assumed both NSs to be non-spinning, their orbital motion should still remain on a plane with the presence of Post-Newtonian correction and tidal effects. As a result, the orbital angular momentum is orthogonal to the orbital plane, so that only the normal component (defined as the $z$ direction) of $P^{J}$ is non-vanishing, and hereafter we drop the vector indices of $P^{J}$.

The orbital energy and angular momentum flux $P^{E}_{\rm orb},P^{J}_{\rm orb}$ are derived in Arun et al. (2008) up to the third Post-Newtonian order. The modal fluxes $P^{E}_{\rm mode},P^{J}_{\rm mode}$ are given by

[TABLE]

and Yang et al. (2018a)

[TABLE]

where we have used Eq. (III.1) and the fact that only the $\ell=2,m=0,\pm 2$ f modes are relevant for the discussions here (The $m=\pm 1$ modes are not excited Press and Teukolsky (1977); Yang et al. (2018a)). These expressions can be made more explicit with the prescription of an unperturbed orbit.

The cross terms $P^{E}_{\rm\times}$ and $P^{J}_{\rm\times}$ comes from the beating between radiation from the orbital motion and the mode excitation.

[TABLE]

and

[TABLE]

In the Newtonian limit, $\ddot{I}_{ij}$ and $\dddot{I}_{ij}$ can be found in Eq.(12.77) of Poisson and Will (2014) for generic eccentric orbits.

II.3 Equilibrium and dynamic tide

The tidal response in the static limit can be formally described by a transformation tensor

[TABLE]

For spherically symmetric object, this transformation tensor reduces to a single tidal Love numbers $\lambda$, such that

[TABLE]

When the frequency for the tidal field $\omega$ is much less then the mode frequencies, we can approximate the solution of Eq. (II.13) as

[TABLE]

which is often referred as the the equilibrium tide approximation. The total induced quadrupole moment is

[TABLE]

Such expression should be compared with Eq. (II.20). While it is not immediately clear why these two expressions are equivalent, they are guaranteed by the spherical symmetry of the star, and $\lambda$ can be obtained as

[TABLE]

where in the second line we have kept the contribution from f-modes. On the other hand, for the case with

[TABLE]

we have

[TABLE]

The tidal response is frequency dependent, which is often referred as the dynamic tide. Certain higher harmonic frequency $\Omega_{\alpha}$ may be comparable or even larger than the f-mode frequency $\omega_{f}$. Resonance occurs when $\Omega_{\alpha}\approx\omega_{f}$.

III Newtonian orbits

In this section, we keep only the leading order Newtonian term in the Hamiltonian, and solve for the orbital evolution in response to the tidal coupling. We also examine the tidal resonances in the low eccentricity limit, and discuss the criteria of detecting such resonances with current and future GW detectors. For simplicity, we only include the mode evolution for one reference star, as it is straightforward to extend the analysis to oscillations of both stars in the BNS system.

The Newtonian Hamiltonian for the orbit is

[TABLE]

with $\mu=M_{1}M_{*}/M$ and $M=M_{1}+M_{*}$. The unperturbed orbit can be characterized as

[TABLE]

with $n=\sqrt{M/a^{3}}$, $u$ being the mean anomaly, $l$ being the mean motion and $\phi$ being the true anomaly, and $t_{p},\phi_{p}$ corresponding to the time and angle at the pericenter passage. We further have energy and eccentricity given by

[TABLE]

Now according to Eq. (II.1), we expand the motion as

[TABLE]

The equations of motion (c.f. Eq. (II.1)) become

[TABLE]

and the solution of the mode excitation is given by Eq. (II.1). For simplicity we neglect the initial oscillation by setting $Q_{\rm free}=0$.

III.1 Dynamic tide

The dynamical-tide excitation corresponds to the tidal response described by Eq. (II.1). For the “background” trajectory described by Eq. (D), the tidal tensor generated by the companion NS with mass $M_{*}$, that acts on the reference NS, is Yang et al. (2018a) (taking $\phi_{p}=0$ here)

[TABLE]

so that (restricting to $\ell=2$ subspace)

[TABLE]

Here $W_{lm}$ are defined in Press and Teukolsky (1977), with relevant components used in this work:

[TABLE]

In order to solve the equations of motion in Eq. (III), let us define Rathore et al. (2005)

[TABLE]

Here $c^{m}_{n,k},s^{m}_{n,k}$ are functions of the eccentricity $e$, which are proportional to the Hansen coefficients Murray and Dermott (1999) (apart from a $(1-e^{2})^{n}$ factor). They can be obtained through

[TABLE]

and

[TABLE]

In terms of the Hansen coefficients, and according to Eq. (II.1), the f-mode excitation becomes ($\gamma\rightarrow 0$)

[TABLE]

where in the last line we have adopted the adiabatic approximation (c.f. Eq. (II.13)), and $q:=a(1-e^{2})$. The star is driven at integer harmonics of the orbital frequency, so that it is possible to have resonance crossing during the inspiral stage. For circular orbits, only the term with $k=2$ survives. In addition, the inspiral usually terminates at the ISCO (Innermost-Stable-Circular-Orbit) frequency or the contact frequency of the two NSs, which are comparable or smaller than the f-mode frequency divided by two. Therefore it is difficult to observe a complete resonance during the inspiral stage for $k=2$. However, it has been shown Steinhoff et al. (2016) that the frequency dependence of the tidal love number is important for describing motion in the late inspiral stage, which is essentially related to the $k=2$ resonance.

We can similarly decompose the driving terms in the equations of motion into a summation of harmonics:

[TABLE]

and

[TABLE]

and

[TABLE]

As a result, we can determine Eq. (III) in an expansion of harmonics. We first write $\delta r$ and $\delta\phi$ as

[TABLE]

There are no $\sin kl$ terms because of the absence of corresponding terms in the driving force. The equations of motion become ($k\geq 0$)

[TABLE]

These equations can be solved in the matrix form. The $k=0$ piece of angular momentum shift $g_{0}$ is zero, because there is no DC angular momentum exchange between the orbit and the stars.

III.2 Small eccentricity limit

In the small eccentricity limit $e\ll 1$, we can take the leading order expansion of the Hansen coefficients in terms of $e$. They can be found in the Appendix A. In particular, for the terms showing up in Eq. (III.1), we always have $m=\pm 2,0$ and $n=3,4$. The nonzero components that are proportional to $e$ are

[TABLE]

The principle part that survives in the circular limit is

[TABLE]

With these results we now try to solve Eq. (III.1) in the small eccentricity limit. The result is

[TABLE]

and we find that for $k=1$, the coefficient of $b_{1}$ becomes zero, which means that the $k=1$ term corresponds to the resonant frequency of the orbit. Any external driving frequency being the same as this frequency will formally make $b_{1}$ diverge. Physically what happens is that the coupling with star’s internal degrees of freedom shifts the radial frequency, and make the orbit precess. In order to fix this problem, we need to assign a different frequency $n+\delta n$ to the radial motion with $\delta n\neq 0,\delta n/n\sim\mathcal{O}(\epsilon^{2})$. The $k=1$ motion should also be absorbed into the background trajectory, with a redefinition of the eccentricity. By expanding the radial equations of motion in $\epsilon$ and keep terms linear in $\epsilon^{2}$, we see that

[TABLE]

with $L_{\rm orb}\approx\mu na^{2}$. As a result, we have

[TABLE]

All other terms ($b_{k},g_{k}$) in the expansion series are zero. With $\delta r$ and $\delta L_{\phi}$ known, it is then straightforward to compute the tidal perturbation of the Hamiltonian. Both $\delta x$ and $\delta L_{\phi}$ are time-dependent, but the Hamiltonian perturbation is not. Therefore we can pick any time to evaluate $\delta\mathcal{H}$. It turns out that $l=\pi/2$ is a convenient choice.

[TABLE]

Similarly for the angular momentum:

[TABLE]

We notice that

[TABLE]

so that $\delta J$ is conserved. At $l=\pi/2$, $\delta L_{\phi}$ is zero (it has no $k=0$ component) and $J_{\rm mode}$ is given by Steinhoff et al. (2016)

[TABLE]

where we have used the identification $\lambda=4\pi Q^{2}_{\xi}/(15\omega^{2}_{f})$ for f mode. In addition, because

[TABLE]

we can evaluate the perturbation of the angular frequency $n+\delta n_{\phi}$ due to the tidal interaction, such that

[TABLE]

Because of the modification of the radial and angular frequencies $\delta n,\delta n_{\phi}$, and the shift of the trajectory ($\delta r,\delta L_{\phi}$), the principle parts of the energy and angular momentum radiation are correspondingly changed.

[TABLE]

The explicit evaluation of these quantities are discussed in detail in Appendix D.

III.2.1 Sample evolution

In order to illustrate the effect of eccentric tidal terms in the flux, energy and angular momentum derived in Sec. III and Appendix. C, we take an equal-mass, $1.3M_{\odot}+1.3M_{\odot}$, binary NS system as an example. The star compactness is assumed to be $M_{*}/R_{*}=0.16$, and the initial eccentricity at $f_{\rm orb}=50{\rm Hz}$ is set to be $0.2$. We use this system to compare two separate evolution schemes. In the first evolution, we adopt Eq. C, Eq. C, Eq. D, Eq. D, Eq. D and Eq. III.2 for a time-domain evolution to obtain a relation between $f_{\rm orb}$ and $t_{\rm full}$. In the second evolution, we drop the terms proportional to $e^{2}\lambda$ in these formulas, and evolve the system again to obtain the function $t_{\rm part}(f_{\rm orb})$. The time difference up to certain orbital frequency, $\delta t$, is defined to be $t_{\rm full}-t_{\rm part}$.

The phase difference accumulated for the dominate $\ell=2,m=2$ mode, can be estimated as

[TABLE]

for this particular evolution if we choose $f_{\rm down}=50{\rm Hz},f_{\rm up}=500{\rm Hz}$. Such phase difference is probably only observable for loud events in the third-generation detector era.

III.3 Tidal resonance

According to Eq. (III.1), when $kn\approx\omega_{f}$ we expect the f mode to be resonantly excited. As the orbit does not stay at resonance for infinite time because of radiative dissipation, we need to account for the total phase shift across the resonance. In this case, we need to consider Eq. (II.1) taking into account the evolution of conserved quantities.

There are two major effects of these tidal resonances, which are intimately related to each other. At first, f mode is resonantly excited across the resonance, and this free oscillation of f modes persists in the post-resonance stage. This f-mode free oscillation also generates GW radiation at frequency $\omega_{f}/(2\pi)$. Secondly, as the orbital energy and angular momentum transfers to the f-mode to support the free oscillation, the orbital motion after the resonance gradually deviates away from the one without the resonance, due to the change of the conserved quantities. In this case, the resulting orbital phase shift can also be determined.

III.3.1 Free oscillation

The orbit evolution across resonance requires that there is a free oscillation piece of $Q_{m}$, excited by the driving force. In other words, the last line of Eq. (III.1) needs to be modified in the near-resonance regime, which introduces a free oscillation piece.

If $kn=\omega_{f}+k\dot{n}(t-t_{f})+\mathcal{O}(t-t_{f})^{2}$, then the free oscillation piece of $Q_{m}$ is

[TABLE]

with the Fresnel functions

[TABLE]

We note that $F_{\rm c}(\infty)=F_{\rm s}(\infty)=\sqrt{\pi/2}$. The idea is that $Q^{\rm ad}_{m}$ smoothly transit to $Q^{\rm free}_{m}$ when $kn\rightarrow\omega_{f}$:

[TABLE]

It can be shown that the transition is smooth (results do not depend sensitively on $t_{1},t_{2}$), for example, using the asymptotic behaviour discussed in Steinhoff $\&$ Hinderer. However, they miss this free oscillation piece in the post-resonance regime ($t>t_{2}$).

Free oscillation of f-modes has been observed in highly eccentric binaries Gold et al. (2012); Yang et al. (2018a) 222For the implusive approximation in Yang et al. (2018a) to hold, the eccentricity needs to be generally larger than $\sim 0.6$. As the eccentricity generally decays with time, the impulsive approximation may break down in the later part of the evolution, when $e<0.6$. . While the orbital timescale may be long, the collective work of many harmonics gives rise to the impulsive interaction near the pericenter. In general, we need to consider $(\omega_{f}/n)$th order harmonic to reproduce this “free oscillation” feather. If the orbit is highly eccentric, this task is computationally expensive, and it is preferable to apply the impulse approximation as in Yang et al. (2018a).

III.3.2 Orbital phase shift

According to Eq. (III.3.1), the amplitude of the f-mode free oscillation gained after the $kn\approx\omega_{f}$ resonance is

[TABLE]

The orbital energy decreases an extra amount across the resonance, due to the excitation of the free oscillation. As a result, the binary should merge faster than those without resonance crossing. In the small eccentricity limit, the only relevant tidal resonances are the associated with $k=2$ or $k=3$. The $k=2$ resonance, although being largest in amplitude, hardly takes place in the inspiral part because the binary merger frequency is generally smaller than half of the f-mode frequency. On the other hand, the mode amplitude for the $k=3$ resonance is proportional to the eccentricity. The orbital frequency associated with this resonance, being one third of the f-mode frequency, could be smaller than the merger frequency.

Therefore, for $k=3$ resonance, the corresponding $(2,2)$ component phase shift of the gravitational waveform is Lai (1994) (for equal mass binaries):

[TABLE]

where the first factor of 2 comes from the mode energy for both NSs. Here $t_{\rm orb}$ is the orbital period at the resonance, $e$ is the eccentricity at the resonance, $t_{D}=f/\dot{f}$ at the resonance and $E_{\rm orb}=M\mu/(2a)$ is the orbital energy at the resonance. The quantity $\mathcal{Q}$ is the dimensionless tidal overlap coefficient defined as

[TABLE]

Notice that the orbital eccentricity $e$ decays with increasing orbital frequency $f$ due to GW radiation, with the scaling being approximately $e\approx f^{-19/18}$. This means that a binary with eccentricity $\sim 0.02$ at $500$ Hz, would have eccentricity $\sim 0.2$ at $50$ Hz and $\sim 0.6$ at $20$ Hz. Such binaries may originate from dynamical captures in globular clusters O’Leary et al. (2009) and multi-body dynamic evolutions. Despite of the recent developments, the rate of these channels are still subject to significant theoretical uncertainties.

The statistical phase error of an event with signal-to-noise ratio (SNR) $\rho$ is approximately $\sqrt{D-1}/\rho$ Chatziioannou et al. (2017), with $D$ being the number of degrees of freedom in the parameter estimation. To resolve a phase $\mathcal{O}(0.5)$ error with $D\sim\mathcal{O}(10)$, events with $\rho\geq 7$ is necessary. Note that this resonance happens in the late-inspiral stage, which is around $1$kHz for the quadrupole GW radiation, so that $\rho$ should represent the SNR of the waveform segment starting from the resonance and ending at merger 333If the post-merger waveform can be accurately modelled, $\rho$ should also include the merger and post-merger SNR.. For a $1.3M_{\odot}+1.3M_{\odot}$ binary NS system at $40{\rm Mpc}$, the SNR for the inspiral waveform from $1$kHz to the plunge frequency ($\sim 1.4{\rm kHz}$) is around 1.5 for Advanced LIGO design sensitivity Martynov et al. (2016) and around $20$ for Cosmic Explorer Abbott et al. (2017e). Therefore, such requirements are more likely satisfied with the third generation GW detectors, such as the Einstein Telescope, Cosmic Explorer, LIGO-HF which also has decent mid-frequency sensitivity Miao et al. (2018); Martynov et al. (2019), etc.

If one (or both) NS(s) is a milli-second pulsar, the rotation frequency $f_{s}$ may be a couple of hundred Hz. As a result, the degeneracy between different f modes with different azimuthal number $m$ is broken. In particular, the frequency split is roughly $0.5mf_{s}$ Lee and Strohmayer (1995). If the NS spin counter-rotates with the binary, the mode with frequency $f_{\rm mode}-f_{s}$ is mostly excited; if the NS spin co-rotates with the binary, the mode with frequency $f_{\rm mode}+f_{s}$ is mostly excited. Therefore for counter-rotating binaries, if the $f_{\rm mode}-f_{s}$ is smaller than the merger frequency, the $k=2$ tidal resonance is present in the inspiral stage Steinhoff et al. (2016)444Private communication with Yanbei Chen.. On the other hand, for the $k=3$ tidal resonance, the mode frequency in Eq. (III.3.2) is $f_{\rm mode}-f_{s}$ and the corresponding orbital frequency will be modified to $(f_{\rm mode}-f_{s})/3$.

IV Tidally Modified QK orbits

In the Post-Newtonian limit, the motion of two gravitationally bounded points masses can be described by the Quasi-Keplerian orbit, as a generalization of the Keplerian orbit in the Newtonian theory:

[TABLE]

with

[TABLE]

Here the Post-Newtonian correction $l_{\rm PN},\phi_{\rm PN}$ are functions of $u$ or $v$. Their detailed expressions, together with the expressions for $e_{r},e_{\phi},e_{t},n^{\prime},K^{\prime}$ can be found in Blanchet (2014). Unlike a Keplerian orbit, the radial and azimuthal frequencies of a QK trajectory are no longer degenerate.

Notice that if we are only interested in obtaining leading-order tidal effect in the Post-Newtonian expansion, we only need to plug in the Newtonian part of $\mathcal{H}_{\rm orb}$ in Eq. (II.1), but sticking to Eq. (IV) for the prescription of the “background” trajectory. In fact, we can further neglect $\phi_{\rm PN}$ and $l_{\rm PN}$ for this purpose because they contribute to oscillatory Post-Newtonian effect, as comparing to the secular Post-Newtonian effect encoded in $e_{r},e_{t},e_{\phi},K^{\prime},n^{\prime}$. Let us write $\delta r,\delta L_{\phi}$ as

[TABLE]

$K_{\phi}=K^{\prime}+\delta K^{\prime}$ is different from $K^{\prime}$ because of the tidal correction. Similarly $n^{\prime}$ needs to be promoted to $n_{t}=n^{\prime}+\delta n^{\prime}$ by the tidal correction as well. In the Newtonian, low-eccentricity limit, such mappings are described by Eq. III.2 and Eq. III.2. Here $n^{\prime}$ can also be determined by the equation of motion for the $k^{\prime}=0,k=1$ component of orbital perturbations.

IV.1 Equations of motion

The mode excitation of the star can still be determined from (in the non-resonant regime)

[TABLE]

where the QK-Hansen coefficients $c^{m}_{n,k,k^{\prime}},s^{m}_{n,k,k^{\prime}}$, which depend on $e_{t,r,\phi},K^{\prime}$, can be defined by

[TABLE]

and

[TABLE]

It it straightforward to check that the QK-Hansen coefficients $c^{m}_{n,k,k^{\prime}},s^{m}_{n,k,k^{\prime}}$ reduce to $c^{m}_{n,k}(1-e^{2})^{-n},s^{m}_{n,k}(1-e^{2})^{-n}$ if we take the limit $e_{r,t,\phi}=e,K^{\prime}=1$.

Strictly speaking, Eq. III is no longer valid as $\mathcal{H}_{\rm orb}$ in Eq. II and Eq. II.1 now contain higher-order Post-Newtonian corrections. However, the leading-Post-Newtonian tidal perturbation in $\delta r,\delta L_{\phi}$ can be still obtained using Eq. III with $Q_{m}$ given in Eq. IV.1. In particular, the equations of motion for $\delta L_{\phi}$ imply that ($k_{0}\equiv k_{1}K^{\prime}+k_{2}$)

[TABLE]

This implies that

[TABLE]

Similarly for $\delta r$ we have

[TABLE]

The $k_{0}$ is selected within a 2-D lattice $(k_{1},k_{2})$, but with the requirement that $k_{0}\geq 0$. In the small eccentricity limit, the relevant components of QK-Hansen coefficients are (also see Appendix A)

[TABLE]

and

[TABLE]

Using these simplified coefficients, it is straight forward to work out the dominant components of $\delta r,\delta L_{\phi}$:

[TABLE]

Similar to the Newtonian case, the $k_{1}=0,k_{2}=1$ (so that $k_{0}=1$) component of Eq. (IV.1) determine the tidally induced orbital precession $\delta n^{\prime}=n_{t}-n^{\prime}$, as it is degenerate with the radial motion of the background trajectory. As emphasized earlier, Eq. (IV.1) already neglects Post-Newtonian terms on its left hand side, which means $n^{\prime}$ is in principle not consistent with this Newtonian equation without the right-hand-side forcing term. However, as we are only interested in the leading Post-Newtonian order of $\delta n^{\prime}$, keeping Newtonian order terms in the principle part should suffice.

[TABLE]

The solution of the above equation, $\delta n^{\prime}$, is given by

[TABLE]

which reduces to Eq. (III.2) in the limit $e_{r,t,\phi}=e,K^{\prime}=1$. On the other hand, the tidal perturbation of the angular frequency, $\delta n^{\prime}_{\phi}$, can be evaluated by

[TABLE]

According to Eq. (IV.1) and Eq. (IV.1), only $k^{\prime}=0$ pieces of $\delta r,\delta L_{\phi}$ would contribute to the right hand side of the above equation. This is because $m$ is zero in these integrations for the QK-Hansen coefficients, so that the product between harmonics with frequency being multiples of $1$ and frequency $k+K^{\prime}k^{\prime}$ should be zero, unless $k^{\prime}=0$. In the small eccentricity limit, it is just

[TABLE]

with $b_{0,0}$ giving in Eq. (IV.1).

V Conclusion

In this work, we have discussed the trajectory model of an eccentric binary neutron star system, that evolves under the the influence of dynamic tidal interaction and gravitational radiation. This formalism is suitable for both Newtonian and Post-Newtonian description of the conserved dynamics. We focus more on the Newtonian description in the present study, as the leading order tidal correction can already be obtained in the Newtonian framework.

Within eccentric orbits, the direction of tidal bulges on the stars generally do not point to the companion star. This is different from circular binaries, where the tidal bulges always point to each other, even with the consideration of dynamic tide. As a result, the stars oscillatorily exchange orbital angular moment and mode angular momentum within orbital timescales. In addition, the energy and angular momentum fluxes are also modified by the beating between the orbital quadrupole moment and the star quadrupole moment. For binaries with eccentricity $e\sim 0.2$ at $f_{\rm orb}=50$ Hz, the eccentric-tidal effect on the $22$ mode radiation is only detectable by third-generation gravitational-wave detectors.

Eccentric tidal interaction also leads to tidal resonances in the inspiral stage. For circular binary, this resonance happens at $\omega_{\rm orb}=\omega_{f}/2$, which is likely higher than the merger frequency of the binary, depending on the star equation of state. For eccentric orbits, the tidal resonances show up at $\omega_{\rm orb}=\omega_{f}/k$ with $k\geq 2$, although the high order resonances are generally weaker for low-eccentricity binaries. We have analyzed the first eccentric tidal resonance, which shows up at $\omega_{\rm orb}=\omega_{f}/3$. We argue that it can be observed in the third-generation detector era. As the GW detectors are continuously improving in sensitivity, there is a growing interest to characterize the gravitational-wave radiation at the late-inspiral stage and the merger/post-merger stage, as a way to probe the neutron star physics beyond the information about tidal love number Miao et al. (2018).

In order to build a Post-Newtonian waveform model for binary neutron stars at arbitrary eccentricity, one needs to solve the matrix equation for the Fourier components of $\delta r,\delta L_{\phi}$, which in turn affects the gravitational-wave radiation of the orbit. For a frequency-domain description, additional difficulty arises in the analytical transformation from time-domain waveform to frequency-domain waveform under the stationary phase approximation, which have been discussed for eccentric binary black holes Moore et al. (2018) in the Newtonian limit. An alternative route is to use the Effective-One-Body framework for the system, for which the time-domain waveform are solved up to $1.5$ Post-Newtonian order Hinderer and Babak (2017). We shall leave the construction and validation of the eccentric binary neutron star waveform to future work.

At last, although the discussion here is presented on binary neutron star systems, the formalism is still valid for eccentric black hole-neutron star binaries. As the mass ratios for these systems are expected to be larger than those of binary neutron stars, except for low-mass black hole and neutron star binaries formed in more exotic astrophysical channels Yang et al. (2018c), it is reasonable to expect very different initial eccentricity distribution at LIGO band even for dynamically formed binaries.

Acknowledgements.

The author thank William East, Vasileios Paschalidis, Frans Pretorius, John Ryan Westernacher-Schneider for interesting discussions. This research was supported by NSERC and in part by the Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science.

Appendix A Hansen coefficients

In general, $c^{m}_{n,k}$ and $s^{m}_{n,k}$ can be obtained through numerical integration, based on Eq. (III.1) and Eq. (III.1). However, in special cases they satisfy the following relations:

[TABLE]

and

[TABLE]

where $J_{k}(x)$ is the Bessel function of the first kind. Some other useful relations are

[TABLE]

and

[TABLE]

In the small eccentricity limit, $c^{m}_{n,k},s^{m}_{n,k}$ become

[TABLE]

and

[TABLE]

Similarly, for the QK-Hansen coefficients we define in Eq. (IV.1) and Eq. (IV.1), the low eccentricity limit is given by (with $\kappa:=k+K^{\prime}k^{\prime}$)

[TABLE]

and

[TABLE]

Appendix B Precession in the large eccentricity case

Consider a 1-D problem with Hamiltonian $\frac{p^{2}}{2\mu}+V(x)$, the period is

[TABLE]

If $V(x)$ is perturbed to $V(x)+\delta V(x)$, $T$ is perturbed as

[TABLE]

Now if the relevant forcing terms in the equation of motion for $\delta r$ (the right hand side of Eq. (IV.1)) can be expressed as $\sum_{k}F_{k}\cos kl$, we have

[TABLE]

so that

[TABLE]

which gives the tidal-induced frequency shift, as $\delta T/T=-\delta n^{\prime}/n^{\prime}$.

Appendix C e-$\varpi$ representation

In the main text (c.f. Section III) the tidal perturbation is evaluated with respect to constant semi-major axis $a$, and the radial and azimuthal frequencies deviate from the Keplerian frequency. In practise, it is more convenient to discuss the modification at fixed azimuthal frequency $n_{\phi}$, as this is more suitable for constructing the frequency-domain waveforms. Therefore we shall rewrite some of the key results in Section III with respect to fixed $n_{\phi}$. For convenience we define $\varpi:=n_{\phi}$ and $x:=\varpi/\omega_{f}$.

In order to ensure constant $\varpi$, the radius has to be further shift by $\delta a$ with (c.f. Eq. (III.2))

[TABLE]

where in the last line we have extended the expression in Eq. (III.2) to include order $\mathcal{O}(e^{2})$ corrections:

[TABLE]

At constant $\varpi$, the tidally induced radius shift is

[TABLE]

and $b^{\prime}_{1}$ is now $-e\delta a$. The radial frequency shift is

[TABLE]

The energy shift at constant $\varpi$ is $\delta\mathcal{H}$ in Eq. (III.2) plus modification due to $\delta a,\delta n_{r}$, with $g_{0}/L_{\phi}=\delta a/(2a)$:

[TABLE]

Similarly the orbital-averaged modification in orbital angular momentum $\delta J_{\rm orb}$ is just $g_{0}$, which is

[TABLE]

Appendix D Evaluation of the fluxes

In the Newtonian limit, we first consider a “background” trajectory, which can be written as

[TABLE]

Adopting the convention of Poisson and Will (2014), we define the vector ${\bf n},{\bf\xi}$, which are the unit vectors $\hat{r},\hat{\phi}$. We also define $\Pi_{1}={\bf n}{\bf n},\,\Pi_{2}={\bf n}{\bf\xi}+{\bf\xi}{\bf n},\,\Pi_{3}={\bf\xi}{\bf\xi}$. We have $I_{ij}=\mu x_{i}x_{j}$, $\mathcal{I}_{ij}=I_{ij}-I\delta_{ij}/3$ and

[TABLE]

and

[TABLE]

and

[TABLE]

Therefore we have the associated fluxes for this orbital trajectory being

[TABLE]

where we have identified $n_{\phi}$ with $\varpi$ and used the expansion of $n_{r}$ in Eq. (C). Similarly, the expansion of angular momentum flux for the background trajectory is

[TABLE]

The physical trajectory is deformed from the background trajectory by $r\rightarrow r+b^{\prime}_{0}$ and $\phi\rightarrow\phi+\delta\phi$, with $\delta\dot{\phi}=g_{1}\cos l/(\mu a^{2})$. As a result, we can apply Eq. III.2 and obtain

[TABLE]

The total energy and angular momentum fluxes are approximately the simulation between orbital flux $P^{E,J}_{\rm orb}$ and $P^{E,J}_{\rm orb}$, as the modal fluxes are higher order in $\epsilon$. In the small eccentricity limit, they are

[TABLE]

and

[TABLE]

Note that Eq. (D), Eq. (D) together with Eq. (D) produce consistent result with Eq. 8 in Flanagan and Hinderer (2008) in the circular limit. In particular, when $e=0$, we have

[TABLE]

Using Eq. (C), Eq. (C) and Eq. (III.2), it is straightforward to check that

[TABLE]

or

[TABLE]

so that circular binaries remain circular even with the inclusion of dynamic tide.