Implementation of the frequency-modulated sideband search method for gravitational waves from low mass X-ray binaries
Letizia Sammut, Christopher Messenger, Andrew Melatos, Benjamin J. Owen
IImplementation of the frequency-modulated sideband search method for gravitational waves fromlow mass X-ray binaries
L. Sammut,
1, 2, ∗ C. Messenger,
3, 2
A. Melatos, and B.J. Owen University Of Melbourne, Victoria 3052, Australia Albert Einstein Institute, Callinstraße 38,D-30167 Hannover,Germany University of Glasgow, Glasgow, G12 8QQ, United Kingdom The Pennsylvania State University, University Park, Pennsylvania 16802, USA (Dated: August 21, 2018)(commitID: 1d41f0b94a6df042af7e99a97529058a81632422)(commitDate: Wed Nov 6 12:08:06 2013 + F -statistic. The sideband search is based onthe incoherent summation of these frequency-modulated F -statistic sidebands. It provides a new detectionstatistic for sources in binary systems, called the C -statistic. The search is well suited to low-mass X-raybinaries, the brightest of which, called Sco X-1, is an ideal target candidate. For sources like Sco X-1, withwell constrained orbital parameters, a slight variation on the search is possible. The extra orbital informationcan be used to approximately demodulate the data from the binary orbital motion in the coherent stage, beforeincoherently summing the now reduced number of sidebands. We investigate this approach and show that itimproves the sensitivity of the standard Sco X-1 directed sideband search. Prior information on the neutronstar inclination and gravitational wave polarization can also be used to improve upper limit sensitivity. Weestimate the sensitivity of a Sco X-1 directed sideband search on 10 days of LIGO data and show that it canbeat previous upper limits in current LIGO data, with a possibility of constraining theoretical upper limits usingfuture advanced instruments. PACS numbers: 95.85.Sz, 97.60.Jd, 97.80.Jp
I. INTRODUCTION
Rotating neutron stars are the main candidates for sourcesof persistent, periodic gravitational radiation detectable byground-based, long-baseline gravitational wave interferom-eters. Time varying quadrupole moments in (and thusgravitational-wave emission from) these sources can resultfrom deformations of the solid crust (and possibly a solidcore) supported by elastic stresses [1–6], deformations of var-ious parts of the star supported by magnetic stresses [7–12],or free precession [13, 14] or long-lived oscillation modes ofthe entire star [15–19].Neutron stars in accreting binary systems are an importantsub-class of periodic gravitational wave sources. Accretionmay trigger or enhance the aforementioned gravitational-waveemission mechanisms, creating or driving the quadrupole mo-ment toward its maximum value through thermal, magnetic, orother e ff ects [1, 16, 20–24]. If a balance is assumed betweenthe gravitational radiation-reaction torque and the accretiontorque [1, 25, 26], then the strongest emitters of continuousgravitational waves are predicted to be sources in low-massX-ray binaries (LMXBs), specifically those accreting at thehighest rate [27, 28].Given the estimated ages ( ∼ yrs) and observed accre-tion rates of LMXBs (reaching near the Eddington limit of ∗ Electronic address: [email protected] ˙ M Edd = × − M (cid:12) yr − ), accretion is expected to spin-up theneutron star beyond the breakup frequency ( ∼ . ff lies well be-low breakup, and suggests the existence of a spin-down torqueto balance the spin-up from accretion. A possible explanation,proposed by Papaloizou and Pringle [25] and advanced byWagoner [26] and Bildsten [1], is gravitational radiation. Thetorque-balance scenario implies a relation between the X-rayflux, spin frequency and gravitational wave strain; the moreluminous the X-ray source, the greater the strain. ScorpiusX-1 (Sco X-1), the brightest LMXB, is therefore a promisingtarget for periodic gravitational wave searches.The global network of kilometer-scale, Michelson-typelaser interferometers are sensitive to gravitational waves in the O (10 − ∼ a r X i v : . [ g r- q c ] N ov three main types of searches for periodic or continuousgravitational wave emission; targeted, directed and all-skysearches. Targeted searches are the most sensitive since theyhave the most tightly constrained parameter space. They tar-get known sources, such as radio pulsars, with very well-constrained sky position, spin frequency and frequency evo-lution, and binary parameters (if any). These searches arefully-coherent, requiring accurately known prior phase infor-mation, making them computationally expensive to performover large regions of parameter space [41–43]. Targeted LIGOand Virgo searches have already set astrophysically interestingupper limits (e.g. beating theoretical indirect limits) on somepulsar parameters such as the gravitational wave strain fromthe Crab pulsar [43, 44] and the Vela pulsar [45]. Directedsearches aim at a particular sky location but search for un-known frequency (and / or frequency evolution). In most casesso far, directed searches have used a fully-coherent approachand approached the limits of computational feasibility. Thesearch directed at the (possible) neutron star in the direction ofthe supernova remnant Cassiopeia A was able to beat indirectlimits [46]. The third type of continuous gravitational wavesearch, the wide parameter-space searches, are also computa-tionally intensive. They can involve searching over the entiresky or any comparably large parameter space, and usually em-ploy semi-coherent approaches, combining short coherentlyanalyzed segments in an incoherent manner. This process istuned to balance the trade-o ff between reduced computationalload and reduced sensitivity. The all-sky searches presentedin [47–51] target isolated neutron star sources (i.e. those notin binary systems). There is also an all-sky search for neutronstars in binary systems currently being run on LIGO’s latestS6 data run [52].Directed searches can also be made for known accretingneutron stars in binaries, and LIGO has previously conductedtwo of these searches for Sco X-1. The first, a coherent analy-sis using data from the second science run (S2), was computa-tionally limited to the analysis of six-hour data segments fromthe LIGO interferometers, and placed 95% upper limits on thewave strain of h ≈ × − for two di ff erent 20 Hz bands[53]. This search utilized a maximum likelihood detectionstatistic based on matched filtering called the F -statistic [54].The second, a directed version of the all-sky, stochastic, cross-correlation analysis, known as the “Radiometer” search, wasfirst conducted on all 20 days of data from the S4 science run[55], and later on the ∼ h > × − (over the range 40 – 1500 Hz, with the minimum around 150Hz) [56].Semi-coherent search methods provide a compromise be-tween sensitivity and the computational cost of a fully coher-ent search. They should be the most sensitive at fixed com-puting cost [57, 58]. A fast and robust search strategy forthe detection of signals from binary systems in gravitationalwave data was proposed in [59]. The signal from a source ina binary system is phase- (or frequency-) modulated due itsperiodic orbital motion, forming “sidebands” in the gravita-tional wave frequency spectrum. In searching detector data,this technique, called the sideband search, uses the same co- herent ( F -statistic) stage as the previous coherent (S2) search.It then combines the frequency-modulated sidebands arisingin F -statistic data in a (computationally inexpensive) incoher-ent stage, reducing the need for a large template bank. Thisapproach is based on a method that has successfully been em-ployed in searches for binary pulsars in radio data [60].Here we develop this sideband technique into a searchpipeline and present a detailed description of how it is ap-plied to gravitational wave detector data, as well as the ex-pected sensitivity. The paper is set out as follows. SectionII briefly describes the astrophysics of LMXBs and their pre-dicted gravitational wave signature. The search algorithm isdescribed in detail in Section III. Section IV outlines the pa-rameter space of the search allowing primary sources to beidentified in Section V. The statistical analysis of the resultsof the search is described in Section VI, along with a defini-tion of the upper limits of the search. The sensitivity of thesearch is discussed in Section VII. A brief summary is pro-vided in Section VIII, with a discussion of the limitations andfuture prospects of the search. II. LOW MASS X-RAY BINARIES
LMXBs are stellar systems where a low-magnetic-field( (cid:46) G) compact object (primary) accretes matter from alower-mass (secondary) companion ( < M (cid:12) ) [27, 61]. Thecompact objects in LMXB systems can be black holes, neu-tron stars or white dwarfs. For gravitational wave emissionwe are interested in LMXBs with neutron stars as the primarymass (typically ∼ . M (cid:12) ), since neutron stars can sustain thelargest quadrupole moment.Observations of thermal X-ray emission from the inner re-gion of the accretion disk provide a measurement of the accre-tion rates in LMXBs. The range of observed accretions ratesis broad, ranging from 10 − M (cid:12) yr − to the Eddington limit,˙ M Edd = × − M (cid:12) yr − [30]. Some LMXBs also exhibitperiodic pulsations or burst type behavior, and so provide ameans of measuring the spin frequency ν s of the neutron starin the system. The measured ν s of these systems lie in therange of 95 ≤ ν s ≤
619 Hz [31–33]. The broad range of ac-cretion rates coupled with the estimated age of these systems( ∼ years implied by evolutionary models [30, 62]) wouldsuggest a greater upper limit on observed spin frequency sinceaccretion exerts substantial torque on the neutron star. How-ever, none of these systems have yet been observed to spinat or near the breakup frequency v b ∼ . v b (cid:38) ν (see Eq. 3 below), a wide range of accre-tion rates then leads to a rather narrow range of equilibriumrotation rates, as observed. A. Gravitational wave emission
Using the torque balancing argument from Wagoner [26]and Bildsten [1], we can estimate the gravitational wave strainamplitude emitted from accreting binary systems from theirobservable X-ray flux. This is a conservative upper limit as itassumes all angular momentum gained from accretion is com-pletely converted into gravitational radiation.The intrinsic strain amplitude h for a system with angularspin frequency Ω s = πν s at a distance d from an observeremitting gravitational waves via a mass quadrupole can be ex-pressed as h = GQc d Ω , (1)where G is the gravitational constant, c is the speed of lightand the quadrupole moment Q = (cid:15) I is a function of the ellip-ticity (cid:15) and moment of inertia I [54]. This can be expressed interms of the spin-down (damping) torque N GW due to gravita-tional radiation giving h = G c d Ω N GW , (2)where N GW = GQ c Ω . (3)The accretion torque N a applied to a neutron star of mass M and radius R accreting at a rate ˙ M is given by N a = ˙ M √ GMR . (4)Assuming that the X-ray luminosity can be written as L X = GM ˙ M / R , the accretion rate ˙ M can be expressed as a functionof the X-ray flux F X , such that˙ M = π Rd GM F X , (5)since L X = π d F X . In equilibrium, where gravitational radi-ation balances accretion torque, N GW = N a . The square of thegravitational wave strain from Eq. 6 can then be expressed interms of the observable X-ray flux such that h = F X ν s (cid:32) GR M (cid:33) / . (6)Selecting fiducial values for the neutron star mass, radius, X-ray flux, and spin frequency (around the middle of the ob-served range), we can express the equilibrium strain upperlimit h EQ0 in terms of ν s and F X via h EQ0 = . × − (cid:32) F x F ∗ (cid:33) / (cid:18) R (cid:19) / (cid:32) . M (cid:12) M (cid:33) / (cid:32) ν s (cid:33) / , (7)where F ∗ = − erg cm − s − . If the system emits gravita-tional waves via current quadrupole radiation instead, as is the case with r -mode oscillations, the relation between grav-itational wave frequency and spin frequency di ff ers. In thiscase the preceding equations are modified slightly, requiringroughly ν s → (2 / ν s [63]. However these expressions, andthe rest of the analysis except where otherwise noted, do notchange if expressed in terms of the gravitational-wave fre-quency.The resulting relation in Eq. 7 implies that LMXBs that ac-crete close to the Eddington limit are potentially strong grav-itational wave emitters. Of these potentially strong sources,Sco X-1 is the most promising due to its observed X-ray flux[28]. III. SEARCH ALGORITHM
In this section we define our detection statistic and showhow it exploits the characteristic frequency-modulation pat-tern inherent to sources in binary systems.Fully-coherent, matched-filter searches for continuousgravitational waves can be described as a procedure that max-imizes the likelihood function over a parameter space. Theamplitude parameters (gravitational wave strain amplitude h ,inclination ι , polarization ψ and reference phase φ ) can, ingeneral, be analytically maximized, reducing the dimensionsof this parameter space [54]. These parameters define thesignal amplitudes in our signal model. Analytic maximizationleaves the phase-evolution parameters (gravitational wave fre-quency f and its derivatives f ( k ) and sky position [ α, δ ]) to benumerically maximized over. Numerical maximization is ac-complished through a scheme of repeated matched-filteringperformed over a template bank of trial waveforms definedby specific locations in the phase parameter space, which istypically highly computationally expensive [28, 41, 53].The method we outline here makes use of the fact that weknow the sky position of our potential sources and, hence, thephase evolution due to the motion of the detector can be accu-rately accounted for. We also know that the phase evolutiondue to the binary motion of the source will result in a specificdistribution of signal power in the frequency domain. Thisdistribution has the characteristics that signal power is dividedamongst a finite set of frequency-modulation sidebands. Thenumber of sidebands and their relative frequency spacing canbe predicted with some knowledge of the binary orbital pa-rameters.In order to avoid the computational limits imposed by afully-coherent parameter space search, we propose a singlefully-coherent analysis stage, that accounts for detector mo-tion only, is followed by a single incoherent stage in whichthe signal power contained within the frequency-modulatedsidebands in summed to form a new detection statistic. Thissumming procedure is accomplished via the convolution of anapproximate frequency domain power template with the out-put of the coherent stage.The three main stages of the search, the F -statistic, side-band template, and C -statistic, are graphically illustratedin Figure 1. In this noise-free example, the frequency-modulation sidebands are clearly visible. The F -statisticis also amplitude-modulated due to the daily variation ofthe detector antenna response, resulting in the amplitude-modulation applied to each frequency-modulation sideband.The second panel represents the approximate frequency do-main template, a flat comb function with unit amplitude teeth(the spikes or delta functions). When convolved with the F -statistic in the frequency domain we obtain the C -statisticshown in the right-hand panel. The maximum power is clearlyrecovered at the simulated source frequency.The following sub-sections discuss each of the search com-ponents in more detail. Section III A presents the phasemodel used to characterize the gravitational wave signal froma source in a binary system. The signal model is introducedin Sec. III B. The F -statistic is introduced in Section III C andits behavior as a function of search frequency is described inSec. III D. Section III E then goes on to describe how match-ing a filter for an isolated neutron star system to a signal from asource in a binary system results in frequency-modulated side-bands appearing in the output of the F -statistic. The detec-tion statistic for gravitational wave sources in binary orbits isfully described as the incoherent sum of frequency-modulated F -statistics in Sec. III F. The simplest, unit amplitude, side-band template, and its justification over a more realistic tem-plate, are discussed in Sec. III G. A more sensitive implemen-tation, incorporating an approximate binary phase model inthe calculation of the F -statistic and reducing the width of thefrequency-modulated sideband pattern by the fractional errorson the semi-major axis parameters, is discussed in Sec. III H.Beginning in this section and continuing in the following sec-tions, we have made a distinction between the intrinsic valuesof a search parameter θ (denoted with a subscript zero) andthe observed values θ (no subscripts). A. Phase model
The phase of the signal at the source can be modeled by aTaylor series such that Φ ( τ ) = π s (cid:88) k = f ( k ) ( k + τ − τ ) k + , (8)where f ( k ) represents the k th time derivative of the gravita-tional wave frequency evaluated at reference time τ . For thepurposes of this work we restrict ourselves to a monochro-matic signal and hence set f ( k ) = k > f = f (0) | τ = f ( τ ) as the intrinsic frequency. We discussthis choice in Sec. IV B. The phase received at the detectoris Φ [ τ ( t ( t (cid:48) ))], where we define the retarded times measuredat the Solar system barycenter (SSB) and detector as t and t (cid:48) respectively. The relation between t and t (cid:48) is a function ofsource sky position relative to detector location and, since weonly concern ourselves with sources of known sky position,we assume that the e ff ects of phase contributions from detec-tor motion can be exactly accounted for. For this reason wework directly within the SSB frame. The relationship betweenthe source and retarded times for a non-relativistic eccentric binary orbit is given by [64] t = τ + a (cid:104) √ − e cos ω sin E + sin ω (cos E − e ) (cid:105) , (9)where a is the light crossing time of the semi-major axis pro-jected along the line of sight. The orbital eccentricity is de-fined by e and the argument of periapse, given by ω , is theangle in the orbital plane from the ascending node to the di-rection of periapse. The variable E is the eccentric anomalydefined by 2 π ( τ − t p ) / P = E − e sin E , where t p is the time ofpassage through periapse (the point in the orbit where the twobodies are at their closest) and P is the orbital period.It is expected that dissipative processes within LMXBsdrive the orbits to near circularity. In the low eccentricity limit e (cid:28)
1, we obtain the following expression t = τ + a (cid:26) sin [ Ω ( τ − t a )] + e cos ω (cid:2) Ω ( τ − t p ) (cid:3) + e sin ω (cid:2) Ω ( τ − t p ) (cid:3)(cid:41) + O ( e ) (cid:41) , (10)where Ω = π/ P and we have used the time of passagethrough the ascending node t a = t p − ω/ Ω as our referencephase in the first term. For circular orbits, where e =
0, theexpression reduces to only this first term. Using t a is sensiblein this case since t p and ω , are not defined in a circular orbit.Note that the additional eccentric terms are periodic at multi-ples of twice the orbital frequency. Using Eq. 10, we wouldexpect to accumulate timing errors of order µ s for the most ec-centric known LMXB systems. We shall return to this featurein Sec. IV E.To write the gravitational wave phase as a function of SSBtime, we invert Eq. 10 to obtain τ ( t ). The binary phase canbe corrected for in a standard matched filter approach, whereEq. 9 is solved numerically. In our method we instead approx-imate the inversion as Φ ( t ) (cid:39) π f ( t − t ) − π f a sin [ Ω ( t − t a )] (11)for circular orbits. Under this approximation the signal phasecan be represented as a linear combination of the phase con-tributions from the spin of the neutron star φ spin and from thebinary orbital motion of the source φ bin , such that φ spin = π f ( t − t ) , (12a) φ bin = − π f a sin [ Ω ( t − t a )] . (12b) B. Signal model
We model the data x ( t ) collected by a detector located atthe SSB as the signal s ( t ) plus stationary Gaussian noise n ( t )so that x ( t ) = s ( t ) + n ( t ) (13) Phase errors caused by this inversion approximation amount to maximumphase o ff sets of ∼ π f a Ω . F − statistic Template C − statistic FIG. 1: Graphical illustration of the sideband search pipeline, showing the frequency-modulated F -statistic(left, red), the sideband template (middle, black), and their convolution, known as the C -statistic (right,blue). In this noise free case, a signal of f =
200 Hz with amplitude h = a = .
005 s, P = . ↵ = . = . = . , cos ◆ = . , = presents the phase model used to characterize the gravitational wave signal from a source in abinary system. The signal model is introduced in Sec. III B. The F -statistic is introduced in Sec-tion III C and its behavior as a function of search frequency is described in Sec. III D. Section III Ethen goes on to describe how matching a filter for an isolated neutron star system to a signal froma source in a binary system results in frequency-modulated sidebands appearing in the output ofthe F -statistic. The detection statistic for gravitational wave sources in binary orbits is fully de-scribed as the incoherent sum of frequency-modulated F -statistics in Sec. III F. The simplest, unitamplitude, sideband template, and its justification over a more realistic template, are discussed inSec. III G. A more sensitive implementation, incorporating an approximate binary phase modelin the calculation of the F -statistic and reducing the width of the frequency-modulated sidebandpattern by the fractional errors on the semi-major axis parameters, is discussed in Sec. III H. Begin-ning in this section and continuing in the following sections, we have made a distinction betweenthe intrinsic values of a search parameter ✓ (denoted with a subscript zero) and the observed values ✓ (no subscripts). 9 FIG. 1: Graphical illustration of the sideband search pipeline, showing the frequency-modulated F -statistic (left, red), the sideband template(middle, black), and their convolution, known as the C -statistic (right, blue). In this noise free case, a signal of f =
200 Hz with amplitude h = a = .
005 s, P = . α = . δ = . ψ = . , cos ι = . , φ = with s ( t ) = A µ h µ ( t ) , (14)where we employ the Einstein summation convention for µ = , , ,
4. The coe ffi cients A µ are independent of time,detector location and orientation. They depend only on thesignal amplitude parameters λ = { h , ψ, ι, φ } , where h isthe dimensionless gravitational wave strain amplitude, ψ is thegravitational wave polarization angle, ι is the source inclina-tion angle and φ is the signal phase at a fiducial referencetime. The coe ffi cients A µ are defined as A = A + cos φ cos 2 ψ − A × sin φ sin 2 ψ (15a) A = A + cos φ sin 2 ψ − A × sin φ cos 2 ψ (15b) A = A + sin φ cos 2 ψ − A × cos φ sin 2 ψ (15c) A = A + sin φ sin 2 ψ − A × cos φ cos 2 ψ (15d)where A + = h (1 + cos ι ) , A × = h cos ι (16)are the polarization amplitudes. The time dependent signalcomponents h µ ( t ) are defined as h = a ( t ) cos Φ ( t ) , h = b ( t ) cos Φ ( t ) , h = a ( t ) sin Φ ( t ) , h = b ( t ) sin Φ ( t ) , (17)where Φ ( t ) is the signal phase at the detector (which we modelas located at the SSB) given by Eq.11 and the antenna patternfunctions a ( t ) and b ( t ) are described by Eqs. (12) and (13)in [54]. C. F -statistic The F -statistic is a matched-filter based detection statisticderived via analytic maximization of the likelihood over un-known amplitude parameters [54]. Let us first introduce the multi-detector inner product( x | y ) = (cid:88) X ( x X | y X ) = (cid:88) X S X ( f ) ∞ (cid:90) −∞ w X ( t ) x X ( t ) y X ( t ) dt , (18)where X indexes each detector and S X ( f ) is the detectorsingle-sided noise spectral density. We modify the definitionsof [54] and [65] to explicitly include gaps in the time-series byintroducing the function w X ( t ) which has value 1 when data ispresent and 0 otherwise. This also allows us to extend the lim-its of our time integration to ( −∞ , ∞ ) since the window func-tion will naturally account for the volume and span of data foreach detector.The F -statistic itself is defined as2 F = x µ M µν x ν , (19)where M µν form the matrix inverse of M µν and we follow theshorthand notation of [65] defining x µ ≡ ( x | h µ ) and M µν ≡ ( h µ | h ν ). Evaluation of M leads to a matrix of the form M = (cid:32) C C (cid:33) , where C = (cid:32) A CC B (cid:33) (20)where the components A = ( a | a ) , B = ( b | b ) , C = ( a | b ) , (21)are antenna pattern integrals. For a waveform with exactlyknown phase evolution Φ ( t ) in Gaussian noise, the F -statisticis a random variable distributed according to a non-central χ -distribution with 4 degrees of freedom. The non-centralityparameter is equal to the optimal signal-to-noise ratio (SNR) ρ = (cid:104) A ( A + A ) + B ( A + A ) + C ( A A + A A ) (cid:105) (22)such that the expectation value and variance of 2 F are givenby E[2 F ] = + ρ (23a)Var[2 F ] = + ρ , (23b)respectively. In the case where no signal is present in the data,the distribution becomes a central χ -distribution with 4 de-grees of freedom. D. F -statistic and mismatched frequency In this section we describe the behavior of the F -statistic asa function of search frequency f for a fixed source frequency f . In this case the inner product that defines x µ becomes x µ = A ν ( h ν | h µ (cid:48) ) + ( n | h µ (cid:48) ) , (24)where h ν are the components of a signal with frequency f ,and h (cid:48) µ is a function of search frequency f . If we focus on the µ = ν = h | h (cid:48) ) (cid:27) (cid:88) X S X ( f ) ∞ (cid:90) −∞ w X ( t ) a X ( t ) cos (2 π f t ) cos (2 π f t ) dt (25)where we note that the product of cosine functions results inan integrand that contains frequencies at f − f and f + f . Sinceboth a X ( t ) and w X ( t ) are functions that evolve on timescales ofhours–days we approximate the contribution from the f + f component as averaging to zero. We are left with( h | h (cid:48) ) (cid:27) Re (cid:88) X S X ( f ) ∞ (cid:90) w X ( t ) a X ( t ) e − π i ( f − f ) t dt (cid:27) (cid:88) X Re (cid:2) A X ( f − f ) (cid:3) (26)where we have defined the result of the complex integral as A X ( f ). This is the Fourier transform of the antenna patternfunctions weighted by the window function and evaluated at f − f . It is equal to the quantity ( a X | a X ) when its argument iszero.The quantity a X ( t ) (and similarly b X ( t ) and a X ( t ) b X ( t )) areperiodic quantities with periods of 12 and 24 hours plus anon-oscillating component. When in a product with a sinu-soidal function and integrated over time they will result indiscrete amplitude-modulated sidebands with frequencies at0 , ± / P ⊕ , ± / P ⊕ , ± / P ⊕ , ± / P ⊕ Hz where P ⊕ represents theorbital period of the Earth (1 sidereal day). We will ignoreall but the zero-frequency components of these functions forthe remainder of this paper. We do note that complicationsregarding the overlap of amplitude-modulated and frequency-modulated sidebands (discussed in the next section) will onlyarise for sources in binary orbits with periods equal to thosepresent in the antenna pattern functions. In addition, the window function describing the gaps in thedata will influence A X ( f ). For a gap-free observation the win-dow function serves to localize signal power to within a fre-quency range ∼ / T where T is the typical observation length.When gaps are present this range is broadened and has a de-terministic shape given by the squared modulus of the Fouriertransform of the window function. We can therefore use afurther approximation that A X ( f − f ) ≈ ( a X | a X ) ˜ w X ( f − f )˜ w X (0) (27)where the Fourier transform of the window function is nor-malized by ˜ w X (0) ≡ (cid:82) dtw X ( t ) such that it has a value of unityat the true signal frequency.We now define the antenna-pattern weighted window func-tion as ˜ W ( f − f ) (cid:27) (cid:88) X ( a X | a X ) A ˜ w X ( f − f )˜ w X (0) (28a) (cid:27) (cid:88) X ( b X | b X ) B ˜ w X ( f − f )˜ w X (0) (28b) (cid:27) (cid:88) X ( a X | b X ) C ˜ w X ( f − f )˜ w X (0) , (28c)which is true for observation times T X (cid:29) days. This complexwindow function has the property that ˜ W (0) is a real quan-tity with maximum absolute value of unity when the templatefrequency matches the true signal frequency.Finally we are able to combine Eqs. 26, 27, and 28 to obtain( h | h (cid:48) ) (cid:27) A (cid:104) ˜ W ( f − f ) (cid:105) (29)which together with similar calculations for the additionalcomponents in Eq. 24 give us( h ν | h µ (cid:48) ) (cid:27) M (cid:40)(cid:32) I 00 I (cid:33) Re (cid:104) ˜ W ( f − f ) (cid:105) + (cid:32) − I 0 (cid:33) Im (cid:104) ˜ W ( f − f ) (cid:105)(cid:41) (30)as the complete set of inner products between frequency-mismatched signal components where I is the 2 × f = f this expression reduces toEq. 19.If we now form the expectation value of the F -statistic(Eq. 19) for mismatched frequencies we find thatE[2 F ( f )] = + ρ | ˜ W ( f − f ) | . (31)Here we see that the fraction of the optimal SNR that con-tributes to the non-centrality parameter of the F -statistic χ distribution is reduced by evaluation of the mod-squared ofthe antenna-pattern weighted window function with a non-zero argument. E. Frequency-modulation and the F -statistic We now consider the computation of the F -statistic in thecase where the data contains a signal from a source in a cir-cular binary orbit but the phase model used in the F -statistictemplate is that of a monochromatic signal of frequency f . Weagain expand x µ as done in Eq. 24 where no prime indicatesthe signal and the prime represents the monochromatic tem-plate. We again focus on the mismatched signal inner-product( h | h (cid:48) ) as an example. Starting with Eq. 25 we discard therapidly oscillating terms inside the integral that will averageto zero. We are then left with( h | h (cid:48) ) (cid:27) (cid:88) X Re (cid:40) S X ( f ) ∞ (cid:90) w x ( t ) a X ( t ) e − π it ( f − f ) exp ( − π f a sin ( Ω ( t − t a ))) (cid:41) (32)where the final term involving the exponential of a sinusoidalfunction can be represented using the Jacobi-Anger expansion e iz sin θ = ∞ (cid:88) n = −∞ J n ( z ) e in θ , (33)where J n ( z ) is the n th order Bessel function of the first kind.This expansion allows us to transform the binary phase terminto an infinite sum of harmonics such that we can now write( h | h (cid:48) ) (cid:27) ∞ (cid:88) n = −∞ J n (2 π f a ) × Re (cid:40) e i φ n (cid:88) X S X ( f ) ∞ (cid:90) w X ( t ) a X ( t ) e − π it ( f − f n ) dt (cid:41) (cid:27) ∞ (cid:88) n = −∞ J n (2 π f a ) A (cid:110) ˜ W ( f − f n ) e i φ n (cid:111) (34)It follows that all of the signal components can be expandedin the same way giving us( h ν | h µ (cid:48) ) (cid:27) M m (cid:88) n = − m J n (2 π f a ) (cid:40)(cid:32) I 00 I (cid:33) Re (cid:104) ˜ W ( f − f n ) e i φ n (cid:105) + (cid:32) − I 0 (cid:33) Im (cid:104) ˜ W ( f − f n ) e i φ n (cid:105)(cid:41) (35)where we have truncated the infinite summation (explainedbelow) and defined the monochromatic modulated sidebandfrequencies and their respective phases as f n = f − n / P , (36a) φ n = n Ω t a . (36b)The Jacobi-Anger expansion has allowed us to represent thecomplex phase of a frequency-modulated signal as an infinitesum of discrete signal harmonics, or sidebands, each sepa-rated in frequency by 1 / P Hz. Each is weighted by the Besselfunction of order n where n indexes the harmonics and has acomplex phase factor determined by the orbital reference time t a . In the limit where the order exceeds the argument, n (cid:29) z ,the Bessel function rapidly approaches zero allowing approx-imation of the infinite sum in Eqs. 33 and 34 as a finite sum over the finite range [ − m , m ] where m = ceil [2 π f a ]. Thesummation format of Eq. 35 highlights the e ff ects of the bi-nary phase modulation. The signal can be represented as thesum of M = m + f n cen-tered on the intrinsic frequency f , where each harmonic peakis separated from the next by 1 / P Hz.Combining Eqs. 19, 24, and 35, we can express the expec-tation value of the F -statistic for a binary signal as a functionof search frequency f as E [2 F ( f )] = + ρ m (cid:88) n = − m J n (2 π f a ) | ˜ W ( f − f n ) | . (37)This expression should be interpreted in the following man-ner. For a given search frequency f the contribution to thenon-centrality parameter (the SNR dependent term) is equal tothe sum of all sideband contributions at that frequency. Eachsideband will contribute a fraction of the total optimal SNRweighted by the n th order Bessel function squared, but willalso be strongly weighted by the window function. The win-dow function will only contribute significantly if the searchfrequency is close to the sideband frequency. Hence, at a givensearch frequency close to a sideband, for observation times (cid:29) P , the sidebands will be far enough separated in frequencysuch that only one sideband will contribute to the F -statistic. F. C -statistic The F -statistic is numerically maximized over the phaseparameters of the signal on a discrete grids. For this searchthe search frequency f is such a parameter and consequentlythe F -statistic is computed over a uniformly spaced set of fre-quency values f j spanning the region of interest. In this sec-tion we describe how this F -statistic frequency-series can beused to approximate a search template that is then used to gen-erate a new statistic sensitive to signals from sources in binarysystems.The expectation value of the F -statistic (Eq. 37) resolvesinto localized spikes at M frequencies separated by 1 / P Hzand centered on the intrinsic gravitational wave frequency f .A template T based on this pattern with amplitude defined by G n , takes the general form T ( f ) = m (cid:48) (cid:88) n = − m (cid:48) G n | ˜ W (cid:0) f − f (cid:48) n (cid:1) | , (38)with m (cid:48) = ceil [2 π f a (cid:48) ] and f (cid:48) n = f (cid:48) − n / P (cid:48) where we makea distinction between the intrinsic (unknown) values of eachparameter (subscript zero) and values selected in the templateconstruction (denoted with a prime). The window function ˜ W is dependent only on the times for which data is present andis, therefore, also known exactly.We define our new detection statistic C as C ( f ) ≡ (cid:88) j F ( f j ) T ( f j − f ) = (2 F ∗ T ) ( f ) , (39)where the sum over the index j indicates the sum over the dis-crete frequency bins f j and f is the search frequency. Sincethe template’s “zero frequency” represents the intrinsic grav-itational wave frequency, f corresponds to the intrinsic fre-quency. We see that the C -statistic is, in fact, the convolutionof the F -statistic with our template, assuming the templateis constant with search frequency (an issue we address in thenext section).The benefit of this approach is that the computation of the F -statistic for a known sky position and without accountingfor binary e ff ects has relatively low computational cost. Simi-larly, the construction of a template on the F -statistic is in-dependent of the orbital phase parameter and only weaklydependent upon the orbital semi-major axis and eccentric-ity. The template is highly dependent upon the orbital period,which, for the sources of interest, is known to high precision.Also, since the C -statistic is the result of a convolution, wecan make use of the convolution theorem and the speed of theFast-Fourier-Transform (FFT). Computing C for all frequen-cies requires only three applications of the FFT. In practice,the C -statistic is computed using C ( f k ) = (2 F ∗ T ) ( f k ) = N − (cid:88) j = e π i jk / N N − (cid:88) p = e − π ip j / N F ( f p ) N − (cid:88) q = e − π i jq / N T ( f q ) , (40)which is simply the inverse Fourier transform of the productof the Fourier transforms of the F -statistic and the template.The F -statistic and the template are both sampled on the sameuniform frequency grid containing N frequency bins. The C -statistic is then also output as a function of the same frequencygrid. G. Choice of F -statistic template Treating the F -statistic as the pre-processed input datasetto the C -statistic computation, it might be assumed that theoptimal choice of template is that which exactly matches theexpected form of 2 F in the presence of a signal. As shown inFig. 2, this approach is highly sensitive to the accuracy withwhich the projected orbital semi-major axis is known.We instead propose the use of a far simpler template: onethat captures the majority of the information contained withinthe F -statistic and, by design, is relatively insensitive to theorbital semi-major axis. We explicitly choose T F ( f k ) = m (cid:88) j = − m δ k l [ j ] , (41)for discrete frequency f k , where δ i j is the Kronecker delta-function. The frequency index l , defined by, l [ j ] ≡ round (cid:34) jP (cid:48) d f (cid:35) , (42) is a function of the best-guess orbital period P (cid:48) and the fre-quency resolution d f . The round [] function returns the inte-ger closest to its argument. The template is therefore com-posed of the sum of M = m + F -statistic we obtain the corresponding C -statistic,which reduces to C ( f k ) = m (cid:88) j = − m F ( f k − l [ j ] ) . (43)Equation 43 is simply the sum of 2 F values taken from dis-crete frequency bins positioned at the predicted locations ofthe frequency-modulation sidebands. An example is shown inthe right hand panel of Fig. 1, where the most significant C -statistic value is located at f = f , and is the point where allsidebands included in the sum contain some signal.From Eq. 37, the expectation value for the C -statistic usingthe flat-template can be expressed asE[ C ( f k )] = M + ρ m (cid:88) n = − m J n (2 π f a ) | ˜ W ( f n − f k + l [ n ] d f ) | , (44)where the argument of the window function is the frequencydi ff erence between the location of the n th signal sideband andthe n th template sideband on the discrete frequency grid. Notethat the C -statistic is a sum of M statistically independentnon-central χ statistics and hence the result is itself a non-central χ M statistic, i.e. with 4 M degrees of freedom, where M = ceil [2 π f a ] + F values. For a flat template with perfectlymatched intrinsic frequency f = f and orbital period P (cid:48) = P ,infinite precision d f →
0, and where the number of sidebandsincluded in the analysis matches or exceeds the true number,the second term in the above equation reduces to ρ . In thiscase we will have recovered all of the power from the signalbut also significantly increased the contribution from the noisethrough the incoherent summation of F -statistic from inde-pendent frequency bins. In general, where the orbital periodis known well, but not exactly, and the frequency resolution isfinite, the signal power recovery will be reduced by imperfectsampling of the window function term in Eq. 44, i.e. evalua-tion at arguments (cid:44) G n =
1. A more sensitive ap-proach could use T B ( f k ) = m (cid:88) j = − m J j (2 π f a ) δ k l [ j ] (45)for the template, following the expected form of the F -statistic given in Eq. 37 and using a subscript B to denoteBessel function weighting. Although this would increasesensitivity for closely matched signal templates (constructedwith well constrained signal parameters), this performance ishighly sensitive to the number of sidebands included in thetemplate and therefore sensitive to the semi-major axis since M = ceil [2 π f a ] +
1. This is mainly due to the “doublehorned” shape of the expected signal (see the left hand panelof Fig. 1). A large enough o ff set between the true and assumedsemi-major axis will significantly change the template’s over-lap with the sidebands in the F -statistic and reduce the sig-nificance of the C -statistic. Considering the semi-major axisis not well constrained for many LMXBs, a search over manytemplates would be necessary, each with incrementally di ff er-ent semi-major axis values.The simpler, flat-template (Eq. 41) has the benefit of be-ing far more robust against the semi-major axis uncertainty.In this case the semi-major axis parameter controls only thenumber of sidebands to use in the template and does not con-trol the weighting applied to each sideband. It also simplifiesthe statistical properties of the C -statistic, making a Bayesiananalysis of the output statistics (as described in Section VI)far easier to apply.The receiver operator characteristic (ROC) curves shown inFig. 2 compare the performance of the sideband search withboth choices of template (flat and Bessel function weighted)for the case of signal with optimal SNR ρ =
20. As seen fromthe figure, the Bessel function weighted template for exactnumber of sidebands provides improved sensitivity over theflat template. However, when considering the possible (andhighly likely) error on the number of sidebands in the tem-plate, the performance of the Bessel template is already dras-tically diminished, even with only a 10% error on the semi-major axis parameter. It is also interesting to note that forthe flat-template the result of an incorrect semi-major axis isasymmetric with respect to an under or over-estimate. Thesensitivity degradation is far less pronounced when the tem-plate has over-estimated the semi-major axis and, therefore,also over-estimated the number of sidebands. This feature isdiscussed in more detail in Sec. IV D.
H. Approximate binary demodulation
When a putative source has a highly localized position inthe sky, the e ff ect of the Earth’s motion with respect to theSSB can be accurately removed from the signal during the cal-culation of the F -statistic. This leaves only the Doppler mod-ulation from the binary orbit. It is also possible to demodulatethe binary orbit (Doppler) modulation in the F -statistic calcu-lation provided the binary orbital parameters ( a , P , t a ) are wellknown. A fully-coherent (sky position- and binary-) demod-ulated F -statistic search would be very sensitive to any er-rors in the sky position of binary orbital parameters. It wouldtherefore be necessary to construct a bank of templates span-ning the parameter space defined by the uncertainties in theseparameters. Adding dimensions to the parameter space in-creases computational costs and the search becomes unfea-sible considering we are already searching over frequency. A −3 −2 −1 False Alarm Probability D e t e c t i on P r obab ili t y χ ( λ )flatbessel − 10% + 10% FIG. 2: ROC curves comparing performance of the flat (blue) andBessel function weighted (red) templates, described by Eqs 41 and 45respectively. The theoretical (fine black) curve is constructed from anon-central χ M ( λ ) distributed statistic with non-centrality parameter λ = ρ =
20 and represents the expected performance of the flat tem-plate. Dashed and dotted curves represent a template with a positiveand negative 10% error on the semi-major axis respectively. Herethe signal parameters were chosen such that the number of sidebandswere M = realizationsof noise. fully coherent search of this type would be possible for knownsources with known emission frequencies, for example puls-ing sources like millisecond pulsars (MSPs).In this section we show how prior information regarding thebinary orbit of a source can be used to increase the sensitivityof our semi-coherent approach, without increasing computa-tional costs. By performing a “best guess” binary phase de-modulation within the F -statistic, we show that the numberof sidebands in the template is reduced by a factor propor-tional to the fractional uncertainty in the orbital semi-majoraxis. Consequently a reduction in the number of sidebandsincreases the sensitivity of the search by reducing the numberof 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 θ = θ + ∆ θ , (46)we can determine the phase o ff set of the binary orbit from theerror in the binary orbital parameters. The o ff set in phase isthe di ff erence between the true and best estimate binary phaseand using Eq. 12b can be approximated by ∆ φ bin (cid:39) − π f (cid:26) ( a ∆ f + f ∆ a ) sin [ Ω ( t − t a )] + (cid:2) f a ( ∆Ω ( t − t a ) − Ω ∆ t a ) (cid:3) cos [ Ω ( t − t a )] (cid:27) + O ( ∆ θ ) (cid:39) − π f ∆ a κ sin (cid:2) Ω ( t − t a ) + γ (cid:3) , (47)0with κ = (cid:115)(cid:32) + a ∆ a ∆ ff (cid:33) + (cid:32) a ∆ a ( ∆Ω ( t − t a ) − Ω ∆ t a ) (cid:33) , (48) γ = tan − (cid:34) ∆Ω ( t − t a ) − Ω ∆ t a ( ∆ a / a ) − ( ∆ f / f ) (cid:35) + (cid:16) ∆ a a (cid:17) − (cid:16) ∆ ff (cid:17) ≥ π if (cid:16) ∆ a a (cid:17) − (cid:16) ∆ ff (cid:17) < . (49)Here we have expanded the binary phase di ff erence to lead-ing order in the parameter uncertainties and obtained a phaseexpression similar in form to the original binary phase. Inthe specific regime where the fractional uncertainty in the or-bital semi-major axis far exceeds the fractional uncertainty inthe intrinsic frequency we see that the first term in Eq. 48 be-comes ≈
1. Similarly, if the fractional uncertainty in the or-bital semi-major axis also far exceeds the deviation in orbitalangular position ∆Ω ( t − t a ) − Ω∆ t a then the second term ≈ κ can be accurately approximated as unity,yielding ∆ φ bin ≈ − π f ∆ a sin[ Ω ( t − t a ) + γ ] . (50)Hence, after approximate binary demodulation, the argumentof the Bessel function and the summation limits in the ex-pected form of the F -statistic (in Eq. 37 for example) canbe replaced with ∆ z = π f ∆ a . The number of frequency-modulated sidebands is now reduced by a factor of ∆ a / a <
1. We must stress that ∆ a is an unknown quantity and is thedi ff erence between the best estimate value of a and the truevalue a . The F -statistic after such a demodulation processwill therefore have a reduced but unknown number of side-bands, although it will still retain the standard sideband fre-quency spacing 1 / P . The sideband phasing φ n will also beunknown due to the presence of the phase term γ but is of noconsequence to the search since the F -statistic is insensitiveto phase. IV. PARAMETER SPACE
In this section we will discuss each of the parameters in-volved in the search and how the search sensitivity dependsupon the uncertainty in these parameters. Demodulation ofthe signal phase due to the Earth’s motion requires accurateknowledge of the source sky position. If the observation timeis long enough, we need to consider the sky position as asearch parameter, as discussed in Sec. IV A. The gravitationalwave frequency is the primary search parameter. In Sec. IV Bwe discuss the limitations on our search strategy due to itsuncertainty. The orbital period and semi-major axis are dis-cussed in Sections IV C and IV D respectively. The e ff ects ofignoring the orbital eccentricity are discussed in Sec. IV E. A. Sky position and proper motion
In order to quantify the allowable uncertainty in sky posi-tion we will define a simplistic model describing the phase Ψ ( t ) received at Earth from a monochromatic source at infin-ity at sky position ( α, δ ). If we neglect the detector motion dueto the spin of the Earth and consider only the Earth’s orbitalmotion then we have Ψ ( t ) = π f [ t + R ⊕ cos δ cos ( Ω ⊕ ( t − t ref ) + α )] , (51)where f is the signal frequency, α and δ are the true rightascension and declination and R ⊕ and Ω ⊕ are the distanceof the Earth from the SSB and the Earth’s orbital angularfrequency respectively. We also define a reference time t ref that represents the time at which the detector passes throughthe vernal equinox. For an observed sky position ( α (cid:48) , δ (cid:48) ) = ( α +∆ α, δ +∆ δ ) the corresponding phase o ff set ∆Ψ ( t , ∆ α, ∆ δ ) =Ψ ( t , α (cid:48) , δ (cid:48) ) − Ψ ( t , α, δ ) amounts to ∆Ψ ≈ − π f R ⊕ (cid:104) ∆ δ sin δ cos ( Ω ⊕ ( t − t ref ) + α ) +∆ α cos δ sin ( Ω ⊕ ( t − t ref ) + α ) (cid:105) , (52)where we have expanded the expression to leading order in thesky position errors. We now make the reasonable assumptionthat our analysis would be unable to tolerate a deviation inphase between the signal and our template of more than O (1)radian over the course of an observation on the same timescaleof the Earth’s orbit.If we also notice that the worst case scenario (smallest al-lowable sky position errors) corresponds to sky positions forwhich the trigonometric terms in the previous expression arelargest, i.e. of order unity, then we have | ∆ α, ∆ δ | ≤ (2 π f R ⊕ ) − . (53)If we consider signals of frequency 1kHz, this gives a max-imum allowable sky position o ff set of | ∆ α, ∆ δ | (cid:39)
100 mas.This expression also validates our model assumption that thesky position sensitivity to the Earth spin would be dominatedby the e ff ect from the Earth orbit for long observation times.A similar argument can be made for the proper motion ofthe source where we would be safe to model the sky positionas fixed if the change ( ∆ α, ∆ δ ) = ( µ α , µ δ ) T s , over the courseof the observation also satisfied Eq. 53. B. Spin frequency
The spin frequency ν s of some LMXBs can be directly mea-sured from X-ray pulsations, believed to originate from a hot-spot on the stellar surface, where accreted material is funneledonto the magnetic pole with the magnetic axis generally mis-aligned with the spin axis. X-ray pulsations have been ob-served in 13 LMXB systems so far, three of which are inter-mittent [66].Some LMXBs exhibit recurrent thermonuclear X-raybursts. Fourier spectra reveal oscillations during the rise andtail of many bursts, which are believed to originate from asym-metric brightness patterns on the stellar surface. In sevenLMXBs which exhibit both pulsations and bursts, the asymp-totic burst oscillation frequency at late times matches the pulse1frequency. Where there are no pulsations, many bursts needto be observed to measure the asymptotic burst oscillation re-liably. The spin frequency of an additional ten systems hasbeen determined from burst oscillations only [67], but due tothe uncertainties involved, are usually quoted to within uncer-tainties of ± (5 −
10) Hz.Another class of LMXBs exhibit high frequency quasi-periodic oscillationss (QPOs) in their persistent X-ray emis-sion. These kHz QPOs usually come in pairs, although sin-gles and triples are occasionally observed and the QPO peakfrequencies usually change over time. In some cases the sep-aration of the QPO peaks is roughly constant, but this is notalways the case [68–70]. For the few QPO systems where ν s can be determined from pulses or burst oscillations therehas been no evidence suggesting consistency with an existingQPO model that links the QPO and spin frequencies. For ourpurposes, ν s is considered unknown in sources without pulsa-tions or confirmed bursts.In addition to potentially broad uncertainties in ν s , we knowlittle about how its value may fluctuate over time due to accre-tion. Changes in the accretion flow will exert a time varyingtorque on the star which will result in a stochastic wanderingof the spin frequency. In this case the signal can no longerbe assumed monochromatic over a given observation time. Toquantify the resulting phase wandering, we assume that thefluctuating component of the torque δ N a flips sign randomlyon the timescale τ s consistent with the inferred variation inaccretion rate. If the mean torque N a = ˙ M ( GMR ) / dueto steady-state disk-fed accretion, then the angular spin fre-quency Ω s = πν s experiences a random walk with step size( δ N a / I ) τ s , where I is the stellar moment of inertia. After time T s , the root-mean-square drift is (cid:104) ( δ Ω ) (cid:105) / = ( T s /τ s ) / δ N a τ s I . (54)This frequency drift will wander outside a Fourier frequencybin width if (cid:104) ( δ Ω ) (cid:105) / > π/ T s . If we choose τ s such that theaccretion rate can vary up to a factor of two in this time, thenthe worst case δ N a = N a leads to the restriction T spins < (2 π ) / ( GMR ) / (cid:18) I ˙ M (cid:19) / (cid:32) τ s (cid:33) / . (55)This is the primary reason why an application of the the ba-sic sideband search, as described here, must be limited in thelength of data it is allowed to analyse. By exceeding this limitit becomes increasingly likely that the spin wandering inher-ent to a true signal will cause signal power to leak betweenadjacent frequency bins. Consequently the assumption that F -statistic signal power is localized in frequency-modulatedsidebands will become invalid and the sensitivity of the C -statistic will deteriorate. C. Orbital Period
The sideband search relies on relatively preciseelectromagnetic (EM) measurements of the orbital pe-riod in order to construct a search template. The duration of the orbit defines the minimum observation time, since T (cid:38) P / P .We will now provide an indication of the sensitivity of thesearch to errors in the orbital period. If our estimate (obser-vation) of the orbital period P is o ff set from the true value P by an amount ∆ P , we would expect the error to seriously af-fect the C -statistic recovered from the search once it is largeenough to shift the outermost “tooth” in the sideband tem-plate by one canonical frequency bin away from the true side-band location. In this case, the o ff set between the templateand true sideband frequency is proportional to the number ofsidebands from the central spike. There will be low mismatchat the center of the template extending to O (100%) mismatchat the edges. It follows that the average signal power recov-ered from such a mismatched template will be O (50%) andtherefore serves as a useful threshold by which to determinethe maximum allowed ∆ P .If we use the measured value of P as our template param-eter, the template centered at frequency f then consists of (cid:39) π f a unit spikes (or teeth) separated by 1 / P . Assumingthat the central spike is exactly equal to the true intrinsic grav-itational wave frequency, any errors in the orbital period willbe propagated along the comb, causing the o ff set between thetrue and template frequency of any particular sideband to growprogressively larger. The frequency di ff erence ∆ f P betweenthe outermost template sideband, at frequency f + π f a / P (cid:48) ,and the outermost signal sideband at f + π f a / P , is given by ∆ f P ≈ π f a (cid:32) | ∆ P | P (cid:33) , (56)for ∆ P (cid:28) P . To satisfy the condition described above wenow require that this frequency shift should be less than thesize of one frequency bin. The true frequency bin size d f isdetermined by the observation time span and is given by d f = rT s , (57)where r is the resolution used in the F -statistic calculation. Using Eqs. 56 and 57 and imposing the condition that ∆ f P < d f provides an estimate for the maximum allowable (orbitalperiod limited) coherent observation timespan, T P s ≈ P π f a | ∆ P | . (58)Given a relatively poorly constrained orbital period uncer-tainty, this restriction may provide too short a duration. Thiscould be because it is then in conflict with the requirementthat T > P / The default resolution factor is r = δ P derived from sim-ply rearranging Eq. 58 to solve for ∆ P . In relative terms thesideband search places very strong constraints on the priorknowledge of the orbital period compared to the other searchparameters. D. Semi-major axis
An error in the value of the orbital semi-major axis resultsin an incorrect choice for the number of sidebands in the tem-plate. As can be seen in Eq. 44, an underestimate results inthe summation of a fraction of the total power in the signalwhereas an overestimate results in a dilution of the total powerby summing additional noise from sideband frequencies con-taining no signal contribution.If we define the true semi-major axis parameter a as themeasured value a and some (unknown) fraction ξ (where ξ ∈ R ) of the measurement error ∆ a (i.e. a = a + ξ ∆ a ),we can investigate the e ff ects of errors on the semi-major axisparameter in terms of this o ff set parameter ξ . We considerthe advantage of using a deliberately o ff set value a (cid:48) instead ofthe observed value a in order to minimize losses in recoveredSNR.The ROC curves shown in Fig. 3 show the e ff ects of theseo ff sets, and clearly illustrate degradation in the performanceof the C -statistic as | ∆ a | ( | ξ | ) increases. The reduction in de-tection confidence at a given false alarm probability is muchfaster for a (cid:48) < a ( ξ < a (cid:48) > a ( ξ > a (cid:48) = a , this asymmetry ismuch better illustrated in Fig. 4 where for di ff erent values ofthe false alarm rate we show the detection probability plottedagainst the o ff set parameter ξ .This plot provides us with a rough scheme by which to im-prove the search performance by exploiting the asymmetry insearch sensitivity with respect to ξ . In general, we are keento probe the low false alarm and high detection probabilityregime in which it is clear that using a template based on anorbital semi-major axis value > the best estimate reduces thepossibility that the bulk of the signal power (in the horns) willbe missed. Based on Fig. 4 we choose a (cid:48) = a + ∆ a (59)as our choice of semi-major axis with which to generate thesearch template. E. Orbital eccentricity
The orbital eccentricity e of the LMXB sources is expectedto be highly circularized ( e < − ) by the time mass transferoccurs within the system. In Eq. 10 we give the first-order −3 −2 −1 False Alarm Probability D e t e c t i on P r obab ili t y −3−2−1 0+1+2+3 Offset Parameter
FIG. 3: ROC curves showing how the performance of the flat tem-plate C -statistic is a ff ected by an o ff set in the orbital semi-major axisassuming it is measured exactly (i.e. a = a ). The thick black curverepresents a zero o ff set ( ξ = ff set in the semi-major axis ( ξ > ff sets ( ξ < χ M ( λ ) dis-tribution with a non-centrality parameter λ = ρ and represents thetheoretical expected behavior of a perfectly matched template. Sig-nal parameters are the same as described in Fig. 2, with ρ = M = realizations of noise. Offset Parameter D e t e c t i on P r obab ili t y −5 −4 −3 −2 −1 0 1 2 3 4 50.10.20.30.40.50.60.70.80.91 Fa −3−2.5−2−1.5−1−0.50 FIG. 4: Performance of the C -statistic with respect to o ff sets in thesemi-major axis at given false alarm probability F a . The o ff set pa-rameter ξ quantifies the semi-major axis error in terms of known pa-rameters a and ∆ a . The color-bar represents the log of the false alarmprobability, ranging from 10 − (bottom, blue) to 0 (top, red). Signalparameters are the same as Figs. 2 and 3, with ρ = M = ∆ a / a =
10% fractional uncertainty on the semi-major axis,using 10 realizations of noise. e ) of the retarded time at the SSB.If we include higher order terms in the expansion, the phase(Eq. 11) can be written as Φ ( t ) (cid:39) π f a ∞ (cid:88) k = c k sin ω cos (cid:2) k Ω (cid:0) t − t p (cid:1)(cid:3) + d k cos ω sin (cid:2) k Ω (cid:0) t − t p (cid:1)(cid:3) (cid:41) , (60)where the first 4 coe ffi cients (expanded to O ( e )) in the sumare given by c = − e + e + O ( e ) (61a) d = − e − e + O ( e ) (61b) c = e − e + O ( e ) (61c) d = e − e + O ( e ) . (61d)Hence the phase for e (cid:44) i (cid:88) k z k sin k θ = ∞ (cid:89) k = ∞ (cid:88) n = −∞ J n ( z k ) e ink θ , (62)where z k corresponds to the k th amplitude term (on the lefthand side) that defines the argument of the Bessel function foreach k in the product (on the right hand side). Equation 62 tellsus that eccentric signals can be thought of in a similar wayto circular orbit cases. The signal can be modeled as beingcomposed of many harmonics all separated by some integernumber of the inverse of the orbital period. In the eccentriccase k is allowed to be > k = ,
2. The k = ≈ z sidebands at frequencies o ff set from the intrin-sic source frequency by integer multiples of θ .In our low eccentricity case we notice that the next to lead-ing order term in the expansion, k =
2, has a correspondingBessel function argument of O ( z e ) and will therefore havefar fewer, ∼ z e , non-negligible terms in sum over n . Takingthe product between the k = k = ≈ z side-bands whereas we now expect the same power to be dividedamongst ≈ z (1 + e ) sidebands. In general, orbital eccentricity causes a redistribution of sig-nal power amongst the existing circular orbit sidebands andwill cause negligible leakage of signal power into additionalsidebands at the boundaries of the sideband structure. Orbitaleccentricity also has the e ff ect of modifying the phase of eachsideband. However, as shown in Section III E, the standardsideband search is insensitive to the phase of individual side-bands. V. PRIMARY SOURCES
The benefit of the sideband search is that it is robustand computationally cheap enough to be run over a widefrequency band [59]. The most suitable targets are thosewith well-measured sky position and orbital periods, reason-ably well constrained semi-major axes, and poorly or uncon-strained spin frequency.The most suitable candidates in terms of these criteria areLMXBs due to their high accretion rate (directly related togravitational wave amplitude) and their visibility in the elec-tromagnetic regime (predominantly X-ray, but optical and ra-dio observations also provide accurate sky position, ephemerisand sometimes orbital information). They are classified intothree main types depending on the behavior of their X-rayemissions: pulsing, bursting or QPO sources. Pulsatingand frequently bursting LMXBs usually have a well deter-mined spin frequency and are better suited to the more sen-sitive, narrow-band techniques, such as LIGO’s known pul-sar pipeline [43, 44] including corrections for the binary mo-tion. Non-pulsing burst sources with irregular or infrequentbursts still have a fairly wide ( O a few Hz) range around thesuspected spin frequency. A convincing relationship betweenQPOs and the spin frequency of the neutron star has not yetbeen determined, so the spin frequency of purely QPO sourcesis considered unknown for our applications.The gravitational wave strain amplitude is directly pro-portional to the square root of the x-ray flux h ∝ F / X (Eq. 7), so the most luminous sources, which are usually thequasi-periodic oscillation (QPO) sources, will be the most de-tectable. In addition, the (already weak) gravitational wavestrain amplitude is proportional to the inverse of the distanceto the source, so closer (i.e. galactic) sources are also favor-able.In this section we present possible sources to which thesideband search can be applied. We start with galacticLMXBs and consider the most detectable sources in terms oftheir parameter constraints. We exclusively consider sourcesrequiring wide frequency search bands ( (cid:38) TABLE I: Target sources for the sideband method. The di ff erent columns list the X-ray flux F X (in units of 10 − erg cm − s − ), distance d , skyposition uncertainty ∆ β , fractional error on the semi-major axis ∆ a / a and orbital period ∆ P / P , and the orbital period limited observation timeat a frequency of 1 kHz T P s | . The horizontal line separates QPO (top) and burst (bottom) sources.Source F X ( F ∗ ) a d (kpc) ∆ β (arc sec) ∆ a / a ( ∆ P / P )( × − ) T P s | Sco X-1 40 2.8 3 × − <
60 0.11 10 17 daysX 1658-298 0.67 12 0.1 0.82 0.1 5 yearsXB 1254-690 0.09 13 0.6 0.12 500 6 hoursEXO 0748-676 0.036 7.4 0.7 0.77 6 40 days4U 1916-053 0.027 8 0.6 0.72 3 2 years a F ∗ = − erg cm − s − A. Galactic LMXBs
The sideband search is best suited to LMXBs with a rela-tively large uncertainty in the spin frequency ( (cid:38) a few Hz),so QPO and poorly constrained burst sources are the best tar-gets. The requirement of a relatively well defined sky positionand orbital period excludes many sources including those thatare considered to be X-ray bright. Table I lists some of thegalactic LMXBs, and their limiting parameters, for which thesideband search is most applicable. The parameters displayedin the table allow us to determine the most suitable targets forthe search.For each source the table lists the bolometric X-ray flux F X , the distance to the source d , the error in the sky posi-tion ∆ β = ( ∆ α, ∆ δ ), the fractional error in the semi-major axis ∆ a / a and the orbital period limited observation time T P s | calculated using Eq. 58 at a frequency of 1 kHz. Althoughwe could expect gravitational wave emission up to ∼ Hz (from the currently measured spin distributions of LMXBsmaxing out at ∼ Hz ), 1 kHz is chosen as an upper boundon the search frequency since the sensitivity of LIGO detec-tors is limited at high frequencies and the amplitudes of thesesystems are not expected to be very strong. Sources withpoorly constrained ( ∆ a / a > .
9) semi-major axis and skyposition ( ∆ β > (cid:48)(cid:48) ) have not been included. The sourcesare listed in order of their bolometric X-ray flux within eachsource group with QPO sources in the top and burst sourcesin the bottom half of the table. The distance is included as areference but is already taken into account in calculation of F X . From these factors alone, Sco X-1 is already the leadingcandidate source. The sky position error ∆ β should be lessthan 100 mas for a source with f = a uncertainties on the scale of 10’s of percent. The final columnlists the orbital period limited observation timespan T P s at afrequency of 1 kHz. Although the spin frequencies of the burstsources are better constrained than QPO sources, the compar- ison of T P s is still made at 1 kHz so that a direct comparisonon the source parameters (rather than search performance) canbe made. This column is included for reference as the orbitalperiod may not be the tightest constraint on the observationtime (c.f. Sec. IV B). It does, however, give an indication ofhow well the orbital period of the source is constrained andspecifically how it a ff ects the search performance. B. Sco X-1
Sco X-1, the first LMXB to be discovered, is also thebrightest extra-solar x-ray source in the sky. The direct rela-tion between gravitational wave strain and x-ray flux given byEq. 7 makes it also the most likely to be a strong gravitationalwave emitter. This, as well as the parameter constraints dis-played in Table I, makes it an ideal candidate for the sidebandsearch. Table II provides a list of Sco X-1 parameters deter-mined from various electromagnetic observations. The tableincludes the parameters required to run the sideband searchtogether with some values used for calculating limits and con-straints on the performance and sensitivity of the search. Thebottom section of the table lists some of the limits and con-straints derived using the above mentioned parameters.Running the standard version of the sideband search re-quires accurate knowledge of the sky position and orbital pe-riod and approximate knowledge of the semi-major axis. Thesky position β = ( α, δ ) listed for Sco X-1 is accurate to within0.3 mas. This error is well within the 100 mas limit definedin Sec. IV A, justifying the assumption of a fixed sky position.The accuracy of measurements of the orbital period requireonly a single sideband template if the observation timespanis within T s < T obss ≈
49 days (for Sco X-1 at 1 kHz). Thesemi-major axis and its measurement error are also requiredfor construction of the sideband template.Estimates of the primary (accreting neutron star) and sec-ondary (donor star) masses, as well measurements of the bolo-metric X-ray flux ( F X ) are required to estimate the indirect,torque balance, gravitational wave strain upper limit h EQ0 usingEq. 7, displayed in the bottom section of the table. The spinfrequency limited observation timespan is also listed here and5
TABLE II: Sco X-1 system parameters required for the sideband search. Directly observable parameters are presented in the top half of thetable. The bottom half, separated by the horizontal line, displays search limits and constraints derived from these.Parameter Symbol Value Units Uncertainty ReferencesRight Ascension α h m s . ± . δ − ◦ (cid:48) .
9” mas ± . µ . − [71]Parallax π β .
36 mas ± .
04 [71]Moment of inertia I kg m [71]Accretion rate ˙ M . × kg s − [71]Bolometric X-ray flux F X × − erg cm − s − [28]Projected semi-major axis light travel time a ± .
18 [72]Orbital Period P ± . ι
44 deg ± ψ
234 deg ± t p ±
400 [72, 75]Strain amplitude (at ν s =
300 Hz) h . × − Eq. 7Spin limited observation timespan T spins
13 days Eq. 55 requires values for the accretion rate ˙ M and moment of inertia I to calculate this value for Sco X-1 using Eq. 55 assuminga spin-wandering timescale τ s = The correspondingvalue of T spins ≈
13 days displayed in the table is more restric-tive than the orbital period limited timespan and is our limitingtime constraint in the search.
VI. STATISTICAL ANALYSIS
Let us first assume that our analysis has yielded no signif-icant candidate signal given a designated significance thresh-old. In this case, with no evidence for detection, we placean upper limit on the possible strength of an underlying sig-nal. In the literature on continuous-wave gravitational sig-nals, it is common to determine these upper limits numeri-cally [53, 78, 79] or semi-analytically [46, 80] using frequen-tist Monte Carlo methods. In these cases simulated signalsare repeatedly added to data over a range of frequencies andrecovered using a localized, computationally cheap, searcharound the point of injection.The sideband algorithm combines signal from many (typi-cally ∼ ) correlated F -statistic frequency bins which mustbe computed over a relatively wide frequency band for eachsimulated signal. Such computations represent a computa-tional cost far in excess of existing methods and are only man-ageable for a small parameter space, e.g. injection studieswhere the signal frequency is known and O (10 ) realizationsare feasible. The computations become daunting for a wide-band search covering more than a few Hz.We choose to optimize the process by calculating upper Assuming the instantaneous accretion rate does not vary more than thex-ray flux, observations of the x-ray variability of Sco X-1 show that theaccretion rate can vary by roughly a factor of two over a timescale τ s = limits within a Bayesian framework. This is an especiallyappealing alternative since the probability density function(PDF) of the C -statistic takes a relatively simple, closed, ana-lytic form. Bayesian upper limits have been computed in time-domain gravitational wave searches targeting known sources(pulsars) [43, 45, 81], and cross-correlation searches for thestochastic background [55, 56, 82]. Comparisons on specificdata sets have shown that Bayesian and frequentist upper lim-its are consistent [42, 45, 83]. A. Bayes Theorem
In the Bayesian framework, the posterior probability den-sity of the hypothesis H given the data D and our backgroundinformation I is defined as p ( θ | D , H , I ) = p ( θ | H , I ) p ( D | θ , H , I ) p ( D | H , I ) , (63)where p ( θ | H , I ) denotes the prior probability distribution ofour model parameters θ given a model H assuming the back-ground information I . The quantity p ( D | θ , H , I ) is the directprobability density (or likelihood function) of the data giventhe parameters, model and background information. The term p ( D | H , I ) is known as the evidence of D given our model andacts as a normalisation constant and does not a ff ect the shapeof the posterior distribution p ( θ | D , H , I ) [84]. The backgroundinformation I (which represents our signal model, assump-tions on Gaussian noise, physicality of parameters etc.) re-mains constant throughout our analysis and will not be men-tioned hereafter. B. Likelihood
When there is no signal in the data, we will say the null hy-pothesis H n , that the data contains any Gaussian noise, is true.6Under these conditions, each C value is drawn from a central χ M distribution. Hence the H n model is parametrized entirelyby M = m +
1, the number of sidebands in the template,where m = ceil [2 π f a ] and depends on the search frequency f and semi-major axis a .The signal hypothesis H GW is true if the data contains Gaus-sian noise plus a signal. The signal is defined by the set ofparameters θ = { h , cos ι, ψ, φ , a , P } . In the case of a sig-nal present in the data, each C -statistic is drawn from a non-central χ M [ λ ( θ )] distribution. The non-centrality parameter λ ( θ ) is defined by the signal parameters θ and is given by λ ( θ ) = ρ m (cid:88) n = − m J n (2 π f a ) | ˜ W ( f n − f k + l ( n ) ∆ f ) | . (64)It represents the total recovered optimal SNR contained withinthe sidebands. The likelihood function (the probability of ourmeasured C value given a parameter set θ ) is then given by p ( C| θ ) =
12 exp (cid:32) −
12 [ C + λ ( θ )] (cid:33) (cid:32) C λ ( θ ) (cid:33) M − I M − (cid:32) (cid:112) C λ ( θ ) (cid:33) . (65)It should be noted that although the quantity M is a functionof the semi-major axis and intrinsic gravitational wave fre-quency, it has been fixed according to the predefined numberof teeth used in the sideband template. It is therefore not afunction of θ . C. Priors
When searching for weak signals, an overly prescriptiveprior is undesirable because it may dominate the posterior.Hence, to be conservative, we adopt a uniform prior on h ≥ h = ff reys prior ∝ / h [81]. The upper limit thus de-rived is consistent with the data, not just a re-iteration of theprior. The same h prior has been adopted in previous searches[42, 53, 78, 85]; the motivation is discussed in more detail in[81].Electromagnetic measurements of the orbital period P andsemi-major axis a are assumed to carry normally distributedrandom errors. Hence we adopt Gaussian priors on the actualvalues P and a . Specifically we take p ( P ) = N ( P , ∆ P ) and p ( a ) = N ( a , ∆ a ), where N ( µ, σ ) denotes a Gaussian (nor-mal) distribution with mean µ given by the electromagneticobservation and standard deviation σ taken as the error in thatobservation.The reference phase φ is automatically maximized overwithin the F stage of the analysis and therefore does not di-rectly a ff ect our (semi-coherent) analysis. The remaining am-plitude parameters serve only to influence the optimal SNR,and therefore also the C . Without prior information fromelectromagnetic observations, we select the least informa-tive (ignorant) physical priors such that p (cos ι ) = / p ( ψ ) = / π on the domains ( − ,
1) and (0 , π ) respectively.Any prior informative measurements (e.g. electromagnetic)on the amplitude parameters can be incorporated into the anal-ysis, and serve to narrow the prior probability distributions. For the Sco X-1 search, we can deduce measurements for cos ι and ψ from observations if we assume the rotation axes of theneutron star and accretion disk are aligned. This implies theneutron star inclination ι is equal to the orbital inclination.We can then set ι = ◦ ± ◦ from the inclination of the orbitalplane suggested from observations of the radio components ofSco X-1 [74]. The same observations measure a position angleof these radio jets of 54 ± ◦ . Under the alignment assump-tion, the position angle is directly related to the gravitationalwave polarization angle ψ , but with a phase shift of 180 ◦ , i.e. ψ = ± ◦ . The above assumes the usual mass-quadrupoleemission; for current-quadrupole emission from r -modes theresults are the same with ψ → ψ + ◦ [63]. D. Posteriors
The probability density function (PDF) on our search pa-rameters given a single C -statistic value is p ( θ |C ) ∝ p ( C| θ ) p ( θ ) , (66)and assuming that the prior PDFs on our parameters are inde-pendent, we can express the posterior PDF as p ( h , cos ι, ψ, P , a |C ) ∝ p ( C| θ ) p ( h ) p (cos ι ) p ( ψ ) p ( P ) p ( a ) . (67)To perform inference on the gravitational wave strain h , wecan marginalize this joint distribution over the other parame-ters leaving us with p ( h |C ) ∝ ∞ (cid:90) −∞ da ∞ (cid:90) ∞ dP π (cid:90) d ψ (cid:90) − d cos ι p ( C| θ ) N ( P , ∆ P ) N ( a , ∆ a ) , (68)where the flat priors on h , cos ι and ψ are absorbed intothe proportionality. Note that the amplitude parameters actthrough the non-centrality parameter λ ( θ ) (Eq. 64) via the op-timal SNR term (Eq. 22), in the likelihood. The orbital param-eters a , P dictate the fraction of recovered SNR based on themismatch in the predicted quantity and location of frequency-modulated sidebands (Eq. 64). E. Detection criteria and upper limits
To determine whether or not a signal is present in the data,we compute a threshold value of the C -statistic such that theprobability of achieving such a value or greater due to noisealone is P a , the false alarm probability. For a single measure-ment of the C -statistic this threshold is computed via P a = ∞ (cid:90) C ∗ p ( C| h = = − P (2 M , C ∗ / , (69)7where the likelihood on C in the noise only case becomesa central χ distribution and P ( k / , x /
2) is the regularizedGamma function with k degrees of freedom (the cumulativedistribution function of a central χ k distribution) defined at x .In the case of N measurements of the C -statistic, assumingstatistically independent trials, the false alarm probability isgiven by P a | N = − (1 − P a ) N = − (cid:2) P (2 M , C ∗ / (cid:3) N . (70)The corresponding threshold C ∗ N such that the probability thatone or more of these values exceeds that threshold is obtainedby solving P (2 M , C ∗ N / = (cid:0) − P a | N (cid:1) / N . (71)This solution is obtained numerically but can be representednotationally by C ∗ N = P − (cid:16) M , (cid:2) − P a | N (cid:3) / N (cid:17) , (72)where P − represents the inverse function of P .In practice the C -statistic values will not be statistically in-dependent as assumed above. The level of independence be-tween adjacent frequency bins will be reduced (i.e. valueswill be become increasingly correlated) as the frequency res-olution of the C -statistic is made finer. Additionally, due to thecomb structure of the signal and template we find that resultsat frequencies separated by an integer number of frequency-modulated sideband spacings j / P Hz for j < m are highly cor-related. This is due to the fact that these results will have beenconstructed from sums of F -statistic values containing manycommon values. This latter e ff ect is dominant over the formerand as an approximation it can be assumed that within thefrequency span of a single comb template there are rT / P in-dependent C -statistic results. The number of templates spansper unit search frequency is ∼ P / rM which leaves us with ∼ T / M independent C -statistic values per unit Hz. This is areduction by a factor of M in the number of statistically inde-pendent results expected.In the event of there being no candidate C -statistic values,the search allows us to compute upper-limits on the amplitudeof gravitational waves from our target source. We define theupper limit on the wave strain h as the value h UL that boundsthe fraction UL of the area of the marginalized posterior distri-bution p ( h |C ). This value is obtained numerically by solvingUL = h UL (cid:90) p ( h |C ) dh . (73)We note that this procedure allows us to compute an upper-limit for each C -statistic value output from a search. The stan-dard practice in continuous gravitational wave data analysis This comes from the number of bins in between each sideband, given bythe sideband separation 1 / P divided by the bin size d f = ( rT s ) − (Eq. 57). −26 −25 −24 −23 −22 frequency (Hz) s t r a i n a m p li t ude h S2S4S5Adv t o r que − ba l an c e li m i t FIG. 5: Sensitivity estimate for a 10 day standard, approximatedemodulated and approximate demodulated with known priors side-band search (fine, medium and bold solid curves respectively) usingLIGO (H1L1) S5 data (upper, purple group), and using the 3-detector(H1L1V) advanced LIGO configuration (lower, red group). Alsoshown are results from the the previous coherent search for Sco X-1in S2 data (solid black dashes) [53] and the maximum upper limitsfor each Hz band of the directed stochastic (radiometer) search in S4and S5 data (light and dark blue dashed curves, respectively) [55, 56].The theoretical torque-balance gravitational wave strain upper limit( h EQ0 from Eq. 7) for Sco X-1 is indicated by the thick gray straightline. is to perform a frequentist upper-limit using computationallyexpensive Monte-Carlo simulations involving repeated sig-nal injections. The results of these injections are then com-pared to loudest detection statistic recovered from the actualsearch [53]. In our approach, by virtue of the fact that weare able to compute upper-limits very e ffi ciently for each C -statistic value and the upper-limit value is a monotonic func-tion of C we naturally also include the worst case (loudestevent) result. The di ff erence in the upper-limits obtained fromboth strategies then becomes an issue of Bayesian versus Fre-quentist interpretation. However, as shown in [83], in the limitof large SNR these upper-limit results become indistinguish-able. When searching wide parameter spaces with large num-bers of templates, as is the case for the sideband search, themost likely largest detection statistic value will be consistentwith large SNR. VII. SENSITIVITY
The sensitivity of a future search can be predicted in a va-riety of ways. We choose to estimate the expected gravita-tional wave strain upper-limits for Initial LIGO data in orderto compare against previous results. We also compare this tothe expected sensitivity of the search with Advanced LIGO.If the search is conducted such that the frequency space issplit into small sub-bands, the sensitivity can be estimated by8computing upper limits on the expected maximum from eachof the sub-bands in Gaussian noise. This is equivalent to as-signing a false alarm probability P a =
50% for N = T s / M tri-als for each, say one Hz frequency sub-band, and using Eq. 72as the expected C -statistic. We can then calculate the posteriordistribution of h from Eqs. 65 and Eq 68.Figure 5 shows the sensitivity estimate of the 90% upperlimit (UL = ff erent modes:standard (described in Section III, represented by the thinsolid curves), binary demodulated (described in Section III H,represented by the medium solid curves), and binary demod-ulated with known priors on cos ι and ψ (described in SectionVI C, represented by the bold solid curves). It compares thesensitivity of the search in two-detector (H1L1) LIGO S5 data(upper, purple group) and three-detector (H1L1V) AdvancedLIGO data (red group) with previous searches for Sco X-1in LIGO S2 (black dashes) [53], S4 and S5 data (light anddark blue dashed curves, respectively) [55, 56]. The h rms up-per limit quoted in the latter two (radiometer) searches is op-timized for the special case of a circularly polarized signaland hence less conservative than the angle averaged h quotedin [53] and commonly used when quoting upper limits forcontinuous gravitational wave searches. Converting detector-strain rms upper limits h rms to source-strain amplitude upperlimit h requires h ∼ . h rms (see [86]). The di ff erent con-fidence on the coherent S2 analysis and S4 ans S5 radiometeranalyses (90 and 95% respectively) also complexify any di-rect comparisons. The theoretical indirect wave strain limit h EQ0 for gravitational waves from LMXBs represented by thethick gray line comes from Eq. 7.The sensitivity curves in Fig. 5 show that the standard side-band search should improve current upper limits on gravita-tional waves from Sco X-1, even though it is limited to only10 days of consecutive data. Running a demodulated searchwith known cos ι and ψ comes close to setting constraints onthe indirect (torque-balance) upper limits in the advanced de-tector era. VIII. DISCUSSION
We have described the sideband algorithm and shown thatit provides a computationally e ffi cient method to search forgravitational waves from sources in binary systems. It re-quires accurate knowledge of the sky position of the sourceand the orbital period of the binary, and less accurate knowl-edge on the semi-major axis. E ff ects of spin wandering can beignored over a short enough coherent integration time.The tolerance on the errors of relevant search parameterswas computed, defining the range over which they can be as-sumed constant (Section IV). In light of these limits, electro-magnetic observations suggest several candidates (Section V).Of these sources, Sco X-1 is identified as the strongest candi- date based on the gravitational wave strain recovered from thetorque-balance argument (Eq. 7). In future, the search canalso be directed at several other of the suitable LMXB candi-dates presented in Section V.A Bayesian upper limit strategy was presented in SectionVI, rather than the frequentist methods commonly employedin frequency-based (LIGO) searches. Knowing the likelihoodfunction in closed analytic form makes the Bayesian approachcomputationally more feasible than Monte Carlo simulations(see VI B). Knowing the gravitational wave polarization angleand inclination leads to additional sensitivity improvementsusing this framework (see VI C).The sensitivity of the search, described in Section VII, isestimated by performing the Bayesian analysis on the designcurves of the S5 and Advanced LIGO noise floors. The sensi-tivity of a 10 day limited Sco X-1 directed sideband searchcompared to previous LIGO searches is shown in Fig. 5.It shows that measurements of an orbital reference time andphase (the time and argument of periapse) can be employedto improve search sensitivity by a factor of 1 . .
5. In its mostsensitive configuration (approximated binary demodulated as-suming known cos ι and ψ in the Advanced detector era), thesideband search brings us closer to testing the theoretical in-direct torque-balance limit.The studies presented here assume pure Gaussian noise.The performance for realistic LIGO-like noise will be pre-sented elsewhere, in a report on the results from the Sco X-1directed search performed on LIGO (S5) data. The searchcould also look forward to running on the next-generation Ad-vanced LIGO data. ACKNOWLEDGEMENTS
This work has benefited from many useful discussions andhelpful comments. In particular we would like to thank KarlWette, Reinhard Prix, Evan Goetz, Holger Pletsch, MiroslavShaltev, Badri Krishnan, Paul Lasky, Mark Bennett, SterlPhinney and the LIGO Scientific Collaboration continuouswaves working group. LS is indebted to Bruce Allen for hissupport and gratefully acknowledges the Max-Planck Soci-ety (Albert-Einstein Institut) for support and hospitality dur-ing periods when significant progress on this work was made.LS and CM acknowledge the Gravitational wave group atCardi ff University for their support. LS was supported by anAustralian Postgraduate Award and a Melbourne UniversityOverseas Research Experience Scholarship. AM and BJO ac-knowledge ARC grant DP-110103347. BJO was supportedby NSF Grants PHY-0855589 and PHY-1206027. [1] L. Bildsten, Astrophys. J. Lett. , L89 (1998), URL http://stacks.iop.org/1538-4357/501/i=1/a=L89 . [2] G. Ushomirsky, C. Cutler, and L. Bildsten, Mon. Not. R. As- tron. Soc. , 902 (2000), URL http://dx.doi.org/10.1046/j.1365-8711.2000.03938.x .[3] B. J. Owen, Phys. Rev. Lett. , 211101 (2005), arXiv:astro-ph / , 231101 (2007), URL http://link.aps.org/doi/10.1103/PhysRevLett.99.231101 .[5] L.-M. Lin, Phys. Rev. D76 , 081502 (2007), 0708.2965.[6] N. K. Johnson-McDaniel and B. J. Owen, Phys. Rev. D ,044004 (2013), 1208.5227.[7] S. Bonazzola and E. Gourgoulhon, Astron. Astrophys. , 675(1996), astro-ph / , 084025 (2002), arXiv:gr-qc / , 1044(2005), URL http://stacks.iop.org/0004-637X/623/i=2/a=1044 .[10] M. Vigelius and A. Melatos, Mon. Not. R. Astron. Soc. ,1294 (2008), 0802.3238.[11] B. Haskell, L. Samuelsson, K. Glampedakis, and N. Andersson,Mon. Not. R. Astron. Soc. , 531 (2008), 0705.1780.[12] A. Mastrano and A. Melatos, Mon. Not. R. Astron. Soc. ,760 (2012), 1112.1542.[13] D. I. Jones and N. Andersson, Mon. Not. R. Astron. Soc. ,203 (2002), arXiv:gr-qc / , 2325 (2012),1107.3503.[15] B. J. Owen, L. Lindblom, C. Cutler, B. F. Schutz, A. Vecchio,et al., Phys. Rev. D58 , 084020 (1998), gr-qc / , 307 (1999), arXiv:astro-ph / D76 , 064019 (2007), 0704.0799.[18] B. Haskell, N. Degenaar, and W. C. G. Ho, Mon. Not. R. As-tron. Soc. , 93 (2012), 1201.2101.[19] R. Bondarescu and I. Wasserman (2013), 1305.2335.[20] M. Nayyar and B. J. Owen, Phys. Rev. D , 084001 (2006),URL http://link.aps.org/doi/10.1103/PhysRevD.73.084001 .[21] A. Melatos, Advances in Space Research , 1472 (2007),ISSN 0273-1177, URL .[22] C. A. van Eysden and A. Melatos, Class. Quant. Grav. , 225020 (2008), URL http://stacks.iop.org/0264-9381/25/i=22/a=225020 .[23] M. Vigelius and A. Melatos, Mon. Not. R. Astron. Soc. ,1972 (2009), 0902.4264.[24] N. Andersson, V. Ferrari, D. I. Jones, K. D. Kokkotas, B. Kr-ishnan, J. S. Read, L. Rezzolla, and B. Zink, General Relativityand Gravitation , 409 (2011), 0912.0384.[25] J. Papaloizou and J. E. Pringle, Mon. Not. R. Astron. Soc. ,501 (1978).[26] R. V. Wagoner, Astrophys. J. , 345 (1984).[27] F. Verbunt, Annu. Rev. Astron. Astrophys. , 93 (1993).[28] A. L. Watts, B. Krishnan, L. Bildsten, and B. F. Schutz, Mon.Not. R. Astron. Soc. , 839 (2008).[29] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astrophys. J. , 823 (1994).[30] G. Ushomirsky, L. Bildsten, and C. Cutler, in American Insti-tute of Physics Conference Series , edited by S. Meshkov (2000),vol. 523 of
American Institute of Physics Conference Series , pp.65–74, arXiv:astro-ph / , 42 (2003), arXiv:astro-ph / , 198 (2007),arXiv:astro-ph / , L148 (2010), 0910.5546.[34] B. P. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen,G. Allen, R. S. Amin, S. B. Anderson, W. G. Anderson, M. A.Arain, et al., Reports on Progress in Physics , 076901 (2009),0711.3041.[35] J. Abadie, B. P. Abbott, R. Abbott, M. Abernathy, C. Adams,R. Adhikari, P. Ajith, B. Allen, G. Allen, E. Amador Ceron,et al., Nuclear Instruments and Methods in Physics Research A , 223 (2010), 1007.3973.[36] URL .[37] URL .[38] T. Accadia et al. (VIRGO Collaboration), JINST , P03012(2012).[39] H. Grote (LIGO Scientific Collaboration), Class. Quant. Grav. , 084003 (2010).[40] URL .[41] P. R. Brady, T. Creighton, C. Cutler, and B. F. Schutz, Phys.Rev. D , 2101 (1998), URL http://link.aps.org/doi/10.1103/PhysRevD.57.2101 .[42] B. Abbott, R. Abbott, R. Adhikari, A. Ageev, B. Allen,R. Amin, S. B. Anderson, W. G. Anderson, M. Araya, H. Ar-mandula, et al. (The LIGO Scientific Collaboration), Phys. Rev.D , 082004 (2004), arXiv:gr-qc / , 671 (2010), 0909.3583.[44] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen,R. Amin, S. B. Anderson, W. G. Anderson, M. A. Arain, et al.(The LIGO Scientific Collaboration), Astrophys. J. Lett. ,L45 (2008), 0805.4758.[45] J. Abadie, B. P. Abbott, R. Abbott, M. Abernathy, T. Accadia,F. Acernese, C. Adams, R. Adhikari, C. A ff eldt, B. Allen, et al.(The LIGO Scientific Collaboration and the Virgo Collabora-tion), Astrophys. J. , 93 (2011), 1104.2712.[46] J. Abadie, B. P. Abbott, R. Abbott, M. Abernathy, C. Adams,R. Adhikari, P. Ajith, B. Allen, G. Allen, E. Amador Ceron,et al. (LIGO Scientific Collaboration), Astrophys. J. , 1504(2010), 1006.2535.[47] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen,R. Amin, S. B. Anderson, W. G. Anderson, M. Arain, et al.(LIGO Scientific Collaboration), Phys. Rev. D , 022001(2008), 0708.3818.[48] B. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen,R. Amin, D. P. Anderson, S. B. Anderson, W. G. Anderson,et al. (The LIGO Scientific Collaboration), Phys. Rev. D ,022001 (2009), 0804.1747.[49] B. P. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen,G. Allen, R. S. Amin, S. B. Anderson, W. G. Anderson, M. A.Arain, et al. (The LIGO Scientific Collaboration), Phys. Rev.Lett. , 111102 (2009), 0810.0283.[50] B. P. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen,G. Allen, R. S. Amin, S. B. Anderson, W. G. Anderson, M. A.Arain, et al. (The LIGO Scientific Collaboration and the VirgoCollaboration), Phys. Rev. D , 042003 (2009), 0905.1705.[51] J. Abadie, B. P. Abbott, R. Abbott, T. D. Abbott, M. Abernathy,T. Accadia, F. Acernese, C. Adams, R. Adhikari, C. A ff eldt,et al. (The LIGO Scientific Collaboration), Phys. Rev. D ,022001 (2012), 1110.0208.[52] E. Goetz and K. Riles, Class. Quant. Grav. , 215006 (2011), URL http://stacks.iop.org/0264-9381/28/i=21/a=215006 .[53] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen,R. Amin, S. B. Anderson, W. G. Anderson, M. Arain, et al.(LIGO Scientific Collaboration), Phys. Rev. D , 082001(2007).[54] P. Jaranowski, A. Kr´olak, and B. F. Schutz, Phys. Rev. D ,063001 (1998).[55] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen,R. Amin, S. B. Anderson, W. G. Anderson, M. Arain, et al.(LIGO Scientific Collaboration), Phys. Rev. D , 082003(2007).[56] J. Abadie, B. P. Abbott, R. Abbott, M. Abernathy, T. Accadia,F. Acernese, C. Adams, R. Adhikari, P. Ajith, B. Allen, et al.(LIGO Scientific Collaboration, Virgo Collaboration), Phys.Rev. Lett. , 271102 (2011), 1109.1809.[57] P. R. Brady and T. Creighton, Phys. Rev. D , 082001 (2000),arXiv:gr-qc / , 084010 (2012),1201.4321.[59] C. Messenger and G. Woan, Class. Quant. Grav. , 469 (2007),arXiv:gr-qc / , 911 (2003), URL http://stacks.iop.org/0004-637X/589/i=2/a=911 .[61] T. M. Tauris and E. P. J. van den Heuvel (Cambridge UniversityPress, 2006), chap. 16, URL http://dx.doi.org/10.1017/CBO9780511536281 .[62] J. van Paradijs and N. White, ”Astrophys. J. Lett.” , L33 + (1995).[63] B. J. Owen, Phys. Rev. D82 , 104002 (2010), 1006.1994.[64] J. H. Taylor and J. M. Weisberg, Astrophys. J. , 434 (1989).[65] R. Prix, Phys. Rev. D , 023004 (2007).[66] K.-H. Lo, F. K. Lamb, S. Boutloukos, and M. C. Miller,in AAS / High Energy Astrophysics Division (2011), vol. 12 of
AAS / High Energy Astrophysics Division , p. 44.05.[67] A. L. Watts, ”Ann. Rev. Astron. Astrophys.” , 609 (2012),1203.2065.[68] M. van der Klis, Advances in Space Research , 925(1998), ISSN 0273-1177, x-Ray Timing and CosmicGamma Ray Bursts, URL .[69] C. M. Zhang, Y. Y. Pan, J. Wang, A. Taani, and Y. H. Zhao, Int.J. Mod. Phys. Conference Series , 414 (2012).[70] M. van der Klis, J. H. Swank, W. Zhang, K. Jahoda, E. H.Morgan, W. H. G. Lewin, B. Vaughan, and J. van Paradijs, As-trophys. J. Lett. , L1 (1996), URL http://stacks.iop. org/1538-4357/469/i=1/a=L1 .[71] C. F. Bradshaw, E. B. Fomalont, and B. J. Geldzahler, Astro-phys. J. Lett. , L121 (1999), URL http://stacks.iop.org/1538-4357/512/i=2/a=L121 .[72] D. Steeghs and J. Casares, The Astrophysical Journal ,273 (2002), URL http://stacks.iop.org/0004-637X/568/i=1/a=273 .[73] D. Galloway, S. Premachandra, D. Steeghs, J. Casares, andC. Messenger, GWPAW Hannover.[74] E. B. Fomalont, B. J. Geldzahler, and C. F. Bradshaw, The As-trophysical Journal , 283 (2001), URL http://stacks.iop.org/0004-637X/558/i=1/a=283 .[75] C. J. Messenger, Ph.D. thesis, The University of Birmingham(2006).[76] H. V. Bradt, G. Moore, L. L. E. Braes, G. K. Miley, W. For-man, E. Kellogg, J. E. Hesser, W. E. Kunkel, W. A. Hiltner,R. Hjellming, et al., Astrophys. J. , 443 (1975).[77] P. Hertz, B. Vaughan, K. S. Wood, J. P. Norris, K. Mitsuda, P. F.Michelson, and T. Dotani, Astrophys. J. , 201 (1992).[78] B. Abbott, R. Abbott, R. Adhikari, J. Agresti, P. Ajith, B. Allen,R. Amin, S. B. Anderson, W. G. Anderson, M. Arain, et al.(The LIGO Scientific Collaboration), Phys. Rev. D , 022001(2008), 0708.3818.[79] J. Aasi, J. Abadie, B. P. Abbott, R. Abbott, T. D. Abbott,M. Abernathy, T. Accadia, F. Acernese, C. Adams, T. Adams,et al. (The LIGO Scientific Collaboration and the Virgo Collab-oration), ArXiv e-prints (2012), 1207.7176.[80] K. Wette, Phys. Rev. D , 042003 (2012), URL http://link.aps.org/doi/10.1103/PhysRevD.85.042003 .[81] R. J. Dupuis and G. Woan, Phys. Rev. D , 102002 (2005),arXiv:gr-qc / , 022001 (2007), arXiv:gr-qc / Bayesian Spectrum Analysis and ParameterEstimation , vol. 48 of