Fast and accurate computational tools for gravitational waveforms from binary stars with any orbital eccentricity

Abstract

The relevance of orbital eccentricity in the detection of gravitational radiation from (steady state) binary stars is emphasized. Computationally effective (fast and accurate) tools for constructing gravitational wave templates from binary stars with any orbital eccentricity are introduced including tight estimation criteria of the pertinent truncation and approximation errors.

gravitational waves, methods: data analysis, binaries: general

1 Introduction

Gravitational wave detection experiments in space including the satellite Doppler-Tracking (Bertotti & Iess 1999) and LISA (http://lisa.jpl.nasa.gov) will hopefully open a window on the low-frequency part of the gravitational wave (henceforth GW) spectrum of cosmic origin. In these frequency bands binary stars are among the most promising continuous detectable source.

A substantial fraction of binaries are expected to have orbits with non-negligible eccentricity (Barone et al. 1988; Hils et al. 1992; Pierro & Pinto 1996c) resulting into the emission of several harmonics of the fundamental orbital frequency. The importance of this fact from the standpoint of signal detection and estimation has been already noted.

For coalescing binaries, Pierro & Pinto (1996b) and Martel & Poisson (1999) pointed out that neglecting residual (albeit very small) orbital eccentricities may seriously deteriorate matched-filter detection performance. Their results, obtained in the frame of the simplest (Newtonian) Peters–Mathews (henceforth PM) model (Peters & Mathews 1963; Peters 1964; Pierro & Pinto 1996c), support the qualitative conclusion that residual orbital eccentricities cannot be bona fide disregarded in building templates for matched-filter detection of gravitational wave chirps from in-spiralling binaries.¹

For steady-state binaries with non-zero orbital eccentricity, on the other hand, using circular-orbit waveform templates, i.e. neglecting higher order harmonics, implies a potentially large loss of signal-to-noise ratio (henceforth SNR), leading to significantly worse detector's performance, as will be shown below.

The main goals of the present paper are

• (i)

  to provide some quantitative hint for validating the applicability of the simple PM model to steady-state binaries;

• (ii)

  to gauge the loss in SNR owing to the simple circular-orbit assumption and, more generally, to set some criteria for spectral waveform truncation;

• (iii)

  to introduce efficient (accurate and fast) computational tools for constructing gravitational waveform templates for (steady state) binary sources with any orbital eccentricity.

The paper is organized as follows. In Section 2 we introduce some (dimensionless) parameters whereby the applicability of the PM model to specific sources can be assessed. In Sections 3.1 and 3.2 we review the GW spectra and waveforms in the frame of the PM model. In Sections 4.1 and 4.2 we show how to evaluate the total harmonic distortion due to spectral waveform truncation, and introduce a modified Carlini–Meissel expansion tool for fast and accurate GW harmonics computation. The results in this section can be readily extended, in principle, to higher order post-Newtonian (henceforth PN) models. As an application, in Section 5 we apply our formalism to some paradigm eccentrical binary sources. Conclusions follow under Section 6. Technical developments are collected in Appendix A to C.

2 Steady-state binaries: the Peters—Mathews model

The PM model for gravitational wave emission from binary systems in a Keplerian orbit was introduced in the 1960s (Peters & Mathews 1963; Peters 1964), and recently re-examined (Pierro & Pinto 1996a). It relies on the following main assumptions: (i) point mass, (ii) weak field, (iii) slow motion, and (iv) adiabatic evolution (negligible change of the orbital parameters over each orbit). These conditions can be checked in terms of the following inequalities (Pierro & Pinto 1996a)

 

(1)

 

(2)

 

(3)

where being the orbital period, is the source gravitational radius, M_1,2 are the companion masses, and e is the eccentricity.

Tidal effects could be neglected provided neither companion star fills its Roche lobe. Following Eggleton (1983), this translates into

 

(4)

where Λ is the ratio between the physical and gravitational companion radius.²

For most steady-state binary systems, i.e. long before coalescence, ξ₁ to ξ₃ above are fairly small (see e.g. Section 5), and the PM model turns out to be perfectly adequate.

3 Steady-state binaries

3.1 Spectra

According to the PM model, the GW power radiated at the n^th harmonic of the orbital frequency by a steady binary source can be conveniently cast into the following universal form (Barone et al. 1988)

 

(5)

where the superscripts +, × refer to the fundamental GW polarization states. The spectral power distribution is embodied in the universal dimensionless functions shown in Fig. 1 for . For circular orbits only the second harmonic is emitted. The function G_max(e) plotted in Fig. 2 (left) is the ratio between the total luminosity (sum over both polarizations) of the brightest GW spectral line, and the total luminosity of a circular-orbit binary having the same χ and Δ. The brightest spectral line is the N_maxth harmonic of the orbital frequency, where N_max is a function of e only, displayed in Fig. 2 (right).

Figure 1.

versus n, for , relevant to equation (5).

Figure 2.

The functions G_max(e) (left) and N_max(e) (right), relevant to equation (5).

It is seen that for non-circular orbits, several spectral lines with comparable intensities are emitted. Thus, use of the circular orbit waveform templates implies a potentially sizeable loss in the available signal power and hence in the SNR, which can spoil the performance of the detector.

3.2 Waveforms

The far-field metric deviation (TT gauge) in the PM model is:³

 

(6)

 

(7)

where the coordinates ϑ and ϕ specify the direction of the observer in a spherical polar system where the orbit lies in the equatorial plane and the binary centre of mass is at the origin.

The metric components in equations (6) and (7) can be expanded into Fourier series under the adiabatic assumption that the orbital parameters do not change appreciably over each orbit. Hence⁴

 

(8)

 

(9)

where is a shorthand for (see Appendix A),

 

(10)

 

(11)

 

(12)

and⁵

 

(13)

For circular orbits one has simply

 

(14)

where δ_pq is the Kronecker symbol.

For steady state binaries the (Robertson) periastron advance⁶ does not produce sensible effects on the waveforms, and is thus deliberately ignored. Inclusion of the periastron advance amounts to splitting each GW spectral line into a doublet at , which cannot be resolved unless the signal is Fourier-transformed over a timespan ∼χ^2/3T s. This time is, e.g., ∼ and ∼ for PSR1534+12 and PSR1913+16, respectively.

4 Computational tools

4.1 Spectral truncation and approximation error

In order to discuss the effect of spectral truncation of equations (8) and (9) on the available SNR it is convenient to introduce the total harmonic distortion (henceforth THD)

 

(15)

where h, h˜ represent the exact and approximate values of the metric tensor, h⁽ⁿ⁾, h˜⁽ⁿ⁾ are the Fourier coefficients of h, h˜, respectively, and the L²-norms are computed by taking the time average over one orbital period of the square of the argument, within the spirit of the adiabatic approximation. If only N_T harmonics are included, then

 

(16)

It is readily recognized² that THD represents the fraction of signal power which is lost as an effect of truncation.⁷

In the most general case, where besides spectral truncation, the Fourier coefficients are computed in approximate form (as e.g. in the next subsection), one has

 

(17)

The harmonic distortions , THD_xy owing to the spectral truncation of equations (8) and (9) can be computed for any given N_T using Kapteyn's theory (Watson 1966, chapter 17) to evaluate in closed form the infinite sums in equation (17). After some lengthy but simple algebra, one obtains (see Appendix B)

 

(18)

 

(19)

 

(20)

The corresponding harmonic distortions for the TT metric components h₊, h_× can be conveniently written as follows.

 

(21)

and

 

(22)

where is the scalar product in . In order to evaluate (22) the further infinite sum

 

(23)

is needed, which is also readily obtained as explained in Appendix B.

The harmonic distortions (21) and (22) can be sensible even at very low eccentricities . Expanding (21) and (22) to lowest order in e yields

 

(24)

 

(25)

The above simple expressions are fairly accurate for , as seen, e.g., from Fig. 3, where the angular averages of the approximate and exact harmonic distortion are drawn, and seen to be almost indistinguishable and non-negligible. The (ϑ, ϕ)-dependent factor in equation (25) is plotted in Fig. 4. Its average value over the sphere is exactly equal to the (ϑ, ϕ)-independent factor in (24).

Figure 3.

Angle-averaged THD_+,× for .

Figure 4.

The angle-dependent factor in equation (25).

The obvious question is how many terms should be included in equations (8) and (9) so as to keep both THD₊ and THD_× below some specified level, for any (ϑ, ϕ).

To answer this question one may resort to the following inequalities

 

(26)

where

 

(27)

The first part of equation (26) follows immediately from equation (21); the second one is obtained from (22) using Schwartz inequality. The supremum of the function Q(ϑ, ϕ, e) occurs at , where (see Fig. 5).

Figure 5.

The function Q(ϑ,ϕ,e), equation (27).

The truncation orders required to keep , deduced from (26) are collected in Table 1.

Table 1.

Truncation orders needed to keep THD_+,× ≤0.01.

4.2 A generalized Carlini–Meissel formula

A key issue for an efficient computation of waveform-templates based on equations (10), (11) and (12) involves clever evaluation of terms such as

 

(28)

It is well known that, in general, whenever the argument and the order are close (here, in fact they are proportional through the orbital eccentricity e), numerical computation of Bessel functions either by series summation (Abramowitz & Stegun, 1968, chapter X), or by (re-normalized, downward) recurrence (Press et al. 1992, section 6.5) is inefficient. As a convenient alternative, we suggest the following generalization of the well-known (see Watson 1976, chapter XVII) Carlini–Meissel (henceforth CM) expansion

 

(29)

where (see Appendix C for the detailed deduction)

 

(30)

 

(31)

Using (31) to evaluate the Fourier coefficients h˜⁽ⁿ⁾ does not significantly spoil the accuracy of the waveforms. Indeed, spectral truncation according to Table 1 still yields THD values below 0.01.

5 Prototype sources

As an application of the above, we wrote a code for waveform template construction, and used it to compute the waveforms for several prototype sources. Taylor et al. (1993) provide data for 24 binary pulsars. In Table 2 we quote PSR 1913+16 and PSR 1534+12, as possible paradigm sources for space detectors, being respectively the most popular and closest known binary pulsars.

Table 2.

Paradigm compact binary sources.

The gravitational waveforms at computed using 8 harmonics for PSR1534+12 and 22 harmonics for PSR1913+16 (consistent with Table 1) are displayed in Figs 6 and 7, respectively. By comparison, the waveforms corresponding to are also drawn.

Figure 6.

PSR1534+12 — a waveform gallery. Circular orbit waveforms are shown as dashed lines.

Figure 7.

PSR1913+16 — a waveform gallery. Circular orbit waveforms are shown as dashed lines.

6 Conclusions

The main results in this paper can be summarized as follows. Orbital eccentricity should not be neglected in detecting gravitational waves from steady-state binaries, for which the simple Peters–Mathews model has been shown to be accurate enough. GW spectral truncation criteria have been discussed, and computationally efficient tools/techniques have been introduced for constructing reliable templates. We stress that the above tools/techniques could be readily extended, to higher order PN models with relative ease.

Acknowledgments

VP has been a Visiting Scientist at the European Space Research & Technology Centre ESTEC-ESA, under a grant from the University of Salerno; ADS, formerly at ESTEC-ESA, was a Visiting Professor at the University of Salerno in 1996. Both wish to express their appreciation to the hosting institutions.

Appendix

Appendix A: Relevant to equations (10) to (16)

In the weak-field slow-motion approximation, the Cartesian far-field harmonic-gauge metric tensor deviation components in equations (8), (9) are simply related to the source quadrupole tensor I_ij through

 

(A1)

where

 

(A2)

ρ being the companion star separation, e the eccentricity, φ the true anomaly and μ the reduced mass.

The relevant terms of the (reduced) quadrupole moment can be conveniently rewritten

 

(A3)

where a is the orbit semimajor axis,

 

(A4)

Then, using the well-known Keplerian equations (see, e.g., Watson 1976, chapter XVII)

 

(A5)

where E is the eccentric anomaly, and the relation between the latter and the mean anomaly M,

 

(A6)

one can expand ξ², η² and ξη into Fourier series of argument M, taking properly into account their parities, namely

 

(A7)

 

(A8)

 

(A9)

The relevant Fourier coefficients are readily found. Hence, using (A5) and (A6):

 

(A10)

 

(A11)

 

(A12)

Upon repeated use of trivial trigonometric identities, and in view of the integral definition of the Bessel function of the 1st kind,

 

(A13)

the Fourier coefficients (A10) to (A12) can be written

 

(A14)

 

(A15)

 

(A16)

Using (A14) to (A16) and (A7) to (A9) in (A1) to (A3) gives equations (10) to (16).

Appendix B: Relevant to equations (22)–(26)

In order to establish equations (22) to (26) one may repeatedly use the recurrency formula

 

(B1)

so as to reduce the sought series to combinations of the following (generalized) Kapteyn's expansions of the second kind

 

(B2)

 

(B3)

 

(B4)

 

(B5)

 

(B6)

 

(B7)

These latter can be summed as follows. From the Fourier analysis of Kepler motion, the following equations are readily established (see, e.g., Watson 1966, chapter 17.2)

 

(B8)

 

(B9)

 

(B10)

Differentiating equation (B8) with respect to M, and using (B10), one obtains

 

(B11)

 

(B12)

Similarly, from (B9)

 

(B13)

 

(B14)

where E is the eccentric anomaly, M the mean anomaly, and e the eccentricity. The following procedure can be then applied to equations (B8) and (B11)–(B14): (i) squaring; (ii) taking the average in M over (0,2π), using equation (B10) again; (iii) using the well-known (Euler) transformations

 

so as to express the sought series as contour integrals on of rational functions of z, which are trivially computed in terms of residues.

As an example, applying the above procedure to equation (B13), one obtains

 

(B15)

The integrand function on the right-hand side of equation (B15) has a double pole at and two simple ones at . Only two poles above fall within , and (B15) gives

 

(B16)

in agreement with Watson, chapter 17.6, equation (2). Similarly, starting from (B8), (B11), (B13) and (14) one obtains, respectively⁸

 

(B17)

 

(B18)

 

(B19)

 

(B20)

The series (B7) can be summed by differentiating both sides of equation (B16) with respect to e. Hence

 

(B21)

Appendix C: Generalized Carlini–Meissel expansions

To obtain the generalized Carlini–Meissel expansion for we start from Bessel equation for

 

(C1)

and let⁹

 

(C2)

On letting equations (C2) into (C1), we obtain

 

(C3)

then, following Carlini and Meissel, we assume that the following asymptotic representation for u_n±k, holds

 

(C4)

Substituting (C4) into (C3), and equating like powers of n (as required by consistency), we obtain

 

(C5)

 

(C6)

 

(C7)

Hence

 

(C8)

 

(C9)

Carrying out the integrations in (C2), and taking into account that we obtain

 

(C10)

 

(C11)

 

(C12)

Plugging the last three equations into equation (C2) we obtain

 

(C13)

The unknown integration constants can be found by enforcing the following obvious asymptotic equality, valid for all n

 

(C14)

Hence

 

(C15)

Hence, from (C13)

 

(C16)

 

(C17)

The right-hand side of equation (C16) above will be henceforth denoted as , and can be more conveniently written as in (29) to (31).

1

The effect of a residual (tiny) orbital eccentricity on the radiation emitted from an in-spiralling binary system was also considered in (Moreno-Garrido et al. 1994), with special emphasis on the possible relevance of periastron advance. The results in Moreno-Garrido et al. (1994) are unfortunately affected by several errors and misprints.

2

Typical values of Λ range from 10⁴ for white dwarfs down to 3 for hadronic stars. Further departures from the standard model are expected due to the possible occurrence of mass-transfer phenomena, which would be present in closely orbiting classical stars, as well as in binaries where one companion is an accreting collapsed object.

3

The GW field can also be obtained by inverting the Keplerian integral of motion relating time to the true anomaly, and exploiting the simple dependance of the radiated waveforms on this latter (Wahlquist 1987). The referred procedure is purely numerical and, to the best of our knowledge, its generalization to higher PN order models is not immediate.

4

The unknown irrelevant phase at t = 0 has been set to zero.

5

Note that for n = 1equations (10) and (11) contain Bessel functions of order −1, for which J₋₁(x)=−J₁(x).

6

The relativistic periastron advance was heuristically (i.e., inconsistently, from the post-Newtonian expansion view point) included in (Moreno-Garrido et al. 1995).

7

The THD is closely related to the fitting factor FF (Apostolatos 1996) between the exact and spectral-truncated (template) waveform. From the very definitions one gets  

8

Note that equation (3) in Watson chapter 17.6, is in error, as seen by comparison with equation (B20), and by direct numerical check. For this (erroneous) result Watson quotes (Schott 1912). The same error appears in Prudnikov et al. (1986, section 5.7.31).

9

This formula is suggested by the well-known McLaurin expansions of Bessel functions.