Model-based cross-correlation search for gravitational waves from Scorpius X-1

Model-based cross-correlation search for gravitational waves from Scorpius X-1

John T Whelan john.whelan@ligo.org School of Mathematical Sciences and Center for Computational Relativity and Gravitation, Rochester Institute of Technology, 85 Lomb Memorial Drive, Rochester, New York 14623, USA Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), D-30167 Hannover, Germany    Santosh Sundaresan Indian Institute of Science Education and Research, Kolkata, Mohanpur Campus, Nadia District, WB 741252, India    Yuanhao Zhang yuanhao.zhang@ligo.org School of Physics and Astronomy and Center for Computational Relativity and Gravitation, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, New York 14623, USA    Prabath Peiris School of Physics and Astronomy and Center for Computational Relativity and Gravitation, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, New York 14623, USA
Wed May 20 17:21:53 2015 +0200
Abstract

We consider the cross-correlation search for periodic gravitational waves and its potential application to the low-mass x-ray binary Sco X-1. This method coherently combines data not only from different detectors at the same time, but also data taken at different times from the same or different detectors. By adjusting the maximum allowed time offset between a pair of data segments to be coherently combined, one can tune the method to trade off sensitivity and computing costs. In particular, the detectable signal amplitude scales as the inverse fourth root of this coherence time. The improvement in amplitude sensitivity for a search with a maximum time offset of one hour, compared with a directed stochastic background search with 0.25-Hz-wide bins is about a factor of 5.4. We show that a search of one year of data from the Advanced LIGO and Advanced Virgo detectors with a coherence time of one hour would be able to detect gravitational waves from Sco X-1 at the level predicted by torque balance over a range of signal frequencies from 30 to 300 Hz; if the coherence time could be increased to ten hours, the range would be 20 to 500 Hz. In addition, we consider several technical aspects of the cross-correlation method: We quantify the effects of spectral leakage and show that nearly rectangular windows still lead to the most sensitive search. We produce an explicit parameter-space metric for the cross-correlation search, in general, and as applied to a neutron star in a circular binary system. We consider the effects of using a signal template averaged over unknown amplitude parameters: The quantity to which the search is sensitive is a given function of the intrinsic signal amplitude and the inclination of the neutron-star rotation axis to the line of sight, and the peak of the expected detection statistic is systematically offset from the true signal parameters. Finally, we describe the potential loss of signal-to-noise ratio due to unmodeled effects such as signal phase acceleration within the Fourier transform time scale and gradual evolution of the spin frequency.

preprint: LIGO-P1200142-v7

I Introduction

The low-mass x-ray binary (LMXB) Scorpius X-1 (Sco X-1)Steeghs and Casares (2002) is one of the most promising potential sources of gravitational waves (GWs) which may be observed by the generation of GW detectors—such as Advanced LIGOAasi et al. (2015a), Advanced VirgoAcernese et al. (2015) and KAGRASomiya (2012)—which will begin operation in 2015 with the first Advanced LIGO observing run, and Advanced Virgo and KAGRA observations expected to follow in the coming years. Sco X-1 is presumed to be a binary consisting of a neutron star which is accreting matter from a low-mass companion; its parameters are summarized in Table 1.

Parameter Value Reference(s)
Right ascension Abbott et al. (2007a) from Bradshaw et al. (1999)
Declination Abbott et al. (2007a) from Bradshaw et al. (1999)
Distance (kpc) Bradshaw et al. (1999)
(sec) Abbott et al. (2007a) from Steeghs and Casares (2002)
(GPS sec) Galloway et al. (2014)
(sec) Galloway et al. (2014)
Table 1: Parameters of the low-mass x-ray binary Scorpius X-1. Since the sky position is determined to microarcsecond or better accuracy, the relevant astrophysical parameters with residual uncertainty are those describing the orbit. Those are the projected semimajor axis of the neutron star’s orbit, the orbital period , and the time at which the neutron star crosses the ascending node (moving away from the observer), measured in the solar-system barycenter. The orbital eccentricity of Sco X-1 is believed to be smallSteeghs and Casares (2002), and the present work presumes the orbit to be circular for simplicity; consideration of eccentric orbits would add two search parameters which are determined by the eccentricity and the argument of periapseMessenger (2011); Leaci and Prix (2015). Note that the observational constraint in Steeghs and Casares (2002) is not on itself, but on the radial velocity amplitude of the primary. We could have formulated the parameter space in terms of and rather than and , but this has no significant impact on the accuracy of the method, since the uncertainty in is dominated by that associated with . Finally, note that the orbital reference time (which we quote as the time of ascension of the compact object, cycle before the time of inferior conjunction of the companion quoted in Galloway et al. (2014)) can be propagated to a later epoch by adding an integer number of periods, at the cost of increasing the uncertainty due to the uncertainty in the period itself.

Nonaxisymmetric deformations in the neutron star can give rise to gravitational radiation, most of which is emitted at twice the rotation frequency of the neutron starJaranowski et al. (1998).111Additionally, unstable rotational modes of the neutron star, or modes Andersson et al. (1999) can lead to GW at 4/3 of the neutron star’s rotational frequency. Such deformations can be maintained by the accretion of matter onto the neutron star. It has been conjectured Bildsten (1998) that the neutron star’s rotation may be in an approximate equilibrium state, where the spin-up torque due to accretion is balanced by the spin-down due to gravitational waves. Scorpius X-1’s high x-ray flux implies a high accretion rate, which makes it the most promising potential source of observable GWs among known LMXBsWatts et al. (2008).

Since Sco X-1 is not seen as a pulsar, its rotation frequency is unknown. There is also residual uncertainty in the orbital parameters which determine the Doppler modulation of the signal, monochromatic in the neutron star’s rest frame, which reaches the solar-system barycenter (SSB). This parameter uncertainty limits the effectiveness of the usual coherent search for periodic gravitational wavesJaranowski et al. (1998). The first search for GW from Sco X-1 with the first generation of interferometric GW detectors, using data from the second LIGO science runAbbott et al. (2007a), was limited to six hours of data for this reason. A subsequent search with data from the fourth LIGO science run Abbott et al. (2007b) used a variant of the cross-correlation method developed to search for stochastic GW backgrounds, treating Sco X-1 as a random unpolarized monochromatic source with a known sky locationBallmer (2006).222Other methods have been developed, specialized to search for LMXBs. These include summing over contributions from sidebands created by Doppler modulationMessenger and Woan (2007); Aasi et al. (2015b), searching for such modulation patterns in doubly-Fourier-transformed dataGoetz and Riles (2011); Aasi et al. (2014), and fitting a polynomial expansion in the Doppler-modulated GW phasevan der Putten et al. (2010).

The stochastic analysis formed the inspiration for a new method to search for periodic gravitational waves with a model-based cross-correlation statistic which takes into account the signal model for continuous GW emission from a rotating neutron starDhurandhar et al. (2008). (This method has also been adapted Chung et al. (2011) to search for young neutron stars in supernova remnants.) The present work further develops some of the details of this method and the specifics of applying it to search for gravitational waves from Sco X-1 and, by extension, other LMXBs.

The paper is organized as follows: Section II reviews the basics of the method and the construction of the combined cross-correlation statistic using a new, streamlined formalism. Section III works out the statistical properties of the cross-correlation statistic, including the first careful determination of the effects of signal leakage and the unknown value of the inclination angle of the neutron star’s axis to the line of sight. It also considers in detail how the sensitivity of the model-based cross-correlation search should compare to the directed unmodeled cross-correlation search for a monochromatic stochastic background. Section IV considers two effects related to the dependence of the statistic on phase-evolution parameters such as frequency and binary orbital parameters: a systematic offset of the maximum in parameter space from the true signal parameters (which depends on the unknown inclination angle), and the quadratic falloff of the signal away from its maximum. The latter is encoded in a parameter space metric, which we construct in general as well as for the LMXB search both in its exact form and in limiting form relevant if the observation time is long compared to the orbital period. In Sec. V we consider limitations to the method from inaccuracies in the signal model, either due to slight variations in frequency (“spin wandering”) arising from an inexact torque-balance equilibrium, or due to phase acceleration during a stretch of data to be Fourier transformed. Finally, in Sec. VI we summarize our results and consider the expected sensitivity of this search to Sco X-1.

Ii Cross-Correlation Method

The cross-correlation method is derived and described in detail in Dhurandhar et al. (2008). In this section, we review the fundamentals, using a more streamlined formalism and including a more careful treatment of signal-leakage issues and nuisance parameters.

ii.1 Short-time Fourier transforms

Because the signal of interest is nearly monochromatic, with slowly varying signal parameters, it is convenient to describe the analysis in the frequency domain by dividing the available data into segments of length and calculating a short-time Fourier transform (SFT) from each. Since the sampling time is typically much less than the SFT duration , we can approximate the discrete Fourier transform of the data by a finite-time continuous Fourier transform. If we use the index to label both the choice of detector and the selected time interval, which has midpoint , the SFT will be333Note that the factor appears in Eq. (2.25) of Dhurandhar et al. (2008) with the wrong sign in the exponent. However, given (2) for integer , this phase correction is simply the sign so the complex conjugate does not change it.

(1)

where the frequency corresponding to the th bin of the SFT is

(2)

In practice, the data are often multiplied by a window function before being Fourier transformed, so that (1) becomes

(3)

In this work we assume that the windowing function is nearly rectangular with some small transition at the beginning and end, so that leakage of undesirable spectral features is suppressed, but the effects of windowing on the signal and noise can otherwise be ignored. The implications of other window choices are considered in Appendix A.

ii.2 Mean and variance of Fourier components

Let the data

(4)

in SFT consist of the signal plus random instrumental noise with one-sided power spectral density (PSD) so that its expectation value is

(5)

and444Strictly speaking, we should allow for data from adjacent SFT intervals in the same detector to be correlated, but we assume that the autocorrelation function falls off quickly compared to , so that we can neglect the correlation between noise in different time intervals.

(6)

If we write the noise contribution to the SFT labeled by as

(7)

then (5) implies and we can use (6) to show that

(8)

(As detailed in Appendix A, this is not the case for nontrivial windowing, where noise contributions from different frequency bins are correlated.) If we can estimate the noise PSD , we can “normalize” the data to define (as in Prix (2011))

(9)

which has mean

(10)

unit covariance

(11)

and zero “pseudocovariance”

(12)

(This is because the real and imaginary parts of each are independent and identically distributed.)

ii.3 Signal contribution to SFT

The signal from a rotating deformed neutron star is determined by various parameters of the system, which can be divided into the following categories Jaranowski et al. (1998).

  1. Amplitude parameters: intrinsic signal amplitude , the angles and which define the orientation of the neutron star’s rotation axis ( is the inclination to the line of sight and is a polarization angle from celestial west to the projection of the rotation axis onto the plane of the sky), and the signal phase at some reference time.

  2. Phase-evolution parameters: intrinsic phase evolution (frequency and frequency derivatives) of the signal, as well as parameters such as sky location and binary orbital parameters which govern the Doppler modulation of the signal.

Those parameters determine the signal received by a gravitational-wave detector at time as

(13)

where and are the antenna pattern functionsJaranowski et al. (1998); Whelan et al. (2010) which change slowly with time as the Earth rotates. The signal contribution to a SFT can be estimated by

(14)

where we have Taylor expanded the phase about the time :

(15)

The validity of this approximation will be one of the limiting factors which determines the choice of SFT duration , as detailed in Sec. V.2.

The form of (14) includes the following parameters and definitions:

  1. and depend on the inclination of the rotation axis to the line of sight.

  2. The antenna patterns and depend on the detector in question, the sidereal time at , the sky position , and the polarization angle .

  3. The relationship between the SSB time and the time at the detector depends on the sky position and time.555Specifically, if is the position of the detector and is the unit vector pointing from the source to the SSB, . Thus the phase depends on time, detector, , , sky position and—in the case of a binary—the binary orbital parameters.

The signal contribution to bin of SFT is

(16)

where we have defined

(17)

in terms of the normalized sinc function . This is plotted in Fig. 1.666Previous sensitivity estimates Dhurandhar et al. (2008); Chung et al. (2011) noted that and therefore replaced each of the finite-time delta functions with the SFT length , but a more careful treatment requires that we keep track of spectral leakage caused by the signal frequency not being centered in a SFT bin.

Figure 1: Plot of defined in (17) which determines the signal contribution to a given frequency bin of a short Fourier transform (SFT) of duration according to (16). Since the spacing between frequency bins is , there will be, for a given signal frequency , one bin whose value of lies between each pair of vertical solid lines.

The signal contribution will be largest in the th Fourier bin, defined by

(18)

whose frequency is closest to . (We have introduced the notation that is the closest integer to .) It will prove useful to define, similarly to Prix (2011),777Note that our definition of differs by a sign from the one used in Prix (2011).

(19)

where

(20)

so that . A simple search would consider, from each SFT , only the Fourier component closest in frequency to the signal frequency at the search parameters. However, as we will see, the sensitivity of the search can be improved by including contributions from additional adjacent bins, so we indicate by the set of bins to be considered from SFT , and we will construct a detection statistic using for all .

We can then write

(21)

which means that, from (10)

(22)

ii.4 Construction of the cross-correlation statistic

For a given choice of signal parameters, which determine for each SFT, and therefore for each Fourier component, it is useful to define888Note that computations can be made more efficient by use of the identity so where only the final factor depends on the bin index .

(23)

This is still normalized so that

(24a)
(24b)

where now

(25)

If we define vectors indexed by SFT number, we can write (24) and (25) in matrix form as

(26a)
(26b)
(26c)

where is the identity matrix, is a matrix of zeros, indicates the matrix transpose and the matrix adjoint (complex conjugate of the transpose).

A real cross-correlation statistic can be constructed by defining a Hermitian matrix and constructing . [Our chosen form of will be defined in (35).] Equation (26) tells us that

(27)

where the second term is a matrix with elements

(28)

where is the difference between the modeled signal phases in the two SFTs and is a geometrical factor which depends on and as follows [compare Eq. (3.10) of Dhurandhar et al. (2008)]:

(29)

where we have used the fact that the dependence of the antenna patterns can be written in terms of the amplitude modulation (AM) coefficients and as

(30)
(31)

The AM coefficientsJaranowski et al. (1998) are determined by the relevant sky position, detector and sidereal time. They can be definedPrix and Whelan (2007) as and where and are a polarization basis defined using one basis vector pointing west along a line of constant declination and one pointing north along a line of constant right ascension. Note that and are properties of the source which do not change for different SFT pairs, while and depend only on the SFT (detector and sidereal time) and sky position. It is also useful to note that the combinations

(32a)
(32b)

are independent of .

Since terms in change signs if we vary and , which are unknown, it is convenient, as proposed in Dhurandhar et al. (2008), to work with the average over those quantities, which picks out the “robust” part:

(33)

Note that is real and non-negative, while is complex. On the other hand, can be factored into , while cannot. If we define (again as in Prix (2011), but with a different overall normalization) “noise-weighted AM coefficients” and by dividing by and construct from those, we can write

(34)

or, as a matrix equation, . Note that Dhurandhar et al. (2008) did not consider issues of spectral leakage responsible for , and used a different convention for the placement of complex conjugates in atomic cross-correlation term, so their would be equal to in the present notation. Similarly, our corresponds to the combination from Dhurandhar et al. (2008).999Note that eq (3.10) of Dhurandhar et al. (2008) is also missing a factor of which should appear in . This omission was pointed out in Chung et al. (2011), but eq (5) of Chung et al. (2011) included the wrong sign in the phase correction and failed to stress that the relevant frequency is rather than .

As noted in Dhurandhar et al. (2008), an “optimal” combination of cross-correlation terms would use a weight proportional to . However, as described above, we work with in order to avoid specifying the parameters and . For reasons of computational cost to be detailed later, we limit the possible set of SFT pairs included in the cross-correlation to some set , in particular by requiring that and . Then we define the Hermitian weighting matrix by

(35)

so that the cross-correlation statistic is

(36)

Since we assume that the list of pairs includes no autocorrelations, the matrix contains no diagonal elements,101010Note that if we analogously constructed the matrix to include only diagonal terms, i.e., constructed a statistic only out of auto-correlations, the statistic would be equivalent to that used in the PowerFlux methodAbbott et al. (2008). which implies . We will later introduce, and use when convenient, the notation that labels a (nonordered) pair of SFTs .

Iii Statistics and Sensitivity

In this section we consider in detail the statistical properties of the cross-correlation statistic which were sketched in a basic form in Dhurandhar et al. (2008). In particular, we consider the impact on the expected sensitivity of spectral leakage and unknown amplitude parameters, and compare the sensitivity of a cross-correlation search to the directed stochastic search by analogy to which it was defined.

iii.1 Mean and variance of cross-correlation statistic

The expectation value of the cross-correlation statistic is

(37)

where we have used the fact that is traceless. The variance is

(38)

The first term can be evaluated by writing ; after some simplification we have

(39)

Ordinarily we would need to know something about the fourth moment of the noise distribution to evaluate the expectation value, but since contains no diagonal elements, and the different elements of are independent of each other, the expectation value can be evaluated using only the variance-covariance matrix of to give

(40)

We choose the normalization constant so that has unit variance in the limit , i.e.,

(41)

i.e.,

(42)

Written in terms of SFT pairs, the expectation value of the statistic is

(43)

Looking at (29) we see that the real part of has a piece proportional to and a piece that depends on :

(44)

The sum over SFT pairs can be broken down as a sum over detector pairs, over time offsets , and over the time stamp halfway between the time stamps of the SFTs in the pair. In an idealized long observing run, if the detector noise is uncorrelated with sidereal time, the sum over means we are averaging the two expressions and (the latter of which depends on the polarization angle ) over sidereal time. Because the former is positive definite and the latter is not, this average tends to suppress the -dependent term. This is in addition to the fact that , possibly substantially, depending on the value of , as illustrated in Fig. 2.

Figure 2: Plot of and , the coefficients of the two contributions to in (44). The factor is also equal to where is the combination of and approximately measured by the cross-correlation statistic, as shown in, e.g., Eq. (45)

If we neglect the second term in (44), (43) becomes

(45)

where

(46)

is the combination of and that we can estimate by filtering with the averaged template.

Since we have normalized the statistic so that for weak signals, the expectation value (45) is an expected signal-to-noise ratio for a signal with a given . This means that if we define a SNR threshold such that corresponds to a detection, the signal will be detectable if

(47)

iii.2 Impact of spectral leakage on estimated sensitivity

Finally, we consider the impact of the leakage factors of the form on the expectation value. Expanding these expressions, we have

(48)

If we choose only the “best bin” from each SFT, defined by (18), we have

(49)

If, instead of the best bin whose frequency is closest to , we take the closest bins to define , the sum becomes

(50)

where and are the integers below and above , respectively. Note that, because of the identity111111This is most easily proved by writing and using . , valid for any , the best we can do by including more bins is and therefore121212Previous sensitivity estimates Dhurandhar et al. (2008); Chung et al. (2011) were missing the factor of and therefore slightly overestimated the sensitivity.

(51)

The sensitivity associated with the inclusion of a finite number of bins from each SFT will depend on the value of corresponding to the signal frequency in each SFT. We can get an estimate of this by assuming that, over the course of the analysis, the Doppler shift evenly samples the range of values, and writing

(52)

with

(53)

where indicates an average over the possible offsets within the bin. We can numerically evaluate

(54)

as shown in Table 2.

1 2 3 4 5 6
Contribution 0.774 0.129 0.028 0.019 0.009 0.007
Cumulative 0.774 0.903 0.931 0.950 0.959 0.966
Table 2: Contributions to , defined in (54), from inclusion of multiple SFT bins. We see that using a single bin from each SFT leads to only around 77.4% of the maximum sensitivity given by (51), but that we can recover over 90% of this sensitivity by using two bins and over 95% by using four bins from each SFT. This table applies for rectangularly windowed data; using other window options further reduces the expected SNR, as described in Appendix A. The table also assumes that the various Doppler modulations move the signal frequency around to accomplish an average over the fractional offset of the signal frequency from the center of the bin. The validity of this approximation is explored in Sundaresan and Whelan (2012).

Since most cross-correlation searches will be computationally limited, the question of how many bins to include from each SFT is one of optimization of resources. The value of for a given , and therefore the sensitivity of the search, can be increased by including more frequency bins from each SFT, but this will involve more computations and therefore more computational resources. If instead those resources were put into a search with a larger , the value of would be higher. Naively, one might expect the computing cost to scale with the number of terms to be combined, and therefore with the square of the number of bins taken from each SFT. So increasing from to could take up to times the computing cost. On the other hand, for a fixed number of bins, we suppose that the cost will scale with the number of SFT pairs to be included times the number of parameter space points to be searched. Typical behavior will be for the density of points in parameter space to scale with for some integer value of ; as described in Sec. IV.2, for a search over frequency and two orbital parameters of an LMXB, as long as is small compared to the binary orbital period, . Since the number of SFT pairs at fixed observation time will also scale like , the overall computing cost will scale like , and quadrupling the computing time would mean multiplying the possible , and thus the number of terms in the sum (52) by . This would increase for a given by a factor of . For , this is , which is very slightly more than the benefit from including a second bin from each SFT. However, the assumption that computing cost scales like is likely an overestimate (since most of the operations can be done once per SFT rather than once per pair), so it is generally advisable to use at least two bins from each SFT.

iii.3 Sensitivity estimate for unknown amplitude parameters

The cross-correlation statistic is normalized so that and, according to (52), and now adopting the notation that refers to an unordered allowed pair of SFTs,

(55)

where is the combination of and given in (46), and is a property of the search which can be determined from noise spectra, AM coefficients, and choices of SFT pairs, without knowledge of signal parameters other than the approximate frequency and orbital parameters. Even if the noise in each data stream is Gaussian distributed, the statistic, which combines the data quadratically, will not be. It was observed in Dhurandhar et al. (2008) that each individual cross correlation between SFTs is Bessel distributed; the of the optimal sum is considered in Appendix B both in its exact form and a numerical approximation. For simplicity, in what follows we assume that the central limit theorem allows us to treat the statistic as approximately Gaussian, with mean and unit variance.131313Note that this approximation is less accurate in the tails of the distribution. Unfortunately, for a search over many independent templates, the most interesting statistic will necessarily be in the tails. For example, with templates, even a false alarm probability for the loudest statistic value would correspond to a single-template false alarm probability of . SeeZhang et al. (2015) for specific examples of this.

We consider the sensitivity estimates in Dhurandhar et al. (2008), which implicitly assume the values of and are known and used to construct the expected cross correlation used in weighting the terms in the statistic. [In our notation this would mean using rather than in the definition (35) of .] Here we perform the analogous calculation, assuming we’re using in the construction of the statistic. Thus the probability of exceeding a threshold will be

(56)

where