Modelbased crosscorrelation search for gravitational waves from Scorpius X1
Abstract
We consider the crosscorrelation search for periodic gravitational waves and its potential application to the lowmass xray binary Sco X1. 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.25Hzwide 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 X1 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 crosscorrelation 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 parameterspace metric for the crosscorrelation 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 neutronstar 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 signaltonoise ratio due to unmodeled effects such as signal phase acceleration within the Fourier transform time scale and gradual evolution of the spin frequency.
I Introduction
The lowmass xray binary (LMXB) Scorpius X1 (Sco X1)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 X1 is presumed to be a binary consisting of a neutron star which is accreting matter from a lowmass 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) 
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).^{1}^{1}1Additionally, 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 spinup torque due to accretion is balanced by the spindown due to gravitational waves. Scorpius X1’s high xray 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 X1 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 solarsystem 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 X1 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 crosscorrelation method developed to search for stochastic GW backgrounds, treating Sco X1 as a random unpolarized monochromatic source with a known sky locationBallmer (2006).^{2}^{2}2Other 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 doublyFouriertransformed dataGoetz and Riles (2011); Aasi et al. (2014), and fitting a polynomial expansion in the Dopplermodulated 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 modelbased crosscorrelation 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 X1 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 crosscorrelation statistic using a new, streamlined formalism. Section III works out the statistical properties of the crosscorrelation 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 modelbased crosscorrelation search should compare to the directed unmodeled crosscorrelation search for a monochromatic stochastic background. Section IV considers two effects related to the dependence of the statistic on phaseevolution 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 torquebalance 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 X1.
Ii CrossCorrelation Method
The crosscorrelation 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 signalleakage issues and nuisance parameters.
ii.1 Shorttime 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 shorttime 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 finitetime 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 be^{3}^{3}3Note 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 onesided power spectral density (PSD) so that its expectation value is
(5) 
and^{4}^{4}4Strictly 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).

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.

Phaseevolution 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 gravitationalwave 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:

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

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

The relationship between the SSB time and the time at the detector depends on the sky position and time.^{5}^{5}5Specifically, 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.^{6}^{6}6Previous sensitivity estimates Dhurandhar et al. (2008); Chung et al. (2011) noted that and therefore replaced each of the finitetime 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.
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),^{7}^{7}7Note 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 .
ii.4 Construction of the crosscorrelation statistic
For a given choice of signal parameters, which determine for each SFT, and therefore for each Fourier component, it is useful to define^{8}^{8}8Note 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 crosscorrelation 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 nonnegative, 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) “noiseweighted 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 crosscorrelation term, so their would be equal to in the present notation. Similarly, our corresponds to the combination from Dhurandhar et al. (2008).^{9}^{9}9Note 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 crosscorrelation 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 crosscorrelation to some set , in particular by requiring that and . Then we define the Hermitian weighting matrix by
(35) 
so that the crosscorrelation statistic is
(36) 
Since we assume that the list of pairs includes no autocorrelations, the matrix contains no diagonal elements,^{10}^{10}10Note that if we analogously constructed the matrix to include only diagonal terms, i.e., constructed a statistic only out of autocorrelations, 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 crosscorrelation 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 crosscorrelation search to the directed stochastic search by analogy to which it was defined.
iii.1 Mean and variance of crosscorrelation statistic
The expectation value of the crosscorrelation 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 variancecovariance 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.
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 signaltonoise 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 identity^{11}^{11}11This is most easily proved by writing and using . , valid for any , the best we can do by including more bins is and therefore^{12}^{12}12Previous 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 
Since most crosscorrelation 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 crosscorrelation 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.^{13}^{13}13Note 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 singletemplate 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