III.3 \mathcal{F}-statistic

Implementation of the frequency-modulated sideband search method for gravitational waves from low mass X-ray binaries


We describe the practical implementation of the sideband search, a search for periodic gravitational waves from neutron stars in binary systems. The orbital motion of the source in its binary system causes frequency-modulation in the combination of matched filters known as the -statistic. The sideband search is based on the incoherent summation of these frequency-modulated -statistic sidebands. It provides a new detection statistic for sources in binary systems, called the -statistic. The search is well suited to low-mass X-ray binaries, the brightest of which, called Sco X-1, is an ideal target candidate. For sources like Sco X-1, with well constrained orbital parameters, a slight variation on the search is possible. The extra orbital information can be used to approximately demodulate the data from the binary orbital motion in the coherent stage, before incoherently summing the now reduced number of sidebands. We investigate this approach and show that it improves the sensitivity of the standard Sco X-1 directed sideband search. Prior information on the neutron star inclination and gravitational wave polarization can also be used to improve upper limit sensitivity. We estimate the sensitivity of a Sco X-1 directed sideband search on 10 days of LIGO data and show that it can beat previous upper limits in current LIGO data, with a possibility of constraining theoretical upper limits using future advanced instruments.

95.85.Sz, 97.60.Jd, 97.80.Jp

I Introduction

Rotating neutron stars are the main candidates for sources of persistent, periodic gravitational radiation detectable by ground-based, long-baseline gravitational wave interferometers. Time varying quadrupole moments in (and thus gravitational-wave emission from) these sources can result from deformations of the solid crust (and possibly a solid core) supported by elastic stresses Bildsten (1998); Ushomirsky et al. (2000a); Owen (2005); Haskell et al. (2007); Lin (2007); Johnson-McDaniel and Owen (2013), deformations of various parts of the star supported by magnetic stresses Bonazzola and Gourgoulhon (1996); Cutler (2002); Melatos and Payne (2005); Vigelius and Melatos (2008); Haskell et al. (2008); Mastrano and Melatos (2012), or free precession Jones and Andersson (2002); Jones (2012) or long-lived oscillation modes of the entire star Owen et al. (1998); Andersson et al. (1999); Bondarescu et al. (2007); Haskell et al. (2012); Bondarescu and Wasserman (2013).

Neutron stars in accreting binary systems are an important sub-class of periodic gravitational wave sources. Accretion may trigger or enhance the aforementioned gravitational-wave emission mechanisms, creating or driving the quadrupole moment toward its maximum value through thermal, magnetic, or other effects Bildsten (1998); Andersson et al. (1999); Nayyar and Owen (2006); Melatos (2007); van Eysden and Melatos (2008); Vigelius and Melatos (2009); Andersson et al. (2011). If a balance is assumed between the gravitational radiation-reaction torque and the accretion torque Papaloizou and Pringle (1978); Wagoner (1984); Bildsten (1998), then the strongest emitters of continuous gravitational waves are predicted to be sources in low-mass X-ray binaries (LMXBs), specifically those accreting at the highest rate Verbunt (1993); Watts et al. (2008).

Given the estimated ages ( yrs) and observed accretion rates of LMXBs (reaching near the Eddington limit of ), accretion is expected to spin-up the neutron star beyond the breakup frequency ( kHz for standard neutron star equations-of-state Cook et al. (1994); Ushomirsky et al. (2000b)). However, measured spin frequencies of LMXB neutron stars (from X-ray pulsations or thermonuclear bursts) so far range only from 95 to 619 Hz Chakrabarty et al. (2003); Bhattacharyya (2007); Watts et al. (2008); Galloway et al. (2010). The spin frequency cut-off lies well below breakup, and suggests the existence of a spin-down torque to balance the spin-up from accretion. A possible explanation, proposed by Papaloizou and Pringle (1978) and advanced by Wagoner (1984) and Bildsten (1998), is gravitational radiation. The torque-balance scenario implies a relation between the X-ray flux, spin frequency and gravitational wave strain; the more luminous the X-ray source, the greater the strain. Scorpius X-1 (Sco X-1), the brightest low-mass X-ray binary (LMXB), is therefore a promising target for periodic gravitational wave searches.

The global network of kilometer-scale, Michelson-type laser interferometers are sensitive to gravitational waves in the Hz frequency band. The Laser Interferometer Gravitational-Wave Observatory (LIGO) detectors achieved design sensitivity during the fifth science run (S5), between November 2005 and October 2007 Abbott et al. (2009a); Abadie et al. (2010a). LIGO consists of three Michelson interferometers (one with 4 km arms at Livingston, Louisiana, and two co-located at Hanford, Washington, with 4 km and 2 km arms) separated by km. Together the LIGO LIG (), Virgo VIR (); Accadia et al. (2012) and GEO600 Grote (2010); GEO () detectors form a world-wide network of broad-band interferometric gravitational wave observatories in an international effort to directly detect gravitational wave emission for the first time.

The LIGO Scientific Collaboration has so far published three main types of searches for periodic or continuous gravitational wave emission; targeted, directed and all-sky searches. Targeted searches are the most sensitive since they have the most tightly constrained parameter space. They target known sources, such as radio pulsars, with very well-constrained sky position, spin frequency and frequency evolution, and binary parameters (if any). These searches are fully-coherent, requiring accurately known prior phase information, making them computationally expensive to perform over large regions of parameter space Brady et al. (1998); Abbott et al. (2004, 2010). Targeted LIGO and Virgo searches have already set astrophysically interesting upper limits (e.g. beating theoretical indirect limits) on some pulsar parameters such as the gravitational wave strain from the Crab pulsar Abbott et al. (2008a); Abbott et al. (2010) and the Vela pulsar Abadie et al. (2011a). Directed searches aim at a particular sky location but search for unknown frequency (and/or frequency evolution). In most cases so far, directed searches have used a fully-coherent approach and approached the limits of computational feasibility. The search directed at the (possible) neutron star in the direction of the supernova remnant Cassiopeia A was able to beat indirect limits Abadie et al. (2010b). The third type of continuous gravitational wave search, the wide parameter-space searches, are also computationally intensive. They can involve searching over the entire sky or any comparably large parameter space, and usually employ semi-coherent approaches, combining short coherently analyzed segments in an incoherent manner. This process is tuned to balance the trade-off between reduced computational load and reduced sensitivity. The all-sky searches presented in Abbott et al. (2008b, 2009b, 2009c, 2009d); Abadie et al. (2012) target isolated neutron star sources (i.e. those not in binary systems). There is also an all-sky search for neutron stars in binary systems currently being run on LIGO’s latest S6 data run Goetz and Riles (2011).

Directed searches can also be made for known accreting neutron stars in binaries, and LIGO has previously conducted two of these searches for Sco X-1. The first, a coherent analysis using data from the second science run (S2), was computationally limited to the analysis of six-hour data segments from the LIGO interferometers, and placed upper limits on the wave strain of for two different 20 Hz bands Abbott et al. (2007a). This search utilized a maximum likelihood detection statistic based on matched filtering called the -statistic Jaranowski et al. (1998). The second, a directed version of the all-sky, stochastic, cross-correlation analysis, known as the “Radiometer” search, was first conducted on all 20 days of data from the S4 science run  Abbott et al. (2007b), and later on the year S5 data reporting a upper limit on the root-mean-square (RMS) strain (over the range 40 – 1500 Hz, with the minimum around 150 Hz)  Abadie et al. (2011b).

Semi-coherent search methods provide a compromise between sensitivity and the computational cost of a fully coherent search. They should be the most sensitive at fixed computing cost Brady and Creighton (2000); Prix and Shaltev (2012). A fast and robust search strategy for the detection of signals from binary systems in gravitational wave data was proposed in Messenger and Woan (2007). The signal from a source in a binary system is phase- (or frequency-) modulated due its periodic orbital motion, forming “sidebands” in the gravitational wave frequency spectrum. In searching detector data, this technique, called the sideband search, uses the same coherent (-statistic) stage as the previous coherent (S2) search. It then combines the frequency-modulated sidebands arising in -statistic data in a (computationally inexpensive) incoherent stage, reducing the need for a large template bank. This approach is based on a method that has successfully been employed in searches for binary pulsars in radio data Ransom et al. (2003).

Here we develop this sideband technique into a search pipeline and present a detailed description of how it is applied to gravitational wave detector data, as well as the expected sensitivity. The paper is set out as follows. Section II briefly describes the astrophysics of LMXBs and their predicted gravitational wave signature. The search algorithm is described in detail in Section III. Section IV outlines the parameter space of the search allowing primary sources to be identified in Section V. The statistical analysis of the results of the search is described in Section VI, along with a definition of the upper limits of the search. The sensitivity of the search is discussed in Section VII. A brief summary is provided in Section VIII, with a discussion of the limitations and future prospects of the search.

Ii Low Mass X-ray Binaries

LMXBs are stellar systems where a low-magnetic-field ( G) compact object (primary) accretes matter from a lower-mass (secondary) companion () Verbunt (1993); Tauris and van den Heuvel (2006). The compact objects in LMXB systems can be black holes, neutron stars or white dwarfs. For gravitational wave emission we are interested in LMXBs with neutron stars as the primary mass (typically ), since neutron stars can sustain the largest quadrupole moment.

Observations of thermal X-ray emission from the inner region of the accretion disk provide a measurement of the accretion rates in LMXBs. The range of observed accretions rates is broad, ranging from to the Eddington limit, Ushomirsky et al. (2000b). Some LMXBs also exhibit periodic pulsations or burst type behavior, and so provide a means of measuring the spin frequency of the neutron star in the system. The measured of these systems lie in the range of Hz Chakrabarty et al. (2003); Bhattacharyya (2007); Galloway et al. (2010). The broad range of accretion rates coupled with the estimated age of these systems ( years implied by evolutionary models van Paradijs and White (1995); Ushomirsky et al. (2000b)) would suggest a greater upper limit on observed spin frequency since accretion exerts substantial torque on the neutron star. However, none of these systems have yet been observed to spin at or near the breakup frequency kHz ( 1 kHz for most equations of state Cook et al. (1994); Ushomirsky et al. (2000b)). The maximum observed spin frequency falls far below the theorized breakup frequency and suggests a competing (damping) mechanism to the spin-up caused by accretion. One explanation for the observed spin frequency distribution of LMXBs is that the spin-up from the accretion torque is balanced by a gravitational wave spin-down torque Wagoner (1984); Bildsten (1998). Since the gravitational wave spin-down torque scales like (see Eq. 3 below), a wide range of accretion rates then leads to a rather narrow range of equilibrium rotation rates, as observed.

ii.1 Gravitational wave emission

Using the torque balancing argument from Wagoner (1984) and Bildsten (1998), we can estimate the gravitational wave strain amplitude emitted from accreting binary systems from their observable X-ray flux. This is a conservative upper limit as it assumes all angular momentum gained from accretion is completely converted into gravitational radiation.

The intrinsic strain amplitude for a system with angular spin frequency at a distance from an observer emitting gravitational waves via a mass quadrupole can be expressed as


where is the gravitational constant, is the speed of light and the quadrupole moment is a function of the ellipticity and moment of inertia Jaranowski et al. (1998). This can be expressed in terms of the spin-down (damping) torque due to gravitational radiation giving




The accretion torque applied to a neutron star of mass and radius accreting at a rate is given by


Assuming that the X-ray luminosity can be written as , the accretion rate can be expressed as a function of the X-ray flux , such that


since . In equilibrium, where gravitational radiation balances accretion torque, . The square of the gravitational wave strain from Eq. 6 can then be expressed in terms of the observable X-ray flux such that


Selecting fiducial values for the neutron star mass, radius, X-ray flux, and spin frequency (around the middle of the observed range), we can express the equilibrium strain upper limit in terms of and via

where . If the system emits gravitational waves via current quadrupole radiation instead, as is the case with -mode oscillations, the relation between gravitational wave frequency and spin frequency differs. In this case the preceding equations are modified slightly, requiring roughly Owen (2010). However these expressions, and the rest of the analysis except where otherwise noted, do not change if expressed in terms of the gravitational-wave frequency.

The resulting relation in Eq. LABEL:eq:hc_Fx_300 implies that LMXBs that accrete close to the Eddington limit are potentially strong gravitational wave emitters. Of these potentially strong sources, Sco X-1 is the most promising due to its observed X-ray flux  Watts et al. (2008).

Iii Search Algorithm

In this section we define our detection statistic and show how it exploits the characteristic frequency-modulation pattern inherent to sources in binary systems.

Fully-coherent, matched-filter searches for continuous gravitational waves can be described as a procedure that maximizes the likelihood function over a parameter space. The amplitude parameters (gravitational wave strain amplitude , inclination , polarization and reference phase ) can, in general, be analytically maximized, reducing the dimensions of this parameter space  Jaranowski et al. (1998). These parameters define the signal amplitudes in our signal model. Analytic maximization leaves the phase-evolution parameters (gravitational wave frequency and its derivatives and sky position ) to be numerically maximized over. Numerical maximization is accomplished through a scheme of repeated matched-filtering performed over a template bank of trial waveforms defined by specific locations in the phase parameter space, which is typically highly computationally expensive Brady et al. (1998); Abbott et al. (2007a); Watts et al. (2008).

The method we outline here makes use of the fact that we know the sky position of our potential sources and, hence, the phase evolution due to the motion of the detector can be accurately accounted for. We also know that the phase evolution due to the binary motion of the source will result in a specific distribution of signal power in the frequency domain. This distribution has the characteristics that signal power is divided amongst a finite set of frequency-modulation sidebands. The number of sidebands and their relative frequency spacing can be predicted with some knowledge of the binary orbital parameters.

In order to avoid the computational limits imposed by a fully-coherent parameter space search, we propose a single fully-coherent analysis stage, that accounts for detector motion only, is followed by a single incoherent stage in which the signal power contained within the frequency-modulated sidebands in summed to form a new detection statistic. This summing procedure is accomplished via the convolution of an approximate frequency domain power template with the output of the coherent stage.

The three main stages of the search, the -statistic, sideband template, and -statistic, are graphically illustrated in Figure 1. In this noise-free example, the frequency-modulation sidebands are clearly visible. The -statistic is also amplitude-modulated due to the daily variation of the detector antenna response, resulting in the amplitude-modulation applied to each frequency-modulation sideband. The second panel represents the approximate frequency domain template, a flat comb function with unit amplitude teeth (the spikes or delta functions). When convolved with the -statistic in the frequency domain we obtain the -statistic shown in the right-hand panel. The maximum power is clearly recovered at the simulated source frequency.

Figure 1: Graphical illustration of the sideband search pipeline, showing the frequency-modulated -statistic (left, red), the sideband template (middle, black), and their convolution, known as the -statistic (right, blue). In this noise free case, a signal of Hz with amplitude in a system with s, s, sky position rad, rad, and phase parameters rad, was simulated for a 10 day observation span. Frequency increases along the horizontal axis, which ranges from 199.998 to 200.002 Hz on each plot. In each case the location of the injected signal at 200 Hz is indicated by the vertical dashed black line.

The following sub-sections discuss each of the search components in more detail. Section III.1 presents the phase model used to characterize the gravitational wave signal from a source in a binary system. The signal model is introduced in Sec. III.2. The -statistic is introduced in Section III.3 and its behavior as a function of search frequency is described in Sec. III.4. Section III.5 then goes on to describe how matching a filter for an isolated neutron star system to a signal from a source in a binary system results in frequency-modulated sidebands appearing in the output of the -statistic. The detection statistic for gravitational wave sources in binary orbits is fully described as the incoherent sum of frequency-modulated -statistics in Sec. III.6. The simplest, unit amplitude, sideband template, and its justification over a more realistic template, are discussed in Sec. III.7. A more sensitive implementation, incorporating an approximate binary phase model in the calculation of the -statistic and reducing the width of the frequency-modulated sideband pattern by the fractional errors on the semi-major axis parameters, is discussed in Sec. III.8. Beginning in this section and continuing in the following sections, we have made a distinction between the intrinsic values of a search parameter (denoted with a subscript zero) and the observed values (no subscripts).

iii.1 Phase model

The phase of the signal at the source can be modeled by a Taylor series such that


where represents the time derivative of the gravitational wave frequency evaluated at reference time . For the purposes of this work we restrict ourselves to a monochromatic signal and hence set for all and define as the intrinsic frequency. We discuss this choice in Sec. IV.2. The phase received at the detector is , where we define the retarded times measured at the Solar system barycenter (SSB) and detector as and respectively. The relation between and is a function of source sky position relative to detector location and, since we only concern ourselves with sources of known sky position, we assume that the effects of phase contributions from detector motion can be exactly accounted for. For this reason we work directly within the SSB frame. The relationship between the source and retarded times for a non-relativistic eccentric binary orbit is given by Taylor and Weisberg (1989)


where is the light crossing time of the semi-major axis projected along the line of sight. The orbital eccentricity is defined by and the argument of periapse, given by , is the angle in the orbital plane from the ascending node to the direction of periapse. The variable is the eccentric anomaly defined by , where is the time of passage through periapse (the point in the orbit where the two bodies are at their closest) and is the orbital period.

It is expected that dissipative processes within LMXBs drive the orbits to near circularity. In the low eccentricity limit , we obtain the following expression


where and we have used the time of passage through the ascending node as our reference phase in the first term. For circular orbits, where , the expression reduces to only this first term. Using is sensible in this case since and , are not defined in a circular orbit. Note that the additional eccentric terms are periodic at multiples of twice the orbital frequency. Using Eq. 10, we would expect to accumulate timing errors of order s for the most eccentric known LMXB systems. We shall return to this feature in Sec. IV.5.

To write the gravitational wave phase as a function of SSB time, we invert Eq. 10 to obtain . The binary phase can be corrected for in a standard matched filter approach, where Eq. 9 is solved numerically. In our method we instead approximate the inversion as


for circular orbits.1 Under this approximation the signal phase can be represented as a linear combination of the phase contributions from the spin of the neutron star and from the binary orbital motion of the source , such that


iii.2 Signal model

We model the data collected by a detector located at the SSB as the signal plus stationary Gaussian noise so that




where we employ the Einstein summation convention for . The coefficients are independent of time, detector location and orientation. They depend only on the signal amplitude parameters , where is the dimensionless gravitational wave strain amplitude, is the gravitational wave polarization angle, is the source inclination angle and is the signal phase at a fiducial reference time. The coefficients are defined as




are the polarization amplitudes. The time dependent signal components are defined as


where is the signal phase at the detector (which we model as located at the SSB) given by Eq.11 and the antenna pattern functions and are described by Eqs. (12) and (13) in Jaranowski et al. (1998).

iii.3 -statistic

The -statistic is a matched-filter based detection statistic derived via analytic maximization of the likelihood over unknown amplitude parameters Jaranowski et al. (1998). Let us first introduce the multi-detector inner product


where indexes each detector and is the detector single-sided noise spectral density. We modify the definitions of Jaranowski et al. (1998) and Prix (2007) to explicitly include gaps in the time-series by introducing the function which has value when data is present and otherwise. This also allows us to extend the limits of our time integration to since the window function will naturally account for the volume and span of data for each detector.

The -statistic itself is defined as


where form the matrix inverse of and we follow the shorthand notation of Prix (2007) defining and . Evaluation of leads to a matrix of the form


where the components


are antenna pattern integrals. For a waveform with exactly known phase evolution in Gaussian noise, the -statistic is a random variable distributed according to a non-central -distribution with 4 degrees of freedom. The non-centrality parameter is equal to the optimal signal-to-noise ratio (SNR)


such that the expectation value and variance of are given by


respectively. In the case where no signal is present in the data, the distribution becomes a central -distribution with 4 degrees of freedom.

iii.4 -statistic and mismatched frequency

In this section we describe the behavior of the -statistic as a function of search frequency for a fixed source frequency . In this case the inner product that defines becomes


where are the components of a signal with frequency , and is a function of search frequency . If we focus on the component as an example we find that

where we note that the product of cosine functions results in an integrand that contains frequencies at and . Since both and are functions that evolve on timescales of hours–days we approximate the contribution from the component as averaging to zero. We are left with


where we have defined the result of the complex integral as . This is the Fourier transform of the antenna pattern functions weighted by the window function and evaluated at . It is equal to the quantity when its argument is zero.

The quantity (and similarly and ) are periodic quantities with periods of 12 and 24 hours plus a non-oscillating component. When in a product with a sinusoidal function and integrated over time they will result in discrete amplitude-modulated sidebands with frequencies at Hz where represents the orbital period of the Earth (1 sidereal day). We will ignore all but the zero-frequency components of these functions for the remainder of this paper. We do note that complications regarding the overlap of amplitude-modulated and frequency-modulated sidebands (discussed in the next section) will only arise for sources in binary orbits with periods equal to those present in the antenna pattern functions.

In addition, the window function describing the gaps in the data will influence . For a gap-free observation the window function serves to localize signal power to within a frequency range where is the typical observation length. When gaps are present this range is broadened and has a deterministic shape given by the squared modulus of the Fourier transform of the window function. We can therefore use a further approximation that


where the Fourier transform of the window function is normalized by such that it has a value of unity at the true signal frequency.

We now define the antenna-pattern weighted window function as


which is true for observation times days. This complex window function has the property that is a real quantity with maximum absolute value of unity when the template frequency matches the true signal frequency.

Finally we are able to combine Eqs. 2627, and 28 to obtain


which together with similar calculations for the additional components in Eq. 24 give us

as the complete set of inner products between frequency-mismatched signal components where I is the identity matrix. Note that when this expression reduces to Eq. 19.

If we now form the expectation value of the -statistic (Eq. 19) for mismatched frequencies we find that


Here we see that the fraction of the optimal SNR that contributes to the non-centrality parameter of the -statistic distribution is reduced by evaluation of the mod-squared of the antenna-pattern weighted window function with a non-zero argument.

iii.5 Frequency-modulation and the -statistic

We now consider the computation of the -statistic in the case where the data contains a signal from a source in a circular binary orbit but the phase model used in the -statistic template is that of a monochromatic signal of frequency . We again expand as done in Eq. 24 where no prime indicates the signal and the prime represents the monochromatic template. We again focus on the mismatched signal inner-product as an example. Starting with Eq. III.4 we discard the rapidly oscillating terms inside the integral that will average to zero. We are then left with


where the final term involving the exponential of a sinusoidal function can be represented using the Jacobi-Anger expansion


where is the order Bessel function of the first kind. This expansion allows us to transform the binary phase term into an infinite sum of harmonics such that we can now write


It follows that all of the signal components can be expanded in the same way giving us


where we have truncated the infinite summation (explained below) and defined the monochromatic modulated sideband frequencies and their respective phases as


The Jacobi-Anger expansion has allowed us to represent the complex phase of a frequency-modulated signal as an infinite sum of discrete signal harmonics, or sidebands, each separated in frequency by Hz. Each is weighted by the Bessel function of order where indexes the harmonics and has a complex phase factor determined by the orbital reference time . In the limit where the order exceeds the argument, , the Bessel function rapidly approaches zero allowing approximation of the infinite sum in Eqs. 33 and 34 as a finite sum over the finite range where . The summation format of Eq. III.5 highlights the effects of the binary phase modulation. The signal can be represented as the sum of discrete harmonics at frequencies centered on the intrinsic frequency , where each harmonic peak is separated from the next by Hz.

Combining Eqs. 1924, and III.5, we can express the expectation value of the -statistic for a binary signal as a function of search frequency as


This expression should be interpreted in the following manner. For a given search frequency the contribution to the non-centrality parameter (the SNR dependent term) is equal to the sum of all sideband contributions at that frequency. Each sideband will contribute a fraction of the total optimal SNR weighted by the order Bessel function squared, but will also be strongly weighted by the window function. The window function will only contribute significantly if the search frequency is close to the sideband frequency. Hence, at a given search frequency close to a sideband, for observation times , the sidebands will be far enough separated in frequency such that only one sideband will contribute to the -statistic.

iii.6 -statistic

The -statistic is numerically maximized over the phase parameters of the signal on a discrete grids. For this search the search frequency is such a parameter and consequently the -statistic is computed over a uniformly spaced set of frequency values spanning the region of interest. In this section we describe how this -statistic frequency-series can be used to approximate a search template that is then used to generate a new statistic sensitive to signals from sources in binary systems.

The expectation value of the -statistic (Eq. 37) resolves into localized spikes at frequencies separated by Hz and centered on the intrinsic gravitational wave frequency . A template based on this pattern with amplitude defined by , takes the general form


with and where we make a distinction between the intrinsic (unknown) values of each parameter (subscript zero) and values selected in the template construction (denoted with a prime). The window function is dependent only on the times for which data is present and is, therefore, also known exactly.

We define our new detection statistic as


where the sum over the index indicates the sum over the discrete frequency bins and is the search frequency. Since the template’s “zero frequency” represents the intrinsic gravitational wave frequency, corresponds to the intrinsic frequency. We see that the -statistic is, in fact, the convolution of the -statistic with our template, assuming the template is constant with search frequency (an issue we address in the next section).

The benefit of this approach is that the computation of the -statistic for a known sky position and without accounting for binary effects has relatively low computational cost. Similarly, the construction of a template on the -statistic is independent of the orbital phase parameter and only weakly dependent upon the orbital semi-major axis and eccentricity. The template is highly dependent upon the orbital period, which, for the sources of interest, is known to high precision. Also, since the -statistic is the result of a convolution, we can make use of the convolution theorem and the speed of the Fast-Fourier-Transform (FFT). Computing for all frequencies requires only three applications of the FFT. In practice, the -statistic is computed using


which is simply the inverse Fourier transform of the product of the Fourier transforms of the -statistic and the template. The -statistic and the template are both sampled on the same uniform frequency grid containing frequency bins. The -statistic is then also output as a function of the same frequency grid.

iii.7 Choice of -statistic template

Treating the -statistic as the pre-processed input dataset to the -statistic computation, it might be assumed that the optimal choice of template is that which exactly matches the expected form of in the presence of a signal. As shown in Fig. 2, this approach is highly sensitive to the accuracy with which the projected orbital semi-major axis is known.

We instead propose the use of a far simpler template: one that captures the majority of the information contained within the -statistic and, by design, is relatively insensitive to the orbital semi-major axis. We explicitly choose


for discrete frequency , where is the Kronecker delta-function. The frequency index , defined by,


is a function of the best-guess orbital period and the frequency resolution . The function returns the integer closest to its argument. The template is therefore composed of the sum of unit amplitude “spikes” positioned at discrete frequency bins closest to the predicted locations of the frequency-modulated sidebands (relative to the intrinsic gravitational wave frequency). The subscript F refers to the constant amplitude, “Flat” template.

If we now convolve this template with the frequency-modulated -statistic we obtain the corresponding -statistic, which reduces to


Equation 43 is simply the sum of values taken from discrete frequency bins positioned at the predicted locations of the frequency-modulation sidebands. An example is shown in the right hand panel of Fig. 1, where the most significant -statistic value is located at , and is the point where all sidebands included in the sum contain some signal.

From Eq. 37, the expectation value for the -statistic using the flat-template can be expressed as


where the argument of the window function is the frequency difference between the location of the signal sideband and the template sideband on the discrete frequency grid. Note that the -statistic is a sum of statistically independent non-central statistics and hence the result is itself a non-central statistic, i.e. with degrees of freedom, where is the number of sidebands in the expected modulation pattern. The non-centrality parameter is equal to the sum of the non-centrality parameters from each of the summed values. For a flat template with perfectly matched intrinsic frequency and orbital period , infinite precision , and where the number of sidebands included in the analysis matches or exceeds the true number, the second term in the above equation reduces to . In this case we will have recovered all of the power from the signal but also significantly increased the contribution from the noise through the incoherent summation of -statistic from independent frequency bins. In general, where the orbital period is known well, but not exactly, and the frequency resolution is finite, the signal power recovery will be reduced by imperfect sampling of the window function term in Eq. 44, i.e. evaluation at arguments .

In terms of the generic template defined in equation Eq. 38, the discrete-frequency flat template is approximately equivalent to the weighting scheme . A more sensitive approach could use


for the template, following the expected form of the -statistic given in Eq. 37 and using a subscript B to denote Bessel function weighting. Although this would increase sensitivity for closely matched signal templates (constructed with well constrained signal parameters), this performance is highly sensitive to the number of sidebands included in the template and therefore sensitive to the semi-major axis since . This is mainly due to the “double horned” shape of the expected signal (see the left hand panel of Fig. 1). A large enough offset between the true and assumed semi-major axis will significantly change the template’s overlap with the sidebands in the -statistic and reduce the significance of the -statistic. Considering the semi-major axis is not well constrained for many LMXBs, a search over many templates would be necessary, each with incrementally different semi-major axis values.

The simpler, flat-template (Eq. 41) has the benefit of being far more robust against the semi-major axis uncertainty. In this case the semi-major axis parameter controls only the number of sidebands to use in the template and does not control the weighting applied to each sideband. It also simplifies the statistical properties of the -statistic, making a Bayesian analysis of the output statistics (as described in Section VI) far easier to apply.

The receiver operator characteristic (ROC) curves shown in Fig. 2 compare the performance of the sideband search with both choices of template (flat and Bessel function weighted) for the case of signal with optimal SNR . As seen from the figure, the Bessel function weighted template for exact number of sidebands provides improved sensitivity over the flat template. However, when considering the possible (and highly likely) error on the number of sidebands in the template, the performance of the Bessel template is already drastically diminished, even with only a error on the semi-major axis parameter. It is also interesting to note that for the flat-template the result of an incorrect semi-major axis is asymmetric with respect to an under or over-estimate. The sensitivity degradation is far less pronounced when the template has over-estimated the semi-major axis and, therefore, also over-estimated the number of sidebands. This feature is discussed in more detail in Sec. IV.4.

Figure 2: ROC curves comparing performance of the flat (blue) and Bessel function weighted (red) templates, described by Eqs 41 and 45 respectively. The theoretical (fine black) curve is constructed from a non-central distributed statistic with non-centrality parameter and represents the expected performance of the flat template. Dashed and dotted curves represent a template with a positive and negative error on the semi-major axis respectively. Here the signal parameters were chosen such that the number of sidebands were and curves were constructed using realizations of noise.

iii.8 Approximate binary demodulation

When a putative source has a highly localized position in the sky, the effect of the Earth’s motion with respect to the SSB can be accurately removed from the signal during the calculation of the -statistic. This leaves only the Doppler modulation from the binary orbit. It is also possible to demodulate the binary orbit (Doppler) modulation in the -statistic calculation provided the binary orbital parameters are well known. A fully-coherent (sky position- and binary-) demodulated -statistic search would be very sensitive to any errors in the sky position of binary orbital parameters. It would therefore be necessary to construct a bank of templates spanning the parameter space defined by the uncertainties in these parameters. Adding dimensions to the parameter space increases computational costs and the search becomes unfeasible considering we are already searching over frequency. A fully coherent search of this type would be possible for known sources with known emission frequencies, for example pulsing sources like millisecond pulsars.

In this section we show how prior information regarding the binary orbit of a source can be used to increase the sensitivity of our semi-coherent approach, without increasing computational costs. By performing a “best guess” binary phase demodulation within the -statistic, we show that the number of sidebands in the template is reduced by a factor proportional to the fractional uncertainty in the orbital semi-major axis. Consequently a reduction in the number of sidebands increases the sensitivity of the search by reducing the number of degrees of freedom (see Sec. VI).

Expressing our current best estimate for each parameter as the sum of the true value and an error , such that


we can determine the phase offset of the binary orbit from the error in the binary orbital parameters. The offset in phase is the difference between the true and best estimate binary phase and using Eq. 12b can be approximated by




Here we have expanded the binary phase difference to leading order in the parameter uncertainties and obtained a phase expression similar in form to the original binary phase. In the specific regime where the fractional uncertainty in the orbital semi-major axis far exceeds the fractional uncertainty in the intrinsic frequency we see that the first term in Eq. 48 becomes . Similarly, if the fractional uncertainty in the orbital semi-major axis also far exceeds the deviation in orbital angular position then the second term . This is generally the case for the known LMXBs (see Table 1) and in this regime can be accurately approximated as unity, yielding


Hence, after approximate binary demodulation, the argument of the Bessel function and the summation limits in the expected form of the -statistic (in Eq. 37 for example) can be replaced with . The number of frequency-modulated sidebands is now reduced by a factor of . We must stress that is an unknown quantity and is the difference between the best estimate value of and the true value . The -statistic after such a demodulation process will therefore have a reduced but unknown number of sidebands, although it will still retain the standard sideband frequency spacing . The sideband phasing will also be unknown due to the presence of the phase term but is of no consequence to the search since the -statistic is insensitive to phase.

Iv Parameter space

In this section we will discuss each of the parameters involved in the search and how the search sensitivity depends upon the uncertainty in these parameters. Demodulation of the signal phase due to the Earth’s motion requires accurate knowledge of the source sky position. If the observation time is long enough, we need to consider the sky position as a search parameter, as discussed in Sec. IV.1. The gravitational wave frequency is the primary search parameter. In Sec. IV.2 we discuss the limitations on our search strategy due to its uncertainty. The orbital period and semi-major axis are discussed in Sections IV.3 and IV.4 respectively. The effects of ignoring the orbital eccentricity are discussed in Sec. IV.5.

iv.1 Sky position and proper motion

In order to quantify the allowable uncertainty in sky position we will define a simplistic model describing the phase received at Earth from a monochromatic source at infinity at sky position . If we neglect the detector motion due to the spin of the Earth and consider only the Earth’s orbital motion then we have


where is the signal frequency, and are the true right ascension and declination and and are the distance of the Earth from the SSB and the Earth’s orbital angular frequency respectively. We also define a reference time that represents the time at which the detector passes through the vernal equinox. For an observed sky position the corresponding phase offset amounts to


where we have expanded the expression to leading order in the sky position errors. We now make the reasonable assumption that our analysis would be unable to tolerate a deviation in phase between the signal and our template of more than radian over the course of an observation on the same timescale of the Earth’s orbit.

If we also notice that the worst case scenario (smallest allowable sky position errors) corresponds to sky positions for which the trigonometric terms in the previous expression are largest, i.e. of order unity, then we have


If we consider signals of frequency kHz, this gives a maximum allowable sky position offset of mas. This expression also validates our model assumption that the sky position sensitivity to the Earth spin would be dominated by the effect from the Earth orbit for long observation times.

A similar argument can be made for the proper motion of the source where we would be safe to model the sky position as fixed if the change , over the course of the observation also satisfied Eq. 53.

iv.2 Spin frequency

The spin frequency of some LMXBs can be directly measured from X-ray pulsations, believed to originate from a hot-spot on the stellar surface, where accreted material is funneled onto the magnetic pole with the magnetic axis generally misaligned with the spin axis. X-ray pulsations have been observed in 13 LMXB systems so far, three of which are intermittent Lo et al. (2011).

Some LMXBs exhibit recurrent thermonuclear X-ray bursts. Fourier spectra reveal oscillations during the rise and tail of many bursts, which are believed to originate from asymmetric brightness patterns on the stellar surface. In seven LMXBs which exhibit both pulsations and bursts, the asymptotic burst oscillation frequency at late times matches the pulse frequency. Where there are no pulsations, many bursts need to be observed to measure the asymptotic burst oscillation reliably. The spin frequency of an additional ten systems has been determined from burst oscillations only Watts (2012), but due to the uncertainties involved, are usually quoted to within uncertainties of Hz.

Another class of LMXBs exhibit high frequency quasi-periodic oscillationss in their persistent X-ray emission. These kHz QPOs usually come in pairs, although singles and triples are occasionally observed and the QPO peak frequencies usually change over time. In some cases the separation of the QPO peaks is roughly constant, but this is not always the case van der Klis (1998); Zhang et al. (2012); van der Klis et al. (1996). For the few QPO systems where can be determined from pulses or burst oscillations there has been no evidence suggesting consistency with an existing QPO model that links the QPO and spin frequencies. For our purposes, is considered unknown in sources without pulsations or confirmed bursts.

In addition to potentially broad uncertainties in , we know little about how its value may fluctuate over time due to accretion. Changes in the accretion flow will exert a time varying torque on the star which will result in a stochastic wandering of the spin frequency. In this case the signal can no longer be assumed monochromatic over a given observation time. To quantify the resulting phase wandering, we assume that the fluctuating component of the torque flips sign randomly on the timescale consistent with the inferred variation in accretion rate. If the mean torque due to steady-state disk-fed accretion, then the angular spin frequency experiences a random walk with step size , where is the stellar moment of inertia. After time , the root-mean-square drift is


This frequency drift will wander outside a Fourier frequency bin width if . If we choose such that the accretion rate can vary up to a factor of two in this time, then the worst case leads to the restriction


This is the primary reason why an application of the the basic sideband search, as described here, must be limited in the length of data it is allowed to analyse. By exceeding this limit it becomes increasingly likely that the spin wandering inherent to a true signal will cause signal power to leak between adjacent frequency bins. Consequently the assumption that -statistic signal power is localized in frequency-modulated sidebands will become invalid and the sensitivity of the -statistic will deteriorate.

iv.3 Orbital Period

The sideband search relies on relatively precise electromagnetic (EM) measurements of the orbital period in order to construct a search template. The duration of the orbit defines the minimum observation time, since is required before sidebands become clearly resolved in the spectrum Ransom et al. (2003). The uncertainty in the orbital period will determine the number of templates required to fully sample the search space, or equivalently, the maximum observation time allowed for a single value of .

We will now provide an indication of the sensitivity of the search to errors in the orbital period. If our estimate (observation) of the orbital period is offset from the true value by an amount , we would expect the error to seriously affect the -statistic recovered from the search once it is large enough to shift the outermost “tooth” in the sideband template by one canonical frequency bin away from the true sideband location. In this case, the offset between the template and true sideband frequency is proportional to the number of sidebands from the central spike. There will be low mismatch at the center of the template extending to mismatch at the edges. It follows that the average signal power recovered from such a mismatched template will be and therefore serves as a useful threshold by which to determine the maximum allowed .

If we use the measured value of as our template parameter, the template centered at frequency then consists of unit spikes (or teeth) separated by . Assuming that the central spike is exactly equal to the true intrinsic gravitational wave frequency, any errors in the orbital period will be propagated along the comb, causing the offset between the true and template frequency of any particular sideband to grow progressively larger. The frequency difference between the outermost template sideband, at frequency , and the outermost signal sideband at , is given by


for . To satisfy the condition described above we now require that this frequency shift should be less than the size of one frequency bin. The true frequency bin size is determined by the observation time span and is given by


where is the resolution used in the -statistic calculation.2 Using Eqs. 56 and 57 and imposing the condition that provides an estimate for the maximum allowable (orbital period limited) coherent observation timespan,


Given a relatively poorly constrained orbital period uncertainty, this restriction may provide too short a duration. This could be because it is then in conflict with the requirement that or simply because more SNR is desired from the signal. In either case, the orbital period space must then be divided into templates with spacing derived from simply rearranging Eq. 58 to solve for . In relative terms the sideband search places very strong constraints on the prior knowledge of the orbital period compared to the other search parameters.

iv.4 Semi-major axis

An error in the value of the orbital semi-major axis results in an incorrect choice for the number of sidebands in the template. As can be seen in Eq. 44, an underestimate results in the summation of a fraction of the total power in the signal whereas an overestimate results in a dilution of the total power by summing additional noise from sideband frequencies containing no signal contribution.

If we define the true semi-major axis parameter as the measured value and some (unknown) fraction (where ) of the measurement error (i.e. ), we can investigate the effects of errors on the semi-major axis parameter in terms of this offset parameter . We consider the advantage of using a deliberately offset value instead of the observed value in order to minimize losses in recovered SNR.

The ROC curves shown in Fig. 3 show the effects of these offsets, and clearly illustrate degradation in the performance of the -statistic as () increases. The reduction in detection confidence at a given false alarm probability is much faster for (), when the template underestimates the width of the sideband structure, than for (). This is natural considering the “horned” shape of the signal (see the left hand panel of Fig. 1 and Sec. III.7). Although it is already clear from this figure that the performance of the search is not symmetric about , this asymmetry is much better illustrated in Fig. 4 where for different values of the false alarm rate we show the detection probability plotted against the offset parameter .

This plot provides us with a rough scheme by which to improve the search performance by exploiting the asymmetry in search sensitivity with respect to . In general, we are keen to probe the low false alarm and high detection probability regime in which it is clear that using a template based on an orbital semi-major axis value the best estimate reduces