Search for gravitational waves from Scorpius X-1 in the first Advanced LIGO observing run with a hidden Markov model

TL;DR

This study conducted a semi-coherent search for continuous gravitational waves from Scorpius X-1 using LIGO O1 data, employing a hidden Markov model to account for spin frequency wandering, but found no evidence of such waves.

Contribution

It introduces a novel semi-coherent search method combining a matched filter with a hidden Markov model for detecting gravitational waves from low-mass X-ray binaries.

Findings

No gravitational waves detected in the 60-650 Hz range.

Upper limits on strain amplitude were established at 106 Hz.

Upper limits are within an order of magnitude of the theoretical torque-balance limit.

Abstract

Results are presented from a semi-coherent search for continuous gravitational waves from the brightest low-mass X-ray binary, Scorpius X-1, using data collected during the first Advanced LIGO observing run (O1). The search combines a frequency domain matched filter (Bessel-weighted $\mathcal{F}$-statistic) with a hidden Markov model to track wandering of the neutron star spin frequency. No evidence of gravitational waves is found in the frequency range 60-650 Hz. Frequentist 95% confidence strain upper limits, $h_0^{95\%} = 4.0\times10^{-25}$, $8.3\times10^{-25}$, and $3.0\times10^{-25}$ for electromagnetically restricted source orientation, unknown polarization, and circular polarization, respectively, are reported at 106 Hz. They are $\leq 10$ times higher than the theoretical torque-balance limit at 106 Hz.

Figures

Tables

Table 1. Table 1: Electromagnetically observed parameters (top half) and search parameters (bottom half) for Sco X-1. The uncertainties are at the 1 ​ σ 1 𝜎 1\sigma confidence level.

Observed parameter	Symbol	Value	Reference
Right ascension	$\alpha$	16h 19m 55.0850s	Bradshaw et al. (1999)
Declination	$\delta$	$-{15}^{\circ }{38}^{\prime }{24.9}^{\prime \prime }$	Bradshaw et al. (1999)
X-ray flux	$F_X$	$4\times {10}^{-7}\mathrm{erg}{\mathrm{cm}}^{-2}{s}^{-1}$	Watts et al. (2008)
Orbital period	$P$	$68023.70496\pm 0.0432$ s	Galloway et al. (2014)
Projected semi-major axis	$a_0$	$1.44\pm 0.18$ s	Steeghs et al. (2002)
Polarization angle	$\psi$	$234\pm 3^{\circ }$	Fomalont et al. (2001)
Inclination angle	$\iota$	$44\pm 6^{\circ }$	Fomalont et al. (2001)
Search parameter	Symbol	Search range	Resolution
Frequency	$f_0$	$60-650$ Hz	$5.787037\times {10}^{-7}$ Hz
Projected semi-major axis	$a_0$	$0.361-3.249$ s	0.01805 s

Table 2. Table 2: Number of candidates surviving each veto.

Veto	Number
First pass	180
After line veto	44
After single IFO veto	6
After half $T_{\mathrm{obs}}$ veto	2
After longer $T_{\mathrm{drift}}$ veto	0

Table 3. Table 3: Actions to be taken for survivors of the known line veto in Section IV.1.1 according to the score S 𝑆 S and the estimated frequency f 0 subscript 𝑓 0 f_{0} from each single detector. S ∪ subscript 𝑆 S_{\cup} and f 0 ∪ subscript subscript 𝑓 0 {f_{0}}_{\cup} stand for the score and estimated frequency yielded by the original search combining two detectors.

Category	Score in one detector $S$	Estimated frequency in one detector $f_0$	Action
A	$S\ge S_{\cup }$ in one detector but in the other	$f_0\approx f_0^{}{}_{\cup }{}^{}$ where $S\ge S_{\cup }$	Veto
B	$S\ge S_{\cup }$ in one detector but in the other	$f_0\ne f_0^{}{}_{\cup }{}^{}$ where $S\ge S_{\cup }$	Keep
C	$S\ge S_{\mathrm{th}}$ in both detectors		Keep
D	in both detectors		Keep

Table 4. Table 4: Candidates surviving both the known line veto and the single interferometer veto. The table lists the sub-band where the candidate is found (column 1), the estimated frequency f 0 subscript 𝑓 0 f_{0} quoted as the arithmetic mean of the minimum and the maximum frequencies ( f min subscript 𝑓 min f_{\rm min} and f max subscript 𝑓 max f_{\rm max} ) in the optimal HMM path (column 2), the number of frequency bins ( Δ ​ f drift Δ subscript 𝑓 drift \Delta f_{\rm drift} ) between f max subscript 𝑓 max f_{\rm max} and f min subscript 𝑓 min f_{\rm min} (column 3), the estimated a 0 subscript 𝑎 0 a_{0} (column 4), the original score S ∪ subscript 𝑆 S_{\cup} yielded by searching the whole data set (column 5), and the scores from searching the first and second half of the data separately (column 6 and 7). The resolutions of f 0 subscript 𝑓 0 f_{0} and a 0 subscript 𝑎 0 a_{0} are 5.787037 × 10 − 7 5.787037 superscript 10 7 5.787037\times 10^{-7}  Hz and 0.01805 s, respectively. The candidates marked with an asterisk survive the manual veto in Section IV.1.3 and require further follow up.

Sub-band (Hz)	$f_0$ (Hz)	$f_{\max }-f_{\min }$	$a_0$ (s)	$S_{\cup }$	$S_{1st\mathrm{half}}$	$S_{2nd\mathrm{half}}$
		($\Delta f_{\mathrm{drift}}$)
395–396	395.8561536	3	0.81	8.05153	6.55545	9.13679
449–450*	449.8116935	3	2.73	7.38701	6.46122	6.50190
459–460	459.5557459	6	0.38	12.76130	14.61070	6.30887
534–535	534.3625717	4	0.36	20.18630	20.53770	6.97788
548–549	548.9457104	7	0.42	16.68650	18.39020	6.46258
593–594*	593.7716675	1	2.98	7.40397	6.17976	5.88553

Table 5. Table 5: Subsequent actions to be taken for survivors of the vetoes in Section IV.1.1 – IV.1.3 according to the amount of spin wandering and S 𝑆 S -versus- T drift subscript 𝑇 drift T_{\rm drift} trend observed by the HMM.

	Higher $S$ with longer $T_{\mathrm{drift}}$	Lower $S$ with longer $T_{\mathrm{drift}}$
Low spin wandering	Follow up with more sensitive method	Veto
High spin wandering	Unlikely to happen	Follow up with more sensitive method
		guided by observed Viterbi path

Table 6. Table 6: Final-step follow-up with longer T drift = 20 subscript 𝑇 drift 20 T_{\rm drift}=20  d in two 1-Hz sub-bands containing the survivors from Section IV.1.3 . The top and bottom halves of the table correspond to T drift = 10 subscript 𝑇 drift 10 T_{\rm drift}=10  d and 20 d, respectively. The estimated f 0 subscript 𝑓 0 f_{0} is quoted as the arithmetic mean of min f 0 ​ ( t ) subscript 𝑓 0 𝑡 f_{0}(t) and max f 0 ​ ( t ) subscript 𝑓 0 𝑡 f_{0}(t) for the optimal Viterbi path. The follow-up score S 𝑆 S with T drift = 20 subscript 𝑇 drift 20 T_{\rm drift}=20  d is always below S th = 7.34 subscript 𝑆 th 7.34 S_{\rm th}=7.34 and lower than the original score. The resolutions of a 0 subscript 𝑎 0 a_{0} and f 0 subscript 𝑓 0 f_{0} are 0.01805 s and 5.787037 × 10 − 7 5.787037 superscript 10 7 5.787037\times 10^{-7}  Hz, respectively for both T drift = 10 subscript 𝑇 drift 10 T_{\rm drift}=10  d and 20 d.

$T_{\mathrm{drift}}$	Quantity	449–450 Hz	593–594 Hz
10 d	$S$	7.38701	7.40397
	$f_0$ (Hz)	449.8116936	593.7716675
	$a_0$ (s)	2.73	2.98
20 d	$S$	6.93366	6.93900
	$f_0$ (Hz)	449.7891863	593.6174193
	$a_0$ (s)	1.70	1.79

Table 7. Table 7: Veto safety verification with synthetic signals. The table lists the data used for the injections (column 1), the injected signal parameters (column 2–5), the original score S ∪ subscript 𝑆 S_{\cup} yielded by searching the whole data set with two interferometers and T obs = 10 subscript 𝑇 obs 10 T_{\rm obs}=10  d (column 6), the scores from searching H1 and L1 separately (column 7 and 8), the scores from searching the first and second half of the data separately (column 9 and 10), and the score with T obs = 20 subscript 𝑇 obs 20 T_{\rm obs}=20  d (column 11). A score is marked with an asterisk if it is above threshold, but the estimated frequency differs significantly from f 0 inj subscript subscript 𝑓 0 inj {f_{0}}_{\rm inj} (i.e. wrong path returned). These twelve injections survive the four veto stages described in Section IV.1.1 – IV.1.4 .

Data	$f_0^{}{}_{\mathrm{inj}}{}^{}$ (Hz)	$a_0^{}{}_{\mathrm{inj}}{}^{}$ (s)	$h_0^{}{}_{\mathrm{inj}}{}^{}$ (${10}^{-25}$)	$\cos {\iota }_{\mathrm{inj}}$	$S_{\cup }$	$S_{\mathrm{H1}}$	$S_{\mathrm{L1}}$	$S_{1st\mathrm{half}}$	$S_{2nd\mathrm{half}}$	$S_{20d}$
S5	64.5774908	0.81	9.58	$-0.5936$	9.12097		7.42935*		7.67254	11.7985
S5	102.2907797	2.47	9.81	$-0.7988$	20.81940	16.17190	12.00540	20.63850	17.38740	25.81390
S5	202.8863982	2.34	11.25	$-0.9205$	15.80950	18.54850	17.46390	19.15680	18.9102	21.4415
S5	254.6697757	3.03	14.55	0.0375	12.50180		9.27111	10.8953	7.74954	15.0849
O1	97.2345635	2.15	4.50	0.71935	9.76216		7.53014*	7.29089	8.91108	9.98727
O1	132.1234568	0.70	4.80	$-0.68154$	16.86500	8.90286	8.63928	13.29010	13.30940	19.54900
O1	185.8094752	1.11	9.90	0.37952	19.05450	14.44080	12.70840	18.07160	17.95120	20.34430
O1	233.9125689	0.46	4.60	0.70917	16.71220		9.18889	12.25070	13.15180	18.02530
O1	345.3456700	1.45	7.00	0.71567	14.09400		9.15852	10.10120	12.83390	14.72410
O1	454.4563891	3.20	7.00	$-0.86725$	9.03162	7.54074*				9.06928
O1	525.7096896	2.81	12.90	0.66578	11.55910	7.83362	8.90156	11.35370	10.04660	13.26430
O1	635.6679700	1.98	10.00	0.72650	10.64010			8.91769	9.13239	11.56240

Equations

$$\begin{array}{rl}P(Q|O)=& L_{o(t_{N_T})q(t_{N_T})}{A}_{q(t_{N_T})q(t_{N_T-1})}\cdots L_{o(t_1)q(t_1)}\\ & \times A_{q(t_1)q(t_0)}{\Pi }_{q(t_0)},\end{array}$$

$$S=\frac{\ln {\delta }_{q^{\ast }}(t_{N_T})-{\mu }_{\ln \delta }(t_{N_T})}{{\sigma }_{\ln \delta }(t_{N_T})}$$

$${\mu }_{\ln \delta }(t_{N_T})=N_Q^{-1}\sum _{i=1}^{N_Q}\ln {\delta }_{q_i}(t_{N_T})$$

$${\sigma }_{\ln \delta }(t_{N_T}{)}^2=N_Q^{-1}\sum _{i=1}^{N_Q}[\ln {\delta }_{q_i}(t_{N_T})-{\mu }_{\ln \delta }(t_{N_T}){]}^2,$$

$$h(t)\propto \sum _{n=-\infty }^{\infty }{J}_n(2\pi f_0{a}_0)\cos [2\pi (f_0+n/P)t],$$

$$G(f)=\mathcal{F}(f)\otimes B(f),$$

$$B(f)=\sum _{n=-(M-1)/2}^{(M-1)/2}[J_n(2\pi f{a}_0){]}^2\delta (f-n/P).$$

$$h_0^{\mathrm{eff}}=h_0{2}^{-1/2}\{[(1+{\cos }^2\iota )/2{]}^2+{\cos }^2\iota {\}}^{1/2},$$

$$h_0^{95\%}\propto S_h^{1/2}{f}_0^{1/4}.$$

$$h_0^{\mathrm{eq}}=$$

$$R_A$$

$$A_{q_j{q}_i}=\mathrm{Pr}[q(t_{n+1})=q_j|q(t_n)=q_i]$$

$$L_{o_j{q}_i}=\mathrm{Pr}[o(t_n)=o_j|q(t_n)=q_i].$$

$$\begin{array}{rl}P(Q|O)=& L_{o(t_{N_T})q(t_{N_T})}{A}_{q(t_{N_T})q(t_{N_T-1})}\cdots L_{o(t_1)q(t_1)}\\ & \times A_{q(t_1)q(t_0)}{\Pi }_{q(t_0)},\end{array}$$

$${\Pi }_{q_i}=\mathrm{Pr}[q(t_0)=q_i]$$

$$\left|{\int }_t^{t+T_{\mathrm{drift}}}d{t}^{\prime }\dot{f_0}(t^{\prime })\right|<\Delta f_{\mathrm{drift}}$$

$$A_{q_{i+1}{q}_i}=A_{q_i{q}_i}=A_{q_{i-1}{q}_i}=\frac{1}{3},$$

$$L_{o(t)q_i}\propto \exp [G({f_0}_i)],$$

$${\Pi }_{q_i}=N_Q^{-1}.$$

$${\delta }_{q_i}(t_k)=L_{o(t_k)q_i}\underset{1\le j\le N_Q}{\max }[A_{q_i{q}_j}{\delta }_{q_j}(t_{k-1})].$$

$${\Phi }_{q_i}(t_k)=\underset{1\le j\le N_Q}{\mathrm{argmax}}[A_{q_i{q}_j}{\delta }_{q_j}(t_{k-1})],$$

$$q^{\ast }(t_k)={\Phi }_{q^{\ast }(t_{k+1})}(t_{k+1})$$

Full text

PACS numbers

95.85.Sz, 97.60.Jd

††preprint: APS/123-QED

I Introduction

Rotating neutron stars are a possible source of persistent, periodic gravitational radiation. The signal is expected at specific multiples of the neutron star spin frequency $f_{\star}$ Riles (2013). Astrophysical models suggest that the radiation may be emitted at levels detectable by ground-based, long-baseline interferometers such as the Laser Interferometer Gravitational Wave Observatory (LIGO) and the Virgo detector Harry et al. (2010); Aasi et al. (2015a); Acernese et al. (2015); Andersson et al. (2011); Riles (2013). A time-varying quadrupole moment can result from thermal Ushomirsky et al. (2000); Johnson-McDaniel and Owen (2013) or magnetic Cutler (2002); Mastrano et al. (2011); Lasky and Melatos (2013) gradients, r-modes Owen et al. (1998); Heyl (2002); Arras et al. (2003); Bondarescu et al. (2009), or nonaxisymmetric circulation in the superfluid interior Peralta et al. (2006); van Eysden and Melatos (2008); Bennett et al. (2010); Melatos et al. (2015).

Accreting neutron stars in binary systems are important search targets, because mass transfer spins up the star to $\gtrsim 10^{2}$ Hz and may simultaneously drive several quadrupole-generating mechanisms Bildsten (1998); Andersson et al. (1999); Nayyar and Owen (2006); Melatos (2007); Vigelius and Melatos (2009). Moreover it is observed that the distribution of spin frequencies of low-mass X-ray binaries (LMXBs) cuts off near 620 Hz Chakrabarty et al. (2003), below the theoretical centrifugal break-up limit $\approx 1.4$ kHz Cook et al. (1994). This has been explained by hypothesizing that the gravitational radiation-reaction torque balances the accretion torque Papaloizou and Pringle (1978); Wagoner (1984); Bildsten (1998), implying a relation between the X-ray flux and gravitational wave strain. Scorpius X-1 (Sco X-1), the most X-ray-luminous LMXB, is therefore a promising target for gravitational wave searches.

Initial LIGO achieved its design sensitivity over a wide band during LIGO Science Run 5 (S5) Abbott et al. (2009) and exceeded it during Science Run 6 (S6) The LIGO Scientific Collaboration and The Virgo Collaboration (2012). The strain sensitivity of the next-generation Advanced LIGO interferometer is expected to improve ten-fold relative to Initial LIGO after several stages of upgrade Harry and the LIGO Scientific Collaboration (2010). In the first observation run (O1), from September 2015 to January 2016, the strain noise is three to four times lower than in S6 across the most sensitive band, between 100 Hz and 300 Hz, and $\sim 30$ times lower around 50 Hz Abbott et al. (2016).

Four types of searches have been conducted for Sco X-1 using data collected by Initial LIGO and Advanced LIGO (O1). None of these searches reported a detection. First, a coherent search, using a maximum likelihood detection statistic called the $\mathcal{F}$-statistic Jaranowski et al. (1998), analysed the most sensitive six-hour data segment from Science Run 2 (S2) and placed a 95% confidence strain upper limit at $h_{0}^{95\%}\approx 2\times 10^{-22}$ for two bands, 464–484 Hz and 604–626 Hz Abbott et al. (2007a). Second, a directed, semi-coherent analysis based on the sideband algorithm was conducted on a 10-day stretch of LIGO S5 data in the band 50–550 Hz and reported median strain upper limits of $1.3\times 10^{-24}$ and $8\times 10^{-25}$ at 150 Hz for arbitrary and electromagnetically restricted source orientations, respectively Aasi et al. (2015b). The sideband method sums incoherently the coherent $\mathcal{F}$-statistic power at frequency-modulated orbital sidebands and generates a new detection statistic called the $\mathcal{C}$-statistic Messenger and Woan (2007); Sammut et al. (2014). Third, a directed version of the all-sky TwoSpect search Goetz and Riles (2011) was applied to S6 data and the second and third Virgo science runs (VSR2 and VSR3, respectively), yielding low-frequency upper limits of $h_{0}^{95\%}\approx 2\times 10^{-23}$ in the band from 20 Hz to 57.25 Hz Aasi et al. (2014). Another search of S6 data was carried out using the subsequently improved TwoSpect method Meadors et al. (2016), spanning frequencies from 40 Hz to 2040 Hz and projected semi-major axis from 0.90 s to 1.98 s. It achieved a 95% confidence level random-polarization upper limit of $h_{0}^{95\%}=1.8\times 10^{-24}$ at 165 Hz Meadors et al. (2017). Fourth, a directed version of the all-sky, radiometer search Ballmer (2005) was conducted on all 20 days of Science Run 4 (S4) data Abbott et al. (2007b), and was later applied to two years of S5 data, yielding a 90% confidence root-mean-square strain upper limit of $7\times 10^{-25}$ at 150 Hz Abadie et al. (2011), which converts to $h_{0}^{90\%}=2\times 10^{-24}$ Messenger (2011). The same method was applied to O1 data, yielding a median frequency-dependent limit of $h_{0}^{90\%}=6.7\times 10^{-25}$ at the most sensitive detector frequencies between 130–175 Hz Abbott et al. (2017).

It is probable that the spin frequency of Sco X-1 wanders stochastically under the fluctuating action of the hydromagnetic torque exerted by the accretion flow de Kool and Anzer (1993); Baykal and Oegelman (1993); Bildsten et al. (1997). Search methods that scan templates without guidance from a measured ephemeris are compromised because of spin wandering; for example the sideband search is restricted to a 10-day stretch of data in Ref. Aasi et al. (2015b), so that the signal power does not leak into adjacent frequency bins. Hidden Markov model (HMM) tracking offers a powerful strategy for detecting a spin-wandering signal Suvorova et al. (2016). An HMM relates a sequence of observations to the most probable Markov sequence of allowed transitions between the states of an underlying, hidden state variable (here the gravitational wave signal frequency $f_{0}$) Quinn and Hannan (2001). It can track $f_{0}$ over the total observation time $T_{\rm obs}$ by incoherently combining segments with duration $T_{\rm drift}=10$ d of the output from a maximum-likelihood, coherent matched filter , improving the sensitivity by a factor $\approx(T_{\rm obs}/T_{\rm drift})^{1/4}$ relative to a single segment.

In this paper, we combine the sideband algorithm with an HMM and apply it to Advanced LIGO O1 data. Specifically we carry out a directed search for Sco X-1 in the band 60–650 Hz. No evidence of a gravitational-wave signal is found. Frequentist 95% confidence strain upper limits of $h_{0}^{95\%}=4.0\times 10^{-25}$, $8.3\times 10^{-25}$, and $3.0\times 10^{-25}$ are derived at 106 Hz, for electromagnetically restricted source orientation, unknown polarization, and circular polarization, respectively. The paper is organized as follows. In Section II, we briefly review the search algorithm. In Section III, we discuss the astrophysical parameters of the source, search procedure, detection threshold and estimated sensitivity. Results of the search, including veto output, candidate follow-up, and gravitational wave strain upper limits are presented in Section IV. We discuss the torque-balance upper limit in Section V and conclude with a summary in Section VI.

II Method

In this section we briefly introduce the HMM formulation of frequency tracking and the Viterbi algorithm for solving the HMM in Section II.1 and Appendices A and B, respectively. A matched filter appropriate for a continuous-wave source in a binary is reviewed in Section II.2. A full description of the method can be found in Ref. Suvorova et al. (2016).

II.1 HMM tracking

A HMM is a finite state automaton, in which a hidden (unobservable) state variable $q(t)$ transitions between values from the set $\{q_{1},\cdots,q_{N_{Q}}\}$ at discrete times $\{t_{0},\cdots,t_{N_{T}}\}$, while an observable state variable $o(t)$ transitions between values from the set $\{o_{1},\cdots,o_{N_{O}}\}$. The probability that $q(t)$ jumps from state $q_{i}$ to state $q_{j}$ is given by the transition matrix $A_{q_{i}q_{j}}$. The likelihood that the hidden state $q_{i}$ gives rise to the observation $o_{j}$ is given by the emission probability $L_{o_{j}q_{i}}$. In this application, we map the discrete hidden states one-to-one to the frequency bins in the output of a frequency-domain estimator $G(f)$ (see Section II.2) computed over an interval of length $T_{\rm drift}$, with bin size $\Delta f_{\rm drift}=1/(2T_{\rm drift})$. The procedure for choosing $T_{\rm drift}$ is described in Appendix A.

For a Markov process, the probability that the hidden path $Q=\{q(t_{0}),\cdots,q(t_{N_{T}})\}$ gives rise to the observed sequence $O=\{o(t_{0}),\cdots,o(t_{N_{T}})\}$ is given by

[TABLE]

where $\Pi_{q_{i}}$ denotes the prior (see Appendix A). The classic Viterbi algorithm Viterbi (1967) provides a recursive, computationally efficient route to computing $Q^{*}(O)$, the path that maximizes $P(Q|O)$. The steps in the algorithm are specified in Appendix B; the number of operations is of order $(N_{T}+1)N_{Q}\ln N_{Q}$ Quinn and Hannan (2001). In this paper, we define a detection score $S$, such that the log likelihood of the optimal Viterbi path equals the mean log likelihood of all paths plus $S$ standard deviations, viz.

[TABLE]

with

[TABLE]

and

[TABLE]

where $\delta_{q_{i}}(t_{N_{T}})$ denotes the maximum probability of the path ending in state $q_{i}$ ($1\leq i\leq N_{Q}$) at step $N_{T}$ (see Appendix B), and $\delta_{q^{*}}{(t_{N_{T}})}$ is the likelihood of the optimal Viterbi path, i.e. $P[Q^{*}(O)|O]$.

II.2 Matched filter: Bessel-weighted $\mathcal{F}$-statistic

The emission probability $L_{o(t)q_{i}}$ is computed from a frequency-domain estimator $G(f)$ as described in Appendix A. In the context of continuous-wave searches, $G(f)$ is a matched filter. The optimal matched filter for a biaxial rotor with no orbital motion is the maximum-likelihood $\mathcal{F}$-statistic Jaranowski et al. (1998), which accounts for the rotation of the Earth and its orbit around the Solar System barycenter (SSB). When the source orbits a binary companion, the gravitational-wave signal frequency is modulated due to the orbital Doppler effect Ransom et al. (2003); Messenger and Woan (2007); Sammut et al. (2014). The $\mathcal{F}$-statistic power is distributed into approximately $M=2m+1$ orbital sidebands with $m=\text{ceil}(2\pi f_{0}a_{0})$, separated in frequency by $1/P$, where $f_{0}$ is the intrinsic gravitational wave frequency, $a_{0}$ is the light travel time across the projected semi-major axis of the orbit, $P$ is the orbital period, and $\text{ceil}(x)$ denotes the smallest integer greater than or equal to $x$. For a Keplerian orbit with zero eccentricity, the gravitational wave strain can be expanded in a Jacobi-Anger series as Abramowitz and Stegun (1964); Suvorova et al. (2016)

[TABLE]

where $J_{n}(z)$ is a Bessel function of order $n$ of the first kind. The mathematical form of (5) suggests a Bessel-weighted $\mathcal{F}$-statistic as the matched filter $G(f)$ for a biaxial rotor in a binary system, which can be expressed as the convolution Suvorova et al. (2016)

[TABLE]

where $B(f)$ is given by

[TABLE]

Compared to the $\mathcal{C}$-statistic, used in a previously published sideband search for Sco X-1 Sammut et al. (2014); Suvorova et al. (2016), where the factor $[J_{n}(2\pi fa_{0})]^{2}$ in (7) is replaced by unity, the Bessel-weighted matched filter recovers approximately $\sqrt{2}$ times more signal. It marshals more power into a single bin, producing a distinct spike with shoulders, instead of the relatively flat onion-dome peak produced by the $\mathcal{C}$-statistic. These characteristics facilitate Viterbi tracking (see Section IV A in Ref. Suvorova et al. (2016) for details). We leverage the existing, efficient, thoroughly tested $\mathcal{F}$-statistic software infrastructure in the LSC Algorithm Library Applications (LALApps)111https://www.lsc-group.phys.uwm.edu/daswg/projects/lalapps/ to compute $\mathcal{F}(f)$ in (6) Prix (2011).

III Implementation

In this section we introduce the electromagnetically measured source parameters of Sco X-1 (Section III.1), and describe the workflow of the pipeline (Section III.2), detection threshold (Section III.3), and search sensitivity (Section III.4).

III.1 Sco X-1 parameters

The sky position ($\alpha,\delta$), orbital elements ($a_{0},P$), and orientation angles ($\iota,\psi$) of Sco X-1 have been measured electromagnetically to various degrees of accuracy. The values and $1\sigma$ (68%) confidence level uncertainties are quoted in the top half of Table 1.

The published uncertainty in the orbital period, $\Delta P=0.0432$ s Galloway et al. (2014), restricts the coherent observation time to $T_{\rm drift}\leq 50$ d Sammut et al. (2014); Aasi et al. (2015b). Hence it is safe to take a single, fixed $P$ value when evaluating the $\mathcal{F}$-statistic, given that the coherent data stretches we analyse are limited to 10 d (20 d for follow up; see Section IV.1.4). The published uncertainty in the projected semi-major axis, inferred from the measured orbital velocity, is $\Delta a_{0}=0.18$ s Steeghs et al. (2002). In the previous S5 sideband search, it is demonstrated that taking a single, fixed $a_{0}$ value does not impact search sensitivity given this published uncertainty Sammut et al. (2014); Aasi et al. (2015b). However recent, unpublished research has revised the range of $a_{0}$ upwards to $0.36\,{\rm s}\leq a_{0}\leq 3.25$ s. This is because the orbital velocity is difficult to measure electromagnetically, and the previous measurement is based on searching for the optimal centre of symmetry in the accretion disk emission, yielding an estimated velocity of $40\pm 5$ km s*-1* Steeghs et al. (2002). The preliminary results from the more recent study, which uses Doppler tomography measurements and Markov Chain Monte-Carlo analysis for the velocity, show that the constraint on the orbital velocity is weaker, corresponding to a range from $10\,{\rm km}\,{\rm s}^{-1}$ to $90\,{\rm km}\,{\rm s}^{-1}$ Wang (2017); Wang et al. . It is shown in Section IV B of Ref. Suvorova et al. (2016) that, if the true value of $a_{0}$ differs from the estimated $a_{0}$ by 10%, it would produce an uncertainty in the estimated frequency of $\approx 0.001$ Hz. Moreover, the log likelihood of the optimal path decreases by $\sim 50\%$, if the true value of $a_{0}$ differs from the estimated $a_{0}$ by 25%. We search over the wider, unpublished range of $a_{0}$ with a resolution of 0.01805 s in order to preserve sensitivity. The orientation angles $\iota$ and $\psi$ are measured from the position angle of the Sco X-1 radio jets on the sky, assuming that the rotation axis of the neutron star is perpendicular to the accretion disk. In the previously published sideband search, two orientation priors are considered: (1) uniform distributions of $\cos\iota$ and $\psi$; and (2) distributions peaked around the observed values in the top half of Table 1.

The parameter space covered by the search is defined in the bottom half of Table 1. We assume uniform priors on both $f_{0}$ and $a_{0}$.

III.2 Workflow

The search is parallelized into 1-Hz sub-bands to assist with managing the relatively large volume of data involved. The sub-bands must be narrow enough, so that we can replace $f$ with the mean value $\bar{f}$ in each sub-band to a good approximation, in order to avoid recalculating $B(f)$ in every frequency bin. The sub-bands must also be wide enough to contain the width of the matched filter. Sub-bands of 1-Hz satisfy both of these requirements, and were also adopted in the S5 sideband search Aasi et al. (2015b).

The flow chart in Figure 1 summarises the procedural steps in the search pipeline. Firstly, the 30-min short Fourier transforms (SFTs) constituting the whole observation are divided into $N_{T}$ blocks, each of duration $T_{\rm drift}=10$ d. In each 1-Hz sub-band, the $\mathcal{F}$-statistic is computed for each block at the known sky location of the source. Next we compute the Bessel-weighted $\mathcal{F}$-statistic $G(f)$ from (6) and (7), taking $a_{0}$ and $P$ as inputs; that is, $G(f)$ is computed in $N_{f_{0}}$ frequency bins for each of the $N_{T}$ blocks. Theoretically the HMM hidden state variable is two-dimensional, because we search over $f_{0}$ and $a_{0}$. In practice $a_{0}$ varies imperceptibly on the time-scale $T_{\rm obs}$, so the algorithm is equivalent to multiple, independent, one-dimensional HMM searches over $f_{0}$ on a grid of $a_{0}$ values. The detection score and corresponding optimal Viterbi path are recorded in each 1-Hz sub-band. We evaluate the detection scores to identify candidates, judge whether or not they come from instrumental artifacts via a well-defined hierarchy of vetoes, and claim a detection or compute strain upper limits for sub-bands without candidates.

III.3 Threshold

We determine the Viterbi score threshold $S_{\rm th}$ for a given false alarm rate $\alpha_{\rm f}$ through Monte-Carlo simulations, such that searching data sets containing pure noise yields a fraction $\alpha_{\rm f}$ of positive detections with $S>S_{\rm th}$. SFTs containing pure Gaussian noise are generated for seven 1-Hz sub-bands, starting at 55 Hz, 155 Hz, 255 Hz, 355 Hz, 455 Hz, 555 Hz, and 650 Hz, with the same single-sided power spectral density (PSD) $S_{h}(f)$ as actual O1 data and with $T_{\rm obs}=130$ d. Searches are repeated for 100 noise realisations in each 1-Hz sub-band following the recipe in Fig. 1. We track 161 $a_{0}$ values from 0.361 s to 3.249 s, with resolution 0.01805 s, as for a real search. We find that the results depend weakly on the sub-bands: the mean $\langle S\rangle$ varies from 6.48 to 6.59, and the standard deviation $\sigma_{S}$ varies from 0.24 to 0.33. Combining the 700 realisations yields $S_{\text{th}}=7.34$ for $\alpha_{\rm f}=1\%$.

To check the influence of non-Gaussian noise on $S_{\rm th}$, we choose three 1-Hz sub-bands, starting at 157 Hz, 355 Hz, and 635 Hz, in O1 interferometer data and repeat the search for real noise. As we have no means of generating multiple, random, real-noise realisations from scratch, we take 100 different sky locations as background noise realisations. We find that $\langle S\rangle$ and $\sigma_{S}$ range from 6.36 to 6.38 and 0.27 to 0.34, respectively. These results match the output from Gaussian noise simulations to better than $\sim 3\%$, as does $S_{\rm th}$. Hence we set $S_{\text{th}}=7.34$ in the forthcoming analysis described in Section IV.

In the follow-up procedures in Section IV, we search a subset of the data either from a single interferometer with $T_{\rm obs}=130\,{\rm d}=13\,T_{\rm drift}$ or two interferometers with $T_{\rm obs}=60\,{\rm d}=6\,T_{\rm drift}$. To check the validity of $S_{\rm th}=7.34$ when searching a subset of the data, we run 400 trials of Gaussian noise simulations using data generated for a single interferometer with $T_{\rm obs}=130\,{\rm d}$ or two interferometers with $T_{\rm obs}=60\,{\rm d}$. The resulting $S_{\rm th}$ remains the same overall, and $\langle S\rangle$ and $\sigma_{S}$ range from 6.44 to 6.50 and 0.27 to 0.30, respectively, matching the output in the simulations with two interferometers and $T_{\rm obs}=130\,{\rm d}$ to better than $\sim 3\%$. Hence we keep $S_{\text{th}}=7.34$ fixed for the follow-up procedures in Section IV.

III.4 Sensitivity

Given the threshold $S_{\rm th}\,(\alpha_{\rm f}=1\%)=7.34$, we evaluate the characteristic wave strain yielding 95% detection efficiency (i.e. 5% false dismissal rate), denoted by $h_{0}^{95\%}$, through Monte-Carlo simulations with signals injected into Gaussian noise. The simulations are performed between 155–156 Hz, where the detectors are most sensitive, with $T_{\rm obs}=130$ d, $T_{\rm drift}=10$ d, $N_{T}=13$, $\sqrt{S_{h}}=1\times 10^{-23}\,{\rm Hz}^{-1/2}$, and source parameters copied from Table 1. We choose $T_{\rm obs}=130$ d to equal the duration of O1. The parameters ${f_{0}}_{\rm inj}$, ${a_{0}}_{\rm inj}$, $\cos\iota_{\rm inj}$, and $\psi_{\rm inj}$ are randomly chosen with a uniform distribution within the ranges 155.34565530–155.3456847 Hz, 0.36–3.25 s, 0.712107–0.726493, and 0–$2\pi$ rad, respectively. We obtain $h_{0}^{95\%}=3\times 10^{-25}$ for electromagnetically restricted orientation by assuming $\iota\approx 44^{\circ}$ Fomalont et al. (2001). In reality, the signal-to-noise ratio scales in proportion to $h_{0}^{\rm eff}$, given by

[TABLE]

rather than $h_{0}$ Jaranowski et al. (1998); Messenger et al. (2015). Hence we can convert the limiting wave strain to $h_{0}^{{\rm eff},95\%}\approx 0.74h_{0}^{95\%}$ using the value $\iota=44^{\circ}$. For $T_{\rm obs}$ fixed, we expect

[TABLE]

The latter scaling is verified by a group of injections in three other frequency bands (55–56 Hz, 355–356 Hz, and 649–650 Hz). Evaluating $S_{h}(f)$ from the O1 PSD, we plot $h_{0}^{95\%}$ versus $f_{0}$ as the blue dashed curve in Figure 2, which represents the 95% detection efficiency curve in Gaussian noise simulations.

In practice interferometer noise is non-Gaussian, and $T_{\rm obs}$ is less than 130 d (duty cycle $\approx 60\%$). To correct for this, we pick 53 1-Hz sub-bands, run 3000 injections in real O1 interferometer data, and compare the resulting $h_{0}^{95\%}$ to the blue dashed curve in Figure 2. The injected signal parameters are chosen in the same way as in the Gaussian noise simulation. In each sub-band tested, the resulting $h_{0}^{95\%}$ values from real O1 injections are plotted as gray stars in Figure 2. The correction factor $\kappa$ in each 1-Hz sub-band is defined as $h_{0}^{95\%}$, as marked by the gray star, divided by the value read off the blue dashed curve. The correction factors in 53 sub-bands fluctuate weakly, with mean $\langle\kappa\rangle=1.56$ and standard deviation $\sigma_{\kappa}=0.03$. We therefore apply the same $\kappa=1.56$ across the full search and adjust the blue dashed curve to give the red solid curve in Figure 2. The latter represents the characteristic wave strain for 95% detection efficiency as a function of frequency in real O1 data. We find that 2846 out of the 3000 O1 injections are detected with $S>S_{\rm th}$, yielding a detection rate of 94.87%, consistent with the targeted detection efficiency.

IV O1 analysis

In this section, we analyse data from the O1 observing run extending from 12 September 2015 to 19 January 2016 UTC (GPS time 1126051217 to 1137254417). The data are divided into 13 blocks, with $T_{\rm drift}=10$ d, and fed into the HMM tracker described in Section II and III.

Narrowband, instrumental noise lines (e.g. power line at 60 Hz, beam splitter violin mode, electronics, mirror suspension, calibration) and their harmonics can obscure astrophysical continuous-wave signals. At low frequencies between 25 Hz and 60 Hz, there are at least six known lines in each 1-Hz sub-band, and $\approx 2/3$ of the sub-bands contain more than 15 lines. Hence we do not search below 60 Hz, because the optimal paths returned by the HMM are dominated by difficult-to-model noise. The sensitivity of the method degrades, as the width $4\pi a_{0}f_{0}/P$ of the matched filter increases (see Section II.2). We terminate the search arbitrarily at $f_{0}=650$ Hz to keep $4\pi a_{0}f_{0}/P$ below $\approx 0.4$ Hz, which is almost half the width of a sub-band.

We record the first-pass candidates identified by the search in Figure 3. We then sift them through a systematic hierachy of vetoes as follows: (1) known instrumental line veto (Section IV.1.1), (2) single interferometer veto (Section IV.1.2), (3) $T_{\rm obs}/2$ veto (Section IV.1.3), and (4) $T_{\rm drift}$ veto (Section IV.1.4). The safety verification of the four-step veto procedure is described in Section IV.2. Table 2 lists the numbers of candidates surviving after each veto. No candidate survives all the vetoes and so we set upper limits on $h_{0}$. The strain upper limits are discussed in Section IV.3.

IV.1 Vetoes

IV.1.1 Known line veto

First-pass candidates with $S>S_{\rm th}=7.34$ (red dots) are plotted in Figure 3 as a function of $f_{0}$ and $a_{0}$ as estimated by the HMM. Each dot stands for a candidate in a 1-Hz sub-band. The colour of a dot indicates its associated $S$ value (higher $S$ in darker shade). The HMM returns an optimal path $f_{0}(t)$ whose wandering is too slight to be discerned visually in Figure 3. We take $f_{0}$ to equal the arithmetic mean of the min $f_{0}(t)$ and max $f_{0}(t)$ in the plot.

A candidate is vetoed, if $f_{0}(t)$ satisfies $|f_{0}(t)-f_{\rm line}|<4\pi a_{0}f_{0}/P$ anywhere on the path, where $f_{\rm line}$ is the frequency of a known instrumental noise line. We find that the line veto excludes 75% of the candidates. The 44 survivors are marked by green circles or blue squares in Figure 3. (The distinction between the green and blue symbols is discussed below.) One immediately notices that most of the red dots appear at $a_{0}\lesssim 0.5$ s for all $f_{0}$. This is because a narrower matched filter produces a higher score when it encounters a narrow noise line. A noise line that produces high $\mathcal{F}$-statistic values concentrated in a handful of frequency bins spreads out when convolved with the matched filter in (7) and contributes to every Bessel-weighted $\mathcal{F}$-statistic bin in the band $|f_{0}(t)-f_{\rm line}|<4\pi a_{0}f_{0}/P$. The Viterbi score computed from the log likelihood of the optimal path is normalized by the standard deviation of all the log likelihoods in a 1-Hz sub-band. It is higher if the $\mathcal{F}$-statistic output containing a noise line is convolved with a narrower matched filter (i.e. smaller $a_{0}$), because the $\mathcal{F}$-statistic-processed noise-line power is dispersed into fewer orbital sidebands. The plot confirms that most vetoed candidates have $a_{0}\lesssim 0.5$ s.

Instrumental lines are picked up readily by the HMM, rendering any astrophysical signal invisible in the relevant 1-Hz sub-band. One might seek to improve the search by notching out the instrumental lines first, before applying the HMM to the rest of the sub-band. However, O1 lines cluster closely below 90 Hz and near 300 Hz and 500 Hz, fragmenting the uncontaminated bands. It is onerous to circumvent the fragmentation, so we postpone this improvement to future searches, when better interferometer sensitivity will warrant the extra effort. In this search, we do not report results in a 1-Hz sub-band, if the optimal path intersects any instrumental line. In total 136 out of 591 1-Hz sub-bands are removed in this way.

IV.1.2 Single interferometer veto

We now examine the 44 candidates surviving the known line veto by searching data from H1 and L1 separately. The sensitivities of the two interferometers during O1 are comparable, implying either that an astrophysical signal should appear in both detectors if it is strong enough or that it cannot be detected in either detector but can be seen after combining data from both. In contrast, a candidate is more likely a noise artifact originating in a single detector if it is detected in one detector with higher $S$ than the original combined score $S_{\cup}$, while the other detector yields $S<S_{\rm th}$.

We can categorize survivors of the known line veto in Section IV.1.1 into four classes presented in Table 3.

Category A: Only one detector yields $S>S_{\rm th}$, equal to or higher than $S_{\cup}$; and the frequency estimated from the detector with $S\geq S_{\cup}$ is approximately equal to that obtained by combining both, with an absolute discrepancy less than $2\pi{a_{0}}_{\cup}{f_{0}}_{\cup}/P$, where ${a_{0}}_{\cup}$ and ${f_{0}}_{\cup}$ are the $a_{0}$ and $f_{0}$ estimated using both detectors. Typically we find that the absolute discrepancy is less than 0.01 Hz, even smaller than $2\pi{a_{0}}_{\cup}{f_{0}}_{\cup}/P$. Any astrophysical signal that is too weak to yield $S>S_{\rm th}$ in one detector is unavoidably obscured by the undocumented noise artifact in the other detector. Hence we veto candidates in category A.

Category B: Only one detector yields $S>S_{\rm th}$, equal to or higher than $S_{\cup}$, but the optimal path from the detector with $S\geq S_{\cup}$ occurs at $f_{0}$ with $|f_{0}-{f_{0}}_{\cup}|\geq 2\pi{a_{0}}_{\cup}{f_{0}}_{\cup}/P$ (denoted by $f_{0}\neq{f_{0}}_{\cup}$ in Table 3). It is possible that a real signal only shows up at ${f_{0}}_{\cup}$ after combining data from two detectors. Hence we keep candidates in category B for follow-up.

Category C: Both detectors yields $S\geq S_{\rm th}$. The candidate may come either from noise or from a real signal registering strongly in both detectors. Hence we keep candidates in category C for follow-up.

Category D: Both detectors yield $S<S_{\rm th}$ even though we have $S_{\cup}\geq S_{\rm th}$. A real signal may be too weak to register in either detector individually but rises above the noise, when the two detectors are combined. Hence we keep candidates in categories D for further examination.

Among the 44 candidates surviving the line veto, 38 in total are vetoed. They are marked by green circles in Figure 3. All of them only appear in H1. The remaining six candidates marked by blue squares need to be examined further manually. Four of them show higher scores in H1 and $S<S_{\rm th}$ in L1, but the estimated $f_{0}$ from H1 is different from that obtained by combining both detectors, falling into category B in Table 3. Two candidates, in the sub-bands 449–450 Hz and 593–594 Hz, fall into category D in Table 3, with $S<S_{\rm th}$ in both H1 and L1.

IV.1.3 $T_{\rm obs}/2$ veto

We now divide the observing run into two halves: 12 September 2015 to 20 November 2015 UTC (GPS time 1126051217 to 1132020365) and 20 November 2015 to 19 January 2016 UTC (GPS time 1132020366 to 1137254417). We search the halves separately in the six 1-Hz sub-bands containing veto survivors listed in Table 4, combining data from two interferometers. Similar to the criteria listed in Section IV.1.2, we veto a candidate, if it appears in one half, with $S\geq S_{\cup}$, but does not appear in the other half, and if the estimated $f_{0}$ value is approximately equal to the original value.

The three candidates near 459 Hz, 534 Hz, and 548 Hz appear in the first half with higher $S$ but not in the second half. The candidate near 395 Hz appears in the second half with higher $S$ but not in the first half. Each one of them is detected in the first or second half at a frequency approximately equal to the original estimated $f_{0}$ with absolute discrepancy less than 0.01 Hz.

In sub-bands 449 Hz and 593 Hz, neither of the two halves yields $S>S_{\rm th}$. These two candidates are marked by an asterisk in Table 4 and require further follow-up.

IV.1.4 $T_{\rm drift}$ veto

In general we can categorize any survivors of the $T_{\rm obs}/2$ veto into four groups with reference to the optimal paths detected in the original search. The groups are defined in Table 5. We expect $S$ to increase, as the block length $T_{\rm drift}$ increases, as long as $T_{\rm drift}$ remains shorter than the intrinsic spin-wandering time-scale. One could therefore imagine vetoing a candidate whose optimal Viterbi path does not wander significantly, if increasing $T_{\rm drift}$ up to the observed wandering time-scale does not increase $S$. However, based on our experience analysing injections (see Section IV.2), we adopt a more conservative approach to reduce the false dismissal rate from this veto step. Specifically, we veto a candidate whose optimal Viterbi path does not wander significantly, if increasing $T_{\rm drift}$ up to the observed wandering time-scale yields $S<S_{\rm th}$ (i.e. $S$ drops below threshold) and the optimal paths returned for the two $T_{\rm drift}$ values do not match. For a candidate whose optimal Viterbi path does wander significantly, we do not expect $S$ to increase with $T_{\rm drift}$, if the intrinsic spin-wandering time-scale is effectively shorter than $T_{\rm drift}$ already. Indeed, it is reasonable for a strongly wandering signal to disappear when tracked with longer $T_{\rm drift}$. On the rare occasion when this does happen, the candidate is likely to be a noise artifact. Candidates surviving the $T_{\rm drift}$ veto need to be followed up with more sensitive search pipelines (e.g. Cross-Correlation Whelan et al. (2015)).

In this search, the two survivors asterisked in Table 4 do not display strong spin wandering; they drift within three and one $f_{0}$ bins (see Figure 5 in Appendix C). Hence we expect $S$ to increase at approximately the same $f_{0}$, as $T_{\rm drift}$ increases all the way up to $T_{\rm obs}$. In fact we find that it suffices to consider $T_{\rm drift}=20$ d. The original and follow-up results are recorded in Table 6. For $T_{\rm drift}=20$ d, no path is detected with $S>S_{\rm th}$ at sub-bands 449 Hz and 593 Hz. The optimal Viterbi paths returned from $T_{\rm drift}=10$ d and 20 d are different in each of the two sub-bands, with an absolute discrepancy $\gtrsim 0.02$ Hz and $\gtrsim 1.03$ s for estimated $f_{0}$ and $a_{0}$, respectively. Normally the absolute uncertainties in the estimated values of $f_{0}$ and $a_{0}$ are less than 0.001 Hz and 0.02 s, respectively (see more details in Section III.1 and Section IV B of Ref. Suvorova et al. (2016)). Hence we do not see any evidence of a real astrophysical signal in these two outliers.

IV.2 Veto safety

The four-step veto procedure is verified with four synthetic signals injected into 120 d of Initial LIGO S5 data recolored to Advanced LIGO O1 noise and 200 signals injected into 130 d of O1 data. The signals feature low spin wandering, drifting within one to four $f_{0}$ bins during the full observation. We do not inject signals into the sub-bands contaminated by known noise lines, so these 204 signals survive the first veto step in Section IV.1.1 automatically. Only two out of the 204 injections are vetoed after the four steps described in Section IV.1.1–IV.1.4, yielding a false dismissal rate $<1\%$ and demonstrating that detectable spin-wandering signals are not commonly rejected. The two vetoed injections are rejected by the $T_{\rm obs}/2$ veto. They return a slightly higher $S$ value than $S_{\cup}$ (one in the first half, the other in the second), with $(S-S_{\cup})/S_{\cup}\leq 3\%$ and $S_{\cup}\lesssim 10$ (i.e. $<50\%$ higher than $S_{\rm th}$). In other words, the two false dismissals happen when both $(S-S_{\cup})/S_{\cup}$ and $S_{\cup}$ are small. By contrast, three out of the four candidates vetoed in Table 4 (Section IV.1.3) return $(S-S_{\cup})/S_{\cup}>10\%$ (with $8<S_{\cup}<16$), and the other returns $S-S_{\cup}=0.35$ with $S_{\cup}>20$ (i.e. $175\%$ higher than $S_{\rm th}$). Hence the four vetoed candidates in Table 4 fail the $T_{\rm obs}/2$ veto more strongly and are unlikely to be false dismissals.

Twelve examples of the synthetic signals surviving the vetoes described in Section IV.1.1–IV.1.4 are listed in Table 7.

IV.3 Strain upper limits

In the absence of a detection, we can place an upper limit on $h_{0}$ at a desired level of confidence (usually 95%) as a function of $f_{0}$.

A Bayesian analytic approach was adopted in the previous S5 sideband search for computing the strain upper limits Aasi et al. (2015b). However, the distribution of Viterbi path probabilities is hard to calculate analytically; Viterbi paths are correlated, and the nonlinear maximization step in the algorithm is hard to handle even within the context of extreme value theory (see Section III C in Ref. Suvorova et al. (2016)). Hence the Bayesian approach is hard to extend to the HMM sideband search. Instead, we adopt an empirical approach to set a frequentist upper limit as follows. We define $h_{0}^{u}$ such that the probability to detect a signal with $h_{0}\geq h_{0}^{u}$ is greater than or equal to $u$, i.e. $\Pr(S\geq S_{\rm th}|h_{0}\geq h_{0}^{u})\geq u$. Hence with no detection we take the $h_{0}^{95\%}$ value plotted in Figure 2 (see Section III.4) as the frequentist 95% confidence upper limit for electromagnetically restricted $\cos\iota$. It can be analytically converted to upper limits for unknown and circular polarizations using the scaling given by Equation (8).

Figure 4 displays the upper limit derived from the O1 search combining data from H1 and L1 as a function of $f_{0}$. Each marker indicates $h_{0}^{95\%}$ in the corresponding 1-Hz sub-band. Bands that do not contain a marker are those containing a candidate vetoed in any of the four veto stages described in Section IV.1.1–IV.1.4. In total 180 out of 591 1-Hz sub-bands contain vetoed candidates (see Table 2). The red dots correspond to assuming $\iota=44^{\circ}$, as inferred from radio observations Fomalont et al. (2001). The blue crosses correspond to assuming unknown polarization and a flat prior on $\cos\iota$. The cyan triangles correspond to assuming circularly polarized signals (i.e. $\cos\iota=\pm 1$). At 106 Hz, the lowest 95% confidence upper limits are $h_{0}^{95\%}=4.0\times 10^{-25}$, $8.3\times 10^{-25}$ and $3.0\times 10^{-25}$ for electromagnetically restricted $\cos\iota$, unknown polarization, and circular polarization, respectively. Hence the electromagnetically restricted prior and circular polarization assumptions improve upon the upper limits for unknown polarization by factors of 2.08 and 2.77, respectively.

As a further check, we compare the frequentist Viterbi upper limit to the frequentist $\mathcal{C}$-statistic upper limit. We run injections in six 1-Hz sub-bands in the best 10-day stretch of the real O1 interferometer data, starting from 110 Hz, 257 Hz, 355 Hz, 454 Hz, 550 Hz and 649 Hz, and search for them with the $\mathcal{C}$-statistic sideband pipeline Sammut et al. (2014); Aasi et al. (2015b). The best 10-day data stretch is selected from O1 as follows Wette (2009); Abadie et al. (2010). A figure of merit, proportional to the signal-to-noise ratio (SNR), is defined by $\sum_{I,J}[S_{h}(f_{I})]_{J}^{-1}$, where $[S_{h}(f_{I})]_{J}$ is the strain noise power spectral density at discrete frequency bin $f_{I}$ in the $J^{\text{th}}$ SFT, and the summation is over all SFTs in each rolling 10-day stretch in O1. The 10-day data stretch with the highest value of this figure over the 60–650 Hz band is selected. We compare the values of $h_{0}^{95\%}$ from the $\mathcal{C}$-statistic to the values plotted in Figure 4. The results show that the frequentist 95% confidence upper limits from the $\mathcal{C}$-statistic are 1.46–1.74 times larger than those achieved from the search described in this paper.

V Torque-balance upper limit

In LMXBs the gravitational wave strain inferred from the torque-balance scenario can be expressed as a function of the spin frequency of the neutron star $f_{\star}$ and the the X-ray flux $F_{X}$ according to Bildsten (1998); Wagoner (1984); Sammut et al. (2014)

[TABLE]

where $R_{\star}$ is the stellar radius and $M_{\star}$ is the stellar mass.222We assume that the system emits gravitational radiation via the mass quadrupole channel. The analogous equation for current quadrupole radiation is given in Ref. Owen (2010). We now ask how $h_{0}^{\rm eq}$ compares to the results of the analysis in Section IV.

Let us take the electromagnetically measured $F_{X}=4\times 10^{-7}$ erg cm*-2* s*-1* Watts et al. (2008) for Sco X-1 and the fiducial values $R_{\star}=10$ km and $M_{\star}=1.4M_{\odot}$. We plot $h_{0}^{\rm eq}$ as a function of $f_{0}=2f_{\star}$ in Figure 4 (green solid curve). Near 106 Hz, where the best $h_{0}^{95\%}$ is reported, we obtain $h_{0}^{\rm eq}\approx 8.3\times 10^{-26}$, which is $4.8$, $10.0$, and $3.6$ times lower than $h_{0}^{95\%}$ for electromagnetically restricted $\cos\iota$, unknown polarization, and circular polarization, respectively. The design sensitivity of Advanced LIGO is expected to improve further about two-fold relative to O1 Abbott et al. (2016). The anticipated $h_{0}^{95\%}$ at the design sensitivity of Advanced LIGO is plotted as a function of $f_{0}$ in Figure 4 as the red curve, assuming an electromagnetically restricted orientation ($\iota=44^{\circ}$) and $T_{\rm obs}=2$ yr. Near 50 Hz, $h_{0}^{95\%}$ reaches $h_{0}^{\rm eq}$.

The green solid curve in Figure 4 is somewhat conservative Aasi et al. (2015b). If we consider the Alfvén radius to be the accretion-torque lever arm, instead of $R_{\star}$ as assumed in (10), then $h_{0}^{\rm eq}$ increases by a factor of a few. The Alfvén radius is given by Bildsten et al. (1997)

[TABLE]

where $B_{\star}$ is the magnetic field of the star, $G$ is Newton’s gravitational constant, and $\dot{M}$ is the accretion rate. The neutron stars in LMXBs have $\dot{M}$ ranging from $\sim 10^{-11}M_{\odot}\,{\rm yr}^{-1}$ to the Eddington limit $2\times 10^{-8}M_{\odot}\,{\rm yr}^{-1}$ Ritter and Kolb (2003); Sammut (2015), and weak magnetic fields in the range $10^{8}\,{\rm G}\lesssim B_{\star}\lesssim 10^{9}$ G Bildsten (1998); Patruno and Watts (2012); Sammut (2015). To estimate the maximum magnitude of the effect, we substitute $\dot{M}=10^{-8}M_{\odot}\,{\rm yr}^{-1}$ and $B_{\star}=10^{9}$ G in Equation (12). The resulting $h_{0}^{\rm eq}$ is shown as the green dashed curve in Figure 4, giving $h_{0}^{95\%}\approx 2h_{0}^{\rm eq}$ for electromagnetically restricted $\cos\iota$. At the design sensitivity of Advanced LIGO, we expect $h_{0}^{95\%}<h_{0}^{\rm eq}$ in the band $30\,{\rm Hz}\lesssim f_{0}\lesssim 250\,{\rm Hz}$.

VI Conclusion

We perform an HMM sideband search for continuous gravitational waves from Sco X-1 in Advanced LIGO O1 data from 60 Hz to 650 Hz. The analysis is computationally efficient, requiring $\lesssim 3\times 10^{3}$ CPU-hr. We see no evidence of gravitational waves. Frequentist 95% confidence upper limits of $h_{0}^{95\%}=4.0\times 10^{-25}$, $8.3\times 10^{-25}$, and $3.0\times 10^{-25}$ are derived at 106 Hz for electromagnetically restricted $\cos\iota$, unknown polarization, and circular polarization, respectively. The upper limits are derived from Monte-Carlo simulations of spin-wandering signals. They are $4.8$, $10.0$, and $3.6$ times larger than the stellar radius torque-balance limit $h_{0}^{\rm eq}$, and approach $h_{0}^{\rm eq}$ more closely, if we treat the Alfvén radius as the accretion-torque lever arm. An analysis of two years of Advanced LIGO data at design sensitivity with this search will be able to constrain the Alfvén radius lever-arm scenario at frequencies below 300 Hz. The best existing Bayesian 90% confidence median strain upper limit from the radiometer O1 search is $h_{0}^{90\%}=6.7\times 10^{-25}$ at 135 Hz Abbott et al. (2017). It converts to 95% confidence median and maximum upper limits $h_{0}^{95\%}=7.8\times 10^{-25}$ and $h_{0}^{95\%}=1.0\times 10^{-24}$, respectively in the sub-band 134–135 Hz Messenger (2011), which are comparable to the results for unknown polarization presented here.333The value of $h_{0}^{95\%}$ from the present search for unknown polarization is 6% higher and 17% lower than the median and maximum $h_{0}^{95\%}$ values from the radiometer search, respectively Abbott et al. (2017). A direct comparison of the best quoted limits from the present search and the radiometer search is complicated by the different approaches of reporting upper limits. The present search returns the optimal Viterbi path (i.e. one upper limit) in each 1-Hz sub-band, while the radiometer search reports a range of upper limits. Although these results are similar in sensitivity, this is the first analysis that searches over the projected semi-major axis of the binary orbit within the uncertainty of the electromagnetic measurement, while taking into account the effects of spin wandering over $T_{\rm obs}$. The spin frequency of Sco X-1 has not been determined conclusively, and could also lie below 60 Hz. In the future, it is hoped that the number of instrumental lines at low frequencies will be reduced, enabling analysis below 60 Hz, where $h_{0}^{\rm eq}$ is higher and hence easier to reach. At the design sensitivity of Advanced LIGO, it is anticipated that $h_{0}^{95\%}$ can be improved further by a factor of 2–3, reaching $h_{0}^{\rm eq}$ near 50 Hz. In addition to Sco X-1, the search can be applied to other X-ray binaries including Cygnus X-3, the next brightest X-ray source after Sco X-1, and sources like XTE J1751-305 and 4U 1636-536, which show periodicities in the X-ray light curves and may indicate r-mode oscillations The LIGO Scientific Collaboration and The Virgo Collaboration (2016); Mahmoodifar and Strohmayer (2013); Haskell (2015).

VII Acknowledgements

The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific and Industrial Research of India, Department of Science and Technology, India, Science & Engineering Research Board (SERB), India, Ministry of Human Resource Development, India, the Spanish Ministerio de Economía y Competitividad, the Vicepresidència i Conselleria d’Innovació, Recerca i Turisme and the Conselleria d’Educació i Universitat del Govern de les Illes Balears, the National Science Centre of Poland, the European Commission, the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the Lyon Institute of Origins (LIO), the National Research Foundation of Korea, Industry Canada and the Province of Ontario through the Ministry of Economic Development and Innovation, the Natural Science and Engineering Research Council Canada, Canadian Institute for Advanced Research, the Brazilian Ministry of Science, Technology, and Innovation, International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), Russian Foundation for Basic Research, the Leverhulme Trust, the Research Corporation, Ministry of Science and Technology (MOST), Taiwan and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, MPS, INFN, CNRS and the State of Niedersachsen/Germany for provision of computational resources. This is LIGO document LIGO-P1700019.

Appendix A Hidden Markov model

An HMM is a finite state automaton defined by a hidden (unobservable) state variable $q(t)$ transitioning between values from the set $\{q_{1},\cdots,q_{N_{Q}}\}$ and an observable state variable $o(t)$ taking values from the set $\{o_{1},\cdots,o_{N_{O}}\}$ at discrete times $\{t_{0},\cdots,t_{N_{T}}\}$. The automaton jumps between hidden states from $t_{n}$ to $t_{n+1}$ with probability

[TABLE]

and is observed in the state $o_{j}$ with emission probability

[TABLE]

For a Markov process, the probability that the hidden path $Q=\{q(t_{0}),\cdots,q(t_{N_{T}})\}$ gives rise to the observed sequence $O=\{o(t_{0}),\cdots,o(t_{N_{T}})\}$ is given by

[TABLE]

where

[TABLE]

is the prior. The most probable path $Q^{*}(O)=\arg\max P(Q|O)$ maximizes $P(Q|O)$ and gives the best estimate of $q(t)$ over the total observation.

In this application, we map the discrete hidden states one-to-one to the frequency bins in the output of a frequency-domain estimator $G(f)$ (see Section II.2) computed over an interval of length $T_{\rm drift}$, with bin size $\Delta f_{\rm drift}=1/(2T_{\rm drift})$. We can always choose an intermediate time-scale $T_{\rm drift}$ in between the duration of one SFT, $T_{\rm SFT}=30$ min, and the total observation time $T_{\rm obs}$ in order to satisfy

[TABLE]

for all $t$.444Frequency-domain, continuous-wave LIGO searches operate on SFTs rather than the time series of the detector output Riles (2013). We assume that the spin wandering caused by accretion noise in Sco X-1 follows an unbiased Wiener process, in which $f_{0}(t)$ experiences a random walk and stays within $\Delta f_{\rm drift}$ for a duration less than a conservatively chosen $T_{\rm drift}=10$ d, based on the assumption that the deviation of the accretion torque from its average value flips sign on the time-scale of observed fluctuations in the X-ray flux Ushomirsky et al. (2000); Aasi et al. (2015b).555For constant spin up or spin down, we are able to track a maximum rate $|\dot{f_{0}}|=\Delta f_{\rm drift}T_{\rm drift}^{-1}=7\times 10^{-13}$ Hz s*-1*. By way of comparison, without considering accretion noise, the secular spin-down (or spin-up) rate of LMXBs satisfies $|\dot{f_{0}}|\lesssim 10^{-14}$ Hz s*-1* Patruno and Watts (2012). Assuming continuous frequency wandering (i.e. no neutron star rotational glitches), equation (13) simplifies to the tridiagonal form

[TABLE]

with all other entries vanishing. The emission probability can be expressed in terms of $G(f)$ as

[TABLE]

where $G({f_{0}}_{i})$ is the log likelihood that the gravitational-wave signal frequency $f_{0}$ (e.g. twice the spin frequency of the star) lies in the frequency bin $[{f_{0}}_{i},{f_{0}}_{i}+\Delta f_{\rm drift}]$ during the interval $[t,t+T_{\rm drift}]$. As we have no advance knowledge of $f_{0}$, we choose a uniform prior, viz.

[TABLE]

Appendix B Viterbi algorithm

The classic Viterbi algorithm Viterbi (1967) provides a recursive, computationally efficient route to computing $Q^{*}(O)$, reducing the number of operations to $(N_{T}+1)N_{Q}\ln N_{Q}$ by binary maximization Quinn and Hannan (2001). At every forward step $k$ ($1\leq k\leq N_{T}$) in the recursion, the algorithm eliminates all but $N_{Q}$ possible state sequences, and stores the $N_{Q}$ maximum probabilities ($1\leq i\leq N_{Q}$)

[TABLE]

It also stores the previous-step states of origin,

[TABLE]

that maximize the probability at that step. The optimal Viterbi path is then reconstructed by backtracking according to

[TABLE]

for $0\leq k\leq N_{T}-1$. A detailed description of the algorithm can be found in Section II D of Ref. Suvorova et al. (2016).

Appendix C $T_{\rm obs}/2$ veto survivors: optimal Viterbi paths

In the $T_{\rm drift}$ veto described in Section IV.1.4, we categorize the two survivors according to their optimal paths detected in the original search. The optimal paths of the two survivors are plotted in Figure 5, showing the estimated frequency $f_{0}$ as a function of time evaluated at the endpoint of each Viterbi step. The paths near 449 Hz and 593 Hz drift within three and one $f_{0}$ bins, respectively over $T_{\rm obs}$. They display low spin wandering.