Gravitational waves from neutron star excitations in binary inspirals

# Gravitational waves from neutron star excitations in binary inspirals

Alessandro Parisi and Riccardo Sturani
ICTP-South American Institute for Fundamental Research, Instituto de Física Teórica (UNESP), 01140-070 São Paulo, Brazil
International Institute of Physics (IIP), Universidade Federal do Rio Grande do Norte (UFRN) CP 1613, 59078-970 Natal-RN, Brazil
Contact e-mail: riccardo@iip.ufrn.br
###### Abstract

In the context of binary inspiral of mixed neutron star - black hole systems, we investigate the excitation of the neutron star oscillation modes by the orbital motion. We study generic eccentric orbits and show that tidal interaction can excite the -mode oscillations of the star by computing the amount of energy and angular momentum deposited into the star by the orbital motion tidal forces via closed form analytic expressions. We study the -mode oscillations of cold neutron stars using recent microscopic nuclear equations of state, and we compute their imprint into the emitted gravitational waves.

###### keywords:
Neutron stars, gravitational wave sources, relativistic star oscillations

## 1 Introduction

After the historical detections of gravitational waves by binary black holes Abbott et al. (2016), it is expected that mixed binaries composed of a neutron star (NS) and a black hole (BH) may be the next, qualitatively different type of source to be detected in the gravitational wave (GW) channel. At first approximation mixed NS-BH can be treated in General Relativity (GR) on equal footing as binary BH systems, however the presence of matter in the GW source may lead to new detectable astrophysical effects in the GW signal that are not expected to appear in the binary BH case like e.g. NS tidal deformations leaving an imprint in the GW signal Bildsten & Cutler (1992); Flanagan & Hinderer (2008) and breaking of the NS giving origin to a gamma ray burst or more general electromagnetic counterpart Lattimer & Schramm (1976), to name only the most studied effects.
Beside their direct phenomenological relevance, these effects carry information on the highly uncertain equation of state of the NS, thus making GW detection an invaluable probe of the internal structure of NSs. In this work we focus on a specific effect in GW signals: NS can be tidally deformed by the orbital motion in generic elliptic orbits, hence setting oscillations of the NS normal modes. The orbit being elliptical can induce resonant oscillations at a frequency much higher than the frequency scale set by the inverse of the orbital period, since in general NS oscillations are much higher than orbital frequency of inspiral binary systems.
Quantifying this phenomenon in light of the exciting prospect of a future GW detection has been the subject of extensive investigations in literature in a number of different contexts. The theoretical setup for studied such tidally induced NS oscillations has been provided in Thorne (1969); Press & Teukolsky (1977). In Fabian et al. (1975) it was originally proposed that tidal encounters between a NS and a main-sequence star might lead to the formation of X-ray binaries in globular clusters. In Shibata (1994) the effects of the tidal resonances for a circular orbital motion has been studied, with the result if the companion of a NS is a BH of mass , the -mode resonance is unimportant, while the -mode resonance may affect the orbital evolution just before the merging. Rathore et al. (2005) considered the energy absorbed by tidal excitations in eccentric orbit (but not their imprint in the GW-form). Reisenegger & Goldreich (1994) compute the effect on the emitted GW phase of resonant mode excitation by the circular inspiral motion. Rotating NS we considered by Ho & Lai (1999) (including g-modes and r-modes) when the spin axis is aligned or anti-aligned with the orbital angular momentum axis. Carter & Luminet (1983) solved for the tidal deformation dynamics of a NS in an external field of a massive object and recently Chirenti et al. (2017) presented a framework for the discussion of binary NS and mixed NS-BH ones oscillation mode excitation and detection via the GWs observed by future GW detector as Einstein Telescope or Cosmic Explorer. Numerical results on the GW emission of tidally excited NS oscillations in the last stages of a coalescence have been given in Gold et al. (2012), and in Steinhoff et al. (2016) the imprint of resonant tidal on the gravitational waveform has been computed within the effective one body description of the two body orbital motion.
In the present paper we consider non-rotating NS with four different equations of states Akmal et al. (1998); Douchin & Haensel (2001); Walecka (1974); Bethe & Johnson (1974) with the goal of translating resonant excitations of various -modes for NSs inspiraling binary NS-BH systems that move in an elliptical orbit into quantitative prediction for the emitted GW-form.
Numerical simulations show that most of the energy released in gravitational waves is indeed transferred into -modes, which are characterized by a wave-function free of nodes along the radial direction. We do not study the possibilities of exciting the g-modes because these modes are related to the presence of density discontinuities in the outer envelopes of NSs, see Finn (1987) and Strohmayer (1993), density discontinuities in the inner core as a consequence of phase transitions at high density, as studied in Sotani et al. (2001), and/or thermal gradients as for a proto-NS, see e.g. Ferrari et al. (2003). In this paper we do not consider the possibility of having discontinuities of the density, moreover we focus on barotropic equations of state where the pressure depends only on the energy density, implying that all g-modes degenerate to zero frequency, hence we focus on the excitations of -modes. Our study is based on the following simplifying assumptions:
(i) we neglect BH rotation, thus we treat the BH as a point particle with mass ; (ii) the hydrodynamic stability of NS is computed using the Oppenheimer-Volkoff equations, but we use Newtonian equations to calculate the oscillation modes, see Appendix A; (iii) the NS does not rotate and we neglect viscous effects.
By implementing the formalism presented in Thorne (1969); Press & Teukolsky (1977) we find generic analytic expressions for the energy and angular momentum deposited into NS oscillations during the elliptic orbital motion, allowing to compute the mass quadrupole which is sourcing GW emission, and eventually comparing it with the orbital quadrupole.
The outline of this paper is as follows: in Sec. 2 we present the setup of the physical system under consideration, and we provide new analytic expressions for the dynamics of tidally induced NS oscillations, which are the main result of this paper. In Sec. 3 we analyze quantitatively their GW emission. Finally, conclusions for future detectability of NS oscillations in the GW channel are drawn in Sec. 4. We set the speed of light throughout this paper.

## 2 Coupling of neutron star oscillation modes to orbital motion

In this section we study the tidal excitation of NS oscillation modes in non-rotating stars in an elliptical orbit. Our analysis will be general, but the astrophysical case we have in mind is that of a binary NS-BH system. The idea to compute the energy deposited in stellar oscillations by the tidal gravitational field is first described by Turner (1977) and Press & Teukolsky (1977).
In this paper we use Newtonian linearized equations to calculate the oscillation modes. The use of Newtonian equations is consistent with our Newtonian description of tidal interactions. For the -mode, general relativistic effects are expected to modify our results of oscillation frequencies by not more than per cent, see Lai (1994), where and are the mass and radius of the NS. We also neglect the spin of the NS. When , the normal modes of the star get more complicated, especially when becomes comparable to the mode frequencies Gaertig & Kokkotas (2008). For the eigenmodes can be adequately approximated by those of a non-rotating spherical star, the basic equations that governing the oscillations of stars are discussed in more detail in Appendix A.
The NS oscillations are excited by tidal forces while the NS is bound in a binary system with black hole in an eccentric orbit whose evolution is driven by gravitational radiation. The distance between two objects in an elliptic orbit can be parametrized by, see e.g. eq. (4.54) of Maggiore (2008),

 D=a(1−e2)1+ecosψ (1)

being the semi-major axis and the eccentricity (with corresponding to the periastron), and the true anomaly is related to the eccentric anomaly and time via, see e.g. eqs. (4.57,58) of Maggiore (2008),

 β≡u−esinu=ω0t,cosψ=cosu−e1−ecosu, (2)

being the orbital period, with the following relationships holding among orbital parameters

 ˙ψ=[GNMa(1−e2)]1/2D2, (3)

(where is the total mass of the binary system and the Newton constant) and the standard definition of the relativistic orbital parameter

 x≡(GNMω0)2/3=GNMa, (4)

the last equality holding only at Newtonian level.
In order to study quantitatively the effect of the gravitational force inducing oscillations into the NS and following the procedure outlined in Press & Teukolsky (1977), it is useful to expand the Newtonian potential in spherical harmonics, see e.g. eq. (3.70) of Jackson (1998), centered at the star as per

 1|D−r|=∞∑ℓ=0ℓ∑m=−ℓ4π2ℓ+1rℓDℓ+1Y∗ℓm(θ,ϕ)Yℓm(π/2,ψ), (5)

being coordinates of the mass elements of the NS, are the spherical harmonic indices and the orbital motion is assumed to be planar (no spin-induced precession). Using eq. (5) for elliptic orbit, it will be useful to expand for generic into a Fourier series of the type

 eimψDℓ+1=1aℓ+1∞∑j=0{c(ℓ,m)j(e)cos(jβ)+is(ℓ,m)j(e)sin(jβ)}. (6)

The detailed calculation of the Fourier coefficients and their analytic expressions are presented in Appendix B.
In order to perform an analytic quantitative analysis we borrow here the framework of Rathore et al. (2005), where NS oscillations are modeled as a series of damped harmonic oscillator displacements driven by external force, that we can take purely monocromatic:

 ¨xn(t)+2˙xn(t)τn+ω2nxn(t)=Cjcos(ωjt)+Sjsin(ωjt), (7)

where is the stellar mode frequency, its damping time111As a possible mechanism for the damping of non-radial NS oscillations we take the gravitational emission, we do not consider neutrino losses, radiative heat leakage, and magnetic damping., is the -th harmonic of the main orbital angular frequency , and the exciting force amplitude.222Note that the time scale of variation is set by the GW radiation and via the Einstein quadrupole formula (as , with ), hence we neglect the time variation of the frequency of “forcing” term in eq. (7). Eq. (7) admits the exact analytic solution

 [(ω2j−ω2n)2+4ω2j/τ2n]xn(t)==(ω2n−ω2j)(Cjcos(ωjt)+Sjsin(ωjt))+2ωj/τn(Cjsin(ωjt)−Sjcos(ωjt)), (8)

the solution to the homogeneous equation being

 x(h)n∝e−t/τncos[(ω2n−1/τ2n)1/2t+ϕ0], (9)

leading to an average absorbed energy per unit of mass per unit of time

 ˙E=(C2j+S2j)ω2j/τn(ω2j−ω2n)2+4ω2j/τ2n. (10)

The NS oscillation vectors satisfy an equation of the type see Kosovichev & Novikov (1992)

 (ρd2dt2+L)→ζ(t,→r)=−ρ→∇U(→r), (11)

where is an operator characterizing the internal restoring force of the star. In order to apply this toy model of a damped harmonic oscillator to the tidally excited NS oscillation, we decompose the oscillation field into normal modes with factorized time and space dependence:

 →ζ(t,→r)=∑n,ℓ,mqnℓm(t)→ξnℓm(→r), (12)

where we have added the spherical harmonics labels and the spatial mode eigenfunctions satisfy

 (L−ρω2n)→ξnℓm=0, (13)

allowing the identification of with the stellar frequency of the eigenmode. The differential equations the oscillation modes fields satisfy are summarized in Appendix A, which are solved for 4 different equations of state and 4 values of the central density of the NS, with the resulting mass, radius, frequency and damping times (the last two depending on ) are reported in Appendix C for .
It is also useful to expand the eigenmodes into a radial () and a poloidal () component

 →ξnℓm(→r)=(ξ(r)nℓ(r)^er+rξ(h)nℓ(r)→∇)Yℓm(θ,ϕ), (14)

and impose the normalization condition333Note that with the normalization chosen have dimension of length, is dimension-less. However the normalization can be arbitrarily chosen without affecting physical results, our choice has the advantage of making following formulae simpler.

 ∫d3xρ(r)→ξ∗nℓm⋅→ξn′ℓ′m′=∫drr2ρ(r)(ξ(r)nℓξ(r)n′ℓ′+ℓ(ℓ+1)ξ(h)nℓξ(h)n′ℓ′)δℓ,ℓ′δm,m′=ρ0R5∗δn,n′δℓ,ℓ′δm,m′, (15)

where are respectively the density, central density and radius of the NS, and we used

 ∫dΩYℓm(θ,ϕ)Y∗ℓ′m′(θ,ϕ)=δℓ,ℓ′δm,m′∫dΩr2→∇Yℓm(θ,ϕ)⋅→∇Y∗ℓ′m′(θ,ϕ)=ℓ(ℓ+1)δℓ,ℓ′δm,m′,∫dΩ→r⋅→∇Yℓm(θ,ϕ)→r⋅→∇Y∗ℓ′m′(θ,ϕ)=δℓ,ℓ′δm,m′, (16)

and the integral of products of spherical harmonics with unequal number of derivatives vanish for any , .
By multiplying both members of eq.(11) by , substituting the expansion in eq. (5), and integrating over the NS volume the mode is singled out and it satisfies an equation of the type (7):

 ¨qnℓm(t)+2τnℓ˙qnℓm(t)+ω2nqnℓm(t)=GNMBHa3(R∗a)ℓ−2QnℓWℓm×∑j(c(ℓ+1,m)j(e)cos(jβ)+is(ℓ+1,m)j(e)sin(jβ)) (17)

where

 Wℓm≡4π2ℓ+1Yℓm(π/2,0),Qnℓ≡1ρ0Rℓ+3∗∫R∗0drr2ρ(r)ℓrℓ−1(ξ(r)nℓ+(ℓ+1)ξ(h)nℓ), (18)

is the black hole mass. Note that the r.h.s of eq.(17) is complex, but given the symmetries of the coefficients: (and if have different parity), , the sum of returns a real quantity. The modes thus satisfy an equation of the type (7) with the coefficients replaced by

 (Cj,Sj)→GNMBHa3(R∗a)ℓ−2QnℓWℓm(c(ℓ,m)j(e),s(ℓ,m)j(e)). (19)

These expressions will be needed in sec. 3 to compute the time varying quadrupole associated to these oscillations, source of GWs.
The rate of energy (per unit of mass, per unit NS radius) absorbed by each oscillation modes can be read from eq. (10) by inserting the above values of , summing over and , the rate of absorbed energy via tidal mechanism being

 ˙E∗=∑j˙Ej=ρ0R∗(R∗a)4(GNMBHa)2∑j,n,ℓ,m(c(ℓ,m)j2+s(ℓ,m)j2)(R∗a)2ℓ−4×Q2nℓW2ℓmω2j/τnℓ(ω2j−ω2nℓ)2+4ω2j/τ2nℓ. (20)

The contribution from individual modes to the rate of energy absorption is plotted in fig. 1 after being divided by the factor

 K≡ρ0R∗(R∗a)4(GNMBHa)2(GNMω0)2ω02≃1.5⋅10−14M⊙sec(x0.01)9(ρ01015gr/cm3)(R∗10Km)5(MBH4M⊙)2(M6M⊙)−6, (21)

where for , . Factorizing the absorbed energy rate by the quantity has the virtue of making dimension-less and independent on the relativistic parameter (as long as the orbital frequency does not hit a resonance with ) and mildly dependent on , .

In fig. 2 we report the absorbed energy rate normalized by

 ˙EGW0≡325GNη2x5 (22)

(being the reduced mass of the orbital system), which is the expression of the leading order in of the GW emission rate at zero eccentricity from a binary inspiral, making visually easier the comparison between GW radiated energy and . For we use the 3PN formula taken from Arun et al. (2008), see also sec. 10.3 of Blanchet (2014).

The absorbed angular momentum can be computed in a similar way, following Lai (1994), where it is noted that the variation of angular momentum

 ˙L∗=−∫d3x(ρ0+δρ)(^z⋅→r×→∇U) (23)

we can derive in our setup

 ˙L∗=∑nℓqnℓm(t)∫d3x→∇⋅(ρ0→ξnℓm)∂U∂ψ=∑nℓmqnℓm(t)∫d3x→∇⋅(ρ0→ξnℓm)GNMBHa(ra)ℓWℓmim×Y∗ℓm(θ,ϕ)∑j(c(ℓ,m)j(e)cos(jβ)+is(ℓ,m)j(e)sin(jβ)), (24)

where in the last passage we have inserted the expansion of eqs. (5,6) and derived by parts inside the integral. In this form the angular momentum absorption rate by NS oscillations can be rewritten as:

 ˙L∗=ρ0R∗(GNMBHa)22∑j,n,ℓ,m>0mc(ℓ,m)js(ℓ,m)j(R∗a)2ℓ×Q2nℓW2ℓmωj/τnℓ(ω2j−ω2nℓ)2+4ω2j/τ2nℓ. (25)

In fig. 3 the absorbed angular momentum rate normalized by the leading order expression in of is reported for various values of the relativistic parameter and the eccentricity . The values of are negligible with respect to and given the typical moment of inertia of a NS ( gr cm, see book of Haensel et al. (2007)), the induced rotation on the NS is also negligibly small.

## 3 Gravitational Wave emission

We have seen in the previous section that the energy absorbed in by the NS is very small compared to the orbital energy at moderate eccentricity values (), hence such absorption will not alter in any significant way the chirping signal. However the energy absorbed will set oscillations in the neutron star that gives rise to a time varying quadrupole, which will in turn generate GWs with a significantly different pattern that the GWs associated to the decaying orbital motion.
The general expression for the GW in the TT gauge is given by, see e.g. eq. (3.275) of Maggiore (2008),

 hTTij(t,r)=1rGN+∞∑ℓ=2ℓ∑m=−ℓ[uℓm(TE2ℓm)ij+vℓm(TB2ℓm)ij] (26)

where () is linearly related to the -th time derivative of the mass (momentum) multipole moments. The leading-order contribution to radiation reaction comes from the mass quadrupole term, for which it is (see e.g. sec. 3 of Maggiore (2008))

 u2m=1615π√3¨QijY2mij∗, (27)

being the tensor spherical harmonics and is the standard quadrupole mass moment in Cartesian coordinates. It will be convenient to express the leading order GW amplitude in terms of the spherical components of the quadrupole, related to their Cartesian counterpart via

 Q2m≡8π15Qij(Y2mij)∗,Qij=∑|m|≤2Q2mY2mij, (28)

leading to (explicit expressions of tensor spherical harmonics are reported in app. D)

 u2m=2√3¨Qm. (29)

We now have all the ingredients to relate the leading GW source to the NS tidal oscillations via

 Q∗2m=8π15∫ρr2Y∗2md3x, (30)

that in terms of the displacement vector introduced in eqs. (11,12) can be expressed as, see Ushomirsky et al. (2000), by

 158πQ∗2m=∫(ρ0+δρ)r2Y∗2md3x=−∑n∫→∇⋅(ρ0→ζn2m)r2Y∗2md3x≃q02m(t)(2∫R∗0ρ0{ξ(r)02+3ξ(h)02}r3dr−ρ0ξ(r)02r4∣∣R∗0), (31)

where an integration by parts has been performed in the last step, the explicit expression of has been inserted and only the contribution has been considered. since we analyzed only the -mode. Observing that the boundary term is numerically smaller than the integral term, substituting the solution of eq. (17) and considering only the resonant contribution for the NS average quadrupole value can be written as

 ⟨Q2∗22⟩1/2≃2√2π15ρ0R5∗Q202W22GNMBHa3τ02ω¯j[(c(2,2)¯j)2+(s(2,2)¯j)2]1/2≃4√2π15(ρ0R5∗τ02)1/2ω02Q02√˙E(ℓ=2)¯j≃10−2M⊙km2(ρ01015gr/cm3)1/2(R∗10km)5/2×⎛⎜⎝˙E(l=2)j10−8M⊙/sec⎞⎟⎠1/2(τ020.1sec)1/2(ω0218kHz)−1. (32)

The quantity directly related to GW emission, , follows straightforwardly via eq. (29). In fig. 4 we report the contribution to the second time derivative of the quadrupole (divided by the reduced mass of the binary system) and as a comparison the (magnified) second derivative of the quadrupole associated to , NS oscillations during an ordinary binary inspiral in which the orbit shrinks due to GW back reactions.

For comparison, we also report in fig. 5 the time evolution of the displacement along the inspiral phase.

## 4 Conclusions

In this paper we have developed and presented a framework able to perform analytic and quantitative study of the excitations of a neutron star in an inspiralling binary system of arbitrary eccentricity. We have computed the energy and the angular momentum deposited into stellar mode oscillations by the tidal field via closed form analytic formulae. The amount of energy absorbed by the neutron star in a given mode depends on the overlap of the tidal force field with the displacement field of the mode, hence it requires solving the equilibrium equations of a neutron star, done here in the Newtonian approximation. We focused our analysis on the fundamental -mode of a non-relativistic star, finding the rate of energy absorbed and angular momentum as a function of eccentricity and of the period of the inspiral orbital, when -mode can be in resonance with higher harmonics of the main orbital frequency.

As a future development of this work, we intend to extend our analysis to the General Relativistic equilibrium equations of a rotating neutron star, with the inclusion of -mode and -modes, and considering a not barotropic equation of state: such modes have lower frequency values than the -mode, and can therefore be excited at resonance in an elliptical orbit earlier in the inspiral phase. The phenomenological impact of the computations presented here relies on the signature that neutron star oscillations will imprint onto the gravitational signals of an inspiral binary system. Despite being sub-dominant with respect to the gravitational wave sourced by the orbital motion, the detailed features of the star oscillation bears invaluable information on its equation of state and density, allowing to make a bridge to the nuclear physics ruling its equilibrium. Since it is expected in the near future that third generation gravitational wave detector could observe signals from binary systems involving neutron star at signal-to-noise ratio of order or more, see e.g. Punturo et al. (2010), and that such detection will involve the observation of hundred of thousand gravitational wave cycles during the inspiral of a binary system for a time stretch of order of several days, the quantitative prediction of the modification of the inspiral signal, even at very low level, will have an impact on the physics outcome of the detection.

## Acknowledgments

The authors wish to thank C. Chirenti for useful discussions. The work of AP has been supported by the FAPESP grant 2016/00096-6, RS has been supported by FAPESP grant 2012/14132-3.

## References

• Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev., X6, 041015
• Akmal et al. (1998) Akmal A., Pandharipande V. R., Ravenhall D. G., 1998, Phys. Rev., C58, 1804
• Arun et al. (2008) Arun K. G., Blanchet L., Iyer B. R., Qusailah M. S. S., 2008, Phys. Rev., D77, 064035
• Balbinski & Schutz (1982) Balbinski E., Schutz B. F., 1982, Mon. Not. R. ast. Soc., 200, 43P
• Bethe & Johnson (1974) Bethe H. A., Johnson M. B., 1974, Nucl. Phys., A230, 1
• Bildsten & Cutler (1992) Bildsten L., Cutler C., 1992, ApJ, 400, 175
• Blanchet (2014) Blanchet L., 2014, Living Rev. Rel., 17, 2
• Carter & Luminet (1983) Carter B., Luminet J.-P., 1983, Astro. Astrophys, 121, 97
• Chirenti et al. (2017) Chirenti C., Gold R., Miller M. C., 2017, Astrophys. J., 837, 67
• Douchin & Haensel (2001) Douchin F., Haensel P., 2001, Astron. Astrophys., 380, 151
• Dziembowski (1971) Dziembowski W. A., 1971, Acta Astronomica, 21, 289
• Fabian et al. (1975) Fabian A. C., Pringle J. E., Rees M. J., 1975, Mon.Not.R.Astr.Soc, 172, 15p
• Ferrari et al. (2003) Ferrari V., Miniutti G., Pons J. A., 2003, Mon. Not. Roy. Astron. Soc., 342, 629
• Finn (1987) Finn L. S., 1987, Mon.Not.R.Astr.Soc, 227, 265
• Flanagan & Hinderer (2008) Flanagan É. É., Hinderer T., 2008, Phys. Rev., D77, 021502
• Gaertig & Kokkotas (2008) Gaertig E., Kokkotas K. D., 2008, Phys. Rev., D78, 064063
• Gold et al. (2012) Gold R., Bernuzzi S., Thierfelder M., Brugmann B., Pretorius F., 2012, Phys. Rev., D86, 121501
• Haensel & Zdunik (2008) Haensel P., Zdunik J. L., 2008, Astron. Astrophys., 480, 459
• Haensel et al. (2007) Haensel P., Potekhin A. Y., Yakovlev D. G., 2007, Neutron stars 1: Equation of state and structure. Springer
• Ho & Lai (1999) Ho W. C. G., Lai D., 1999, Mon. Not. Roy. Astron. Soc., 308, 153
• Jackson (1998) Jackson D., 1998, Classical Electrodynamics. Wiley
• Kosovichev & Novikov (1992) Kosovichev A. G., Novikov I. D., 1992, MNRAS, 258, 715
• Lai (1994) Lai D., 1994, "Mon. Not. Roy. Astron. Soc.", 270, 611
• Lattimer & Schramm (1976) Lattimer J. M., Schramm D. N., 1976, ApJ, 210, 549
• M.Abramowitz & Stegun (1964) M.Abramowitz Stegun I., 1964, Handbook of mathematical functions. Dover
• Maggiore (2008) Maggiore M., 2008, Gravitational Waves. Oxford University Press
• Press & Teukolsky (1977) Press W., Teukolsky S., 1977, Astrophys. J., 213, 183
• Punturo et al. (2010) Punturo M., et al., 2010, Class. Quant. Grav., 27, 194002
• Rathore et al. (2005) Rathore Y., Blandford R. D., Broderick A. E., 2005, Mon. Not. Roy. Astron. Soc., 357, 834
• Reisenegger & Goldreich (1994) Reisenegger A., Goldreich P., 1994, Astrophys. J.
• Shibata (1994) Shibata M., 1994, Prog. Theor. Phys., 91, 871
• Sotani et al. (2001) Sotani H., Tominaga K., Maeda K.-I., 2001, Phys. Rev. D, 65, 024010
• Steinhoff et al. (2016) Steinhoff J., Hinderer T., Buonanno A., Taracchini A., 2016, Phys. Rev., D94, 104028
• Strohmayer (1993) Strohmayer T. E., 1993, The Astrophysical Journal, 417, 273
• Thorne (1969) Thorne K. S., 1969, Astrophys. J., 158, 997
• Turner (1977) Turner M., 1977, ApJ, 216, 914
• Ushomirsky et al. (2000) Ushomirsky G., Cutler C., Bildsten L., 2000, Mon. Not. Roy. Astron. Soc., 319, 902
• Vallisneri (2000) Vallisneri M., 2000, Phys. Rev. Lett., 84, 3519
• Walecka (1974) Walecka J. D., 1974, Annals Phys., 83, 491

## Appendix A Four first-order linear differential equations of non-radial oscillations

The normal modes of a spherical star can be labeled by spherical harmonic indices and , and by a “radial quantum number" . In spherical coordinates the Lagrangian displacement of a fluid element is given by

 ξnℓm=⎡⎣ξ(r)nℓ(r),ξ(h)nℓ(r)∂∂θ,ξ(h)nℓ(r)sinθ∂∂ϕ⎤⎦Yℓm(θ,ϕ)eiσt (33)

where denotes a spherical harmonic; and denotes the pulsation angular frequency. The oscillation is assumed to be adiabatic, we ignore the thermal evolution of the NS, for simplicity we use the Newtonian description in the Dziembowski (1971) formulation, in this case the equations reduce to a system of four first-order differential equations with four dimensionless variables, given by:

 y1=ξ(r)nℓr,y2=1gr(p′ρ+Φ′)=σ2gξ(h)nℓ, (34)
 y3=Φ′gr,y4=1gdΦ′dr, (35)

Here, the meanings of the symbols are as follows: and are the radial part of the Eulerian perturbation to the pressure and the gravitational potential , respectively; is the distance from the center of the star, is the density, and is the local acceleration due to gravity. The system of differential equations that governs the linear adiabatic oscillations of stars is then given by:

 rdy1dr = (Vg−1−ℓ)y1+[ℓ(ℓ+1)c1ω2−Vg]y2+Vgy3 (36) rdy2dr = (c1ω2−A∗)y1+(3−U+A∗−ℓ)y2−A∗y3 (37) rdy3dr = (3−U−ℓ)y3+y4 (38) rdy4dr = A∗Uy1+UVgy2+[ℓ(ℓ+1)−UVg]y3−(U+ℓ−2)y4 (39)

Where

 Vg=−1Γ1dlnpdlnr=grc2s,A∗=1Γ1dlnpdlnr−dlnρdlnr,U≡dlnm(r)dlnr=4πρr3m(r), (40)
 c1≡r3R3∗M∗m(r),Γ1=(∂lnp∂lnρ)S,ω2=R3∗GNM∗σ2. (41)

Here is the first adiabatic exponent, is the sound speed, is the concentric mass, and are the total mass and radius of the star, respectively, and is the gravitational constant. There are four boundary conditions, the inner boundary conditions at are:

 {c1ω2y1−ℓy2=0ℓy3−y4=0,

the outer boundary conditions at are:

 {y1−y2+y3=0(ℓ+1)y3+y4=0.

The two central boundary conditions require that the two divergences involved, , , remain finite. At the surface we require to be finite and , the gravitational force per unit mass, to be continuous across the perturbed surfaces. The above equations and boundary conditions constitute an eigenvalue problem for the eigenvalue .
The expression for the damping time due to emission of gravitational waves in the Newtonian case see Thorne (1969); Balbinski & Schutz (1982) ) is given by:

 τnℓ≡(ℓ−1)[(2ℓ+1)!!]2ℓ(ℓ+1)(ℓ+2)(σ2πG)(cσ)2ℓ+1∫R∗0drρr2[ξ(r)nℓ(r)2+ℓ(ℓ+1)ξ(h)nℓ(r)2]{∫R∗0drρrℓ+1[ξ(r)nℓ(r)+(ℓ+1)ξ(h)nℓ(r)]}2 (42)

where .

## Appendix B Expansion of the Fourier coefficients

Expanding in eq. (6) we have

 c(ℓ,m)j(e)=cjπ(1−e2)ℓ+1∫π−πcos(mψ)(1+ecosψ)ℓ+1cos(jβ)dβ, (43)

for , where we used that is an odd function of time, hence ()is an even (odd) function of time, for , and .
In order to expand into sums of terms of the type it is useful to express it in terms of powers of via M.Abramowitz & Stegun (1964)

 cos(mθ)=Tm(cos(θ)), (44)

where is the Chebyshev polynomial of order and it has the form

 Tm(x)=[m/2]∑k=0t(s)rxm−2k, (45)

begin the integer part of . Using the standard relationships between eccentric anomaly , true anomaly and time , see sec. 2, one finds

 1+ecosψ=1−e21−ecosu,dβ=(1−ecosu)du, (46)

to obtain

 c(ℓ,m)j=2cnπ∫π0[m/2]∑k=0t(k)m(cosu−e)m−2k(1−ecosu)m−2k+ℓcos(ju−jesinu)du. (47)

In order to perform this integral we use the standard Taylor-expansions

 (1−x)n=n∑k=0(−1)kn!k!(n−k)!xk,1(1−x)n=∞∑k=0(n+k−1)!k!(n−1)!xk, (48)

to write

 c(ℓ,m)j(e)=2cjπ∫π0[m/2]∑k=0t(k)m(−e)m−2km−2k∑p=0(m−2k)!p!(m−2k−p)!(−1)p(cosue)p∞∑n=0(m−2k+ℓ+n−1)!n!(m−2k+ℓ−1)!(ecosu)ncos(ju−jesinu)du=2cjπ∫π0[m/2]∑k=0m−2k∑p=0∞∑n=0(−1)p+mt(k)m(m−2k)!p!(m−2k−p)!(m−2k+ℓ+n−1)!n!(m−2k+ℓ−1)!em−2k+n−p(cosu)p+ncos(ju−jesinu)du, (49)

and then we use the DeMoivre formula

 cosn(u)=12nn∑k=0n!k!(n−k)!