A directed search for gravitational waves from Scorpius X-1 with initial LIGO
We present results of a search for continuously-emitted gravitational radiation, directed at the brightest low-mass X-ray binary, Scorpius X-1. Our semi-coherent analysis covers 10 days of LIGO S5 data ranging from 50–550 Hz, and performs an incoherent sum of coherent -statistic power distributed amongst frequency-modulated orbital sidebands. All candidates not removed at the veto stage were found to be consistent with noise at a false alarm rate. We present Bayesian 95% confidence upper limits on gravitational-wave strain amplitude using two different prior distributions: a standard one, with no a priori assumptions about the orientation of Scorpius X-1; and an angle-restricted one, using a prior derived from electromagnetic observations. Median strain upper limits of and are reported at 150 Hz for the standard and angle-restricted searches respectively. This proof of principle analysis was limited to a short observation time by unknown effects of accretion on the intrinsic spin frequency of the neutron star, but improves upon previous upper limits by factors of for the standard, and 2.3 for the angle-restricted search at the sensitive region of the detector.
Recycled neutron stars are a likely source of persistent, quasi-monochromatic gravitational waves detectable by ground-based interferometric detectors. Emission mechanisms include thermocompositional and magnetic mountains Bonazzola and Gourgoulhon (1996); Ushomirsky et al. (2000); Melatos and Payne (2005); Vigelius and Melatos (2009); Priymak et al. (2011); Haskell et al. (2006), unstable oscillation modes Owen et al. (1998) and free precession Jones (2002). If the angular momentum lost to gravitational radiation is balanced by the spin-up torque from accretion, the gravitational wave strain can be estimated independently of the microphysical origin of the quadrupole and is proportional to the observable X-ray flux and spin frequency Wagoner (1984); Bildsten (1998) via . Given the assumption of torque balance, the strongest gravitational wave sources are those that are most proximate with the highest accretion rate and hence X-ray flux, such as low-mass X-ray binary (LMXB) systems. In this sense the most luminous gravitational wave LMXB source is Scorpius X-1 (Sco X-1).
The plausibility of the torque-balance scenario is strengthened by observations of the spin frequencies of pulsating or bursting LMXBs, which show them clustered in a relatively narrow band from Hz Chakrabarty et al. (2003), even though their ages and accretions rates imply that they should have accreted enough matter to reach the centrifugal break-up limit Hz Cook et al. (1994) of the neutron star. The gravitational wave spin-down torque scales as , mapping a wide range of accretion rates into a narrow range of equilibrium spins, so far conforming with observations. Alternative explanations for the clustering of LMXB spin periods involving disc accretion physics have been proposed Haskell et al. (2014). Although this explanation suggests that gravitational radiation is not required to brake the spin-up of the neutron star, it does not rule out gravitational emission from these systems. The gravitational-wave torque-balance argument is used here as an approximate bound.
The initial instruments installed in the Laser Interferometer Gravitational Wave Observatory (LIGO) consisted of three Michelson interferometers, one with 4-km orthogonal arms at Livingston, LA, and two collocated at Hanford, CA, with 4 km and 2 km arms. Initial LIGO achieved its design sensitivity during its fifth science run (between November 2005 and October 2007 ) Abbott et al. (2009); Abadie et al. (2010) and is currently being upgraded to the next-generation Advanced LIGO configuration, which is expected to improve its sensitivity ten-fold in strain Harry and the LIGO Scientific Collaboration (2010).
Three types of searches have previously been conducted with LIGO data for Sco X-1. The first, a coherent analysis using data from LIGO’s second science run (S2), was computationally limited to six-hour data segments. It placed a wave-strain upper limit at confidence of for two 20 Hz bands between 464 – 484 Hz and 604 – 626 Hz Abbott et al. (2007). The second, employing a radiometer technique Ballmer (2006), was conducted using all 20 days of LIGO S4 data Abbott et al. (2007). It improved the upper limits on the previous (S2) search by an order of magnitude in the relevant frequency bands but did not yield a detection. The same method was later applied to S5 data and reported roughly a five-fold sensitivity improvement over the S4 results Abadie et al. (2011). The S5 analysis returned a median confidence root-mean-square strain upper limit of at 150 Hz, the most sensitive detector frequency (this converts to str (); rad (); Messenger (2010)). Thirdly, an all-sky search for continuous gravitational waves from sources in binary systems, which looks for patterns caused by binary orbital motion doubly Fourier-transformed data (TwoSpect), was adapted to search the Sco X-1 sky position, and returned results in the low frequency band from Hz Aasi et al. (2014).
Here we implement a new search for gravitational waves from sources in known binary systems, with unknown spin frequency, initially directed at Sco X-1 on LIGO S5 data to demonstrate feasibility. Values of the coherent, matched-filtered -statistic Jaranowski et al. (1998) are incoherently summed at the locations of frequency-modulated sidebands. This multi-stage, semi-coherent, analysis yields a new detection statistic, denoted the -statistic Messenger and Woan (2007); Sammut et al. (2014). A similar technique was first employed in electromagnetic searches for radio pulsars Ransom et al. (2003). We utilise this technique to efficiently deal with the large parameter space introduced by the orbital motion of a source in a binary system.
A brief description of the search is given in Sec. II, while the astrophysical target source and its associated parameter space are discussed in Sec. III. Section IV outlines the search method, reviews the pipeline, discusses the selection and preprocessing of LIGO S5 data, and explains the post-processing procedure. Results of the search, including upper limits of gravitational wave strain, are presented and discussed in Sections V and VI, respectively and restated in Sec. VII.
Ii Search Method
For a gravitational wave source in a binary system, the frequency of the signal is Doppler modulated by the orbital motion of the source with respect to the Earth Ransom et al. (2003); Messenger and Woan (2007); Sammut et al. (2014). The semi-coherent sideband search method involves the incoherent summation of frequency modulated sidebands of the coherent -statistic Jaranowski et al. (1998); Messenger and Woan (2007); Sammut et al. (2014).
The first step in the sideband search is to calculate the coherent -statistic as a function of frequency, assuming only a fixed sky position. Knowing the sky position, one can account for the phase evolution due to the motion of the detector. For sources in binary systems, the orbital motion splits the signal contribution to the -statistic into approximately sidebands separated by in frequency, where 111The function rounds up to the nearest integer., is the intrinsic gravitational wave frequency, is the light travel time across the semi-major axis of the orbit, and is the orbital period. Knowledge of and , allows us to construct an -statistic sideband template.
The second stage of the sideband pipeline is the calculation of the -statistic, where we convolve the sideband template with the coherent -statistic. The result is an incoherent sum of the signal power at each of the potential sidebands as a function of intrinsic gravitational wave frequency. For our template we use a flat comb function with equal amplitude teeth (see Fig. 1 of Sammut et al. (2014)), and hence, for a discrete frequency bin , the template is given by
where depends on search frequency and the semi-major axis used to construct the template (see Sec. III.4) 222The frequency is the central frequency of the search sub-band. Although the pipeline is designed for a wide-band search, the total band should be broken up into sub-bands narrow enough that the number of sidebands in the template, , does not change significantly from the lower to the upper edge of the sub-band.. The index of the Kronecker delta-function is defined as
for a frequency bin width , where round() returns the closest integer, and denotes our best guess at the orbital period. The following convolution then yields the -statistic,
where the form of is described in Jaranowski et al. (1998); Prix (2007); Sammut et al. (2014). We therefore obtain this final statistic as a function of frequency evaluated at the same discrete frequency bins on which the input -statistic is computed.
Iii Parameter space
Sco X-1 is the brightest LMXB, and the first to be discovered in 1962 Giacconi et al. (1962), located 2.8 kpc away Bradshaw et al. (1999), in the constellation Scorpius. Source parameters inferred from a variety of electromagnetic measurements are displayed in Table 1. Assuming that the gravitational radiation and accretion torques balance, we obtain an indirect upper limit on the gravitational wave strain amplitude for Sco X-1 as a function of 333For a mass quadrupole , while for a current quadrupole .. Assuming fiducial values for the mass , radius km Steeghs and Casares (2002), and moment of inertia kg m Bradshaw et al. (1999) gives
Equation 5 assumes that all the angular momentum due to accretion is transferred to the star and converted into gravitational waves, providing an upper limit on the gravitational wave strain 444The torque-balance argument implies no angular momentum is lost from the system through the radio jets for example..
|Parameter (Name and Symbol)||Value [Reference]|
|X-ray flux||erg cm s||Watts et al. (2008)|
|Distance||kpc||Bradshaw et al. (1999)|
|Right ascension||16h 19m 55.0850s||Bradshaw et al. (1999)|
|Declination||38’ 24.9”||Bradshaw et al. (1999)|
|Sky Position angular resolution||0.3 mas||Bradshaw et al. (1999)|
|Proper motion||14.1 mas yr||Bradshaw et al. (1999)|
|Orbital period||s||Galloway et al. (2014)|
|Projected semi-major axis||1.44 s||Steeghs and Casares (2002)|
|Polarization angle||Fomalont et al. (2001)|
|Inclination angle||Fomalont et al. (2001)|
Optical observations of Sco X-1 have accurately determined its sky position and orbital period and, less accurately, the semi-major axis Gottlieb et al. (1975); Bradshaw et al. (1999); Steeghs and Casares (2002); Galloway et al. (2014). The rotation period remains unknown, since no X-ray pulsations or bursts have been detected. Although twin kHz quasi-periodic oscillations have been observed in the contiuous X-ray flux with separations in the range Hz, there is no consistent and validated method that supports a relationship between the QPO frequencys and the spin frequency of the neutron star (see Watts (2012) for a review). We therefore assume the spin period is unknown and search over a range of . We also assume a circular orbit, which is expected by the time mass transfer occurs in LMXB systems. In general, orbital eccentricity causes a redistribution of signal power amongst the existing circular orbit sidebands and will cause negligible leakage of signal power into additional sidebands at the boundaries of the sideband structure. Orbital eccentricity also has the effect of modifying the phase of each sideband. However, the standard sideband search is insensitive to the phase of individual sidebands.
This section defines the parameter space of the sideband search, quantifying the accuracy with which each parameter is and/or needs to be known. The parameters and their uncertainties are summarised in Table 2.
|Spin limited observation time||13 days|
|Period limited observation time 555at kHz||50 days|
|Maximum sky position error||300 mas|
|Maximum proper motion ††footnotemark: 666for days||3000 mas yr|
|Neutron star inclination 777from EM observations|
|- EM independent|
|Gravitational wave polaristion ††footnotemark:|
|- EM independent|
iii.1 Spin frequency
The (unknown) neutron star spin period is likely to fluctuate due to variations in the accretion rate . The coherent observation time span determines the size of the frequency bins in the calculation of the -statistic, along with an over-resolution factor defined such that a frequency bin is Hz wide. To avoid sensitivity loss due to the signal wandering outside an individual frequency bin, we restrict the coherent observation time to less than the spin limited observation time so that the signal is approximately monochromatic. Conservatively, assuming the deviation of the accretion torque from the mean flips sign randomly on the timescale days Ushomirsky et al. (2000), experiences a random walk which would stay within a Fourier frequency bin width for observation times less than given by
where is the moment of inertia of a neutron star with mass and radius , is the mean accretion torque and is the gravitational constant.
For Sco X-1, with fiducial values for , , and as described earlier, and assuming day (comparable to the timescale of fluctuations in X-ray flux Bildsten et al. (1997)), days. We choose observation time span days to fit safely within this restriction.
iii.2 Orbital period
The orbital period sets the frequency spacing of the sidebands. Uncertainties in this parameter will therefore translate to offsets in the spacing between the template and signal sidebands. The maximum coherent observation timespan allowed for use with a single template value of is determined by the uncertainty and can be expressed via
where a frequency bin in the -statistic has width as explained above, and and are the intrinsic gravitational wave frequency and light crossing time of the projected semi-major axis respectively. For a Sco X-1 search with at kHz, one finds days, longer than the maximum duration allowed by spin wandering (i.e. ). Choosing kHz gives a conservative limit for since lower frequencies will give higher values, and we only search up to 550 Hz. Thus we can safely assume the orbital period is known exactly for a search spanning 50 days.
iii.3 Sky position and proper motion
Knowledge of the source sky position is required to demodulate the effects of detector motion with respect to the barycentre of the source binary system (due to the Earth’s diurnal and orbital motion) when calculating the -statistic. We define an approximate worst-case error in sky position as that which would cause a maximum gravitational wave phase offset of 1 rad, giving us
where is the Earth-Sun distance (1 AU). Additionally, the proper motion of the source also needs to be taken into account. If the motion is large enough over the observation time it will contribute to the phase error in the same way as the sky position error. The worst case proper motion can therefore be determined similarly, viz.
For a 10-day observation at kHz, one finds mas and .
The sky position of Sco X-1 has been measured to within 0.3 mas, with a proper motion of 14.1 mas yr Bradshaw et al. (1999). These are well within the allowed constraints, validating the approximation that the sky position can be assumed known and fixed within our analysis.
iii.4 Semi-major axis
The semi-major axis determines the number of sidebands in the search template. Its uncertainty affects the sensitivity of the search independently of the observation time. To avoid the template width being underestimated, we construct a template using a semi-major axis given by the (best guess) observed value and its uncertainty such that
thus minimising signal losses. For a justification of this choice for , see Section IV. D. in Sammut et al. (2014).
iii.5 Inclination and polarisation angles
The inclination angle of the neutron star is the angle the spin axis makes with respect to the line of sight. Without any observational prior we would assume that the orientation of the spin axis is drawn from an isotropic distribution, and therefore comes from a uniform distribution within the range . The polarisation angle describes the orientation of the gravitational wave polarisation axis with respect to the equatorial coordinate system, and can be determined from the position angle of the spin axis, projected on the sky. Again, with no observational prior we assume that comes from a uniform distribution within the range .
The orientation angles and affect both the amplitude and phase of the incident gravitational wave. The phase contribution can be treated separately from the binary phase and the uncertainty in both and are analytically maximised within the construction of the -statistic. However, electromagnetic observations can be used to constrain the prior distributions on and . This information can be used to improve search sensitivity in post-processing when assessing the response of the pipeline to signals with parameters drawn from these prior distributions.
In this paper, we consider two scenarios for and : (i) uniform distributions within the previously defined ranges; and (ii) prior distributions based on values and uncertainties obtained from electromagnetic observations. From observations of the radio jets from Sco X-1 Fomalont et al. (2001) we can take , assuming the rotation axis of the neutron star is perpendicular to the accretion disk. The same observations yield a position angle of the radio jets of . Again, assuming alignment of the spin and disk normal, the position angle is directly related to the gravitational wave polaristion angle with a phase shift of , such that . For these observationally motivated priors we adopt Gaussian distributions, with mean and variance given by the observed values and their errors, respectively, as quoted above.
iv.1 Data selection
LIGO’s fifth science run (S5) took place between November 4, 2005 and October 1, 2007. During this period the three LIGO detectors (L1 in Livingston, LA; H1 and H2 collocated in Hanford, WA) achieved approximately one year of triple coincidence observation, operating near their design sensitivity Abadie et al. (2010). The amplitude spectral density of the strain noise of the two 4-km detectors (H1 and L1) was a minimum at 140 Hz and over the 100-300 Hz band.
Unkown effects of accretion on the rotation period of the neutron star (spin wandering) restricts the sideband analysis to a 10-day coherent observation time span (Section III.1). A 10-day data-stretch was selected from S5 as follows Abadie et al. (2010). A figure of merit, proportional to the signal-to-noise-ratio (SNR) and defined by , where is the strain noise power spectral density at frequency in the short Fourier transform (SFT), was assigned to each rolling 10-day stretch. The highest value of this quantity over the 100–300 Hz band (the region of greatest detector sensitivity) was achieved in the interval August 21–31, 2007 (GPS times 871760852–872626054) with duty factors of 91 % in H1 and 72% in L1. This data stretch was selected for the search. We search a 500 Hz band, ranging from 50-550 Hz, chosen to include the most sensitive region of the detector. The power spectral density for this stretch of data is shown in Fig. 1. The most prominent peaks in the noise spectrum are due to power line harmonics at 60 Hz and thermally excited violin modes from 330–350 Hz caused by the mirror suspension wires in the interferometer spe ().
Science data (data that excludes detector down-time and times flagged with poor data quality) are calibrated to produce a strain time series , which is then broken up into shorter segments of equal length. Some data are discarded, as not every continuous section of covers an integer multiple of segments. The segments are high-pass filtered above 40 Hz and Fourier transformed to form SFTs. For this search, 1800 sec SFTs are fed into the -statistic stage of the pipeline.
A flowchart of the multi-stage sideband pipeline is depicted in Fig. 2. After data selection, the first stage of the pipeline is the computation of the -statistic Prix (2011); lal () 888Note that there is no thresholding applied to the -statistic.. For the sideband search only the sky position is required at the -statistic stage, where the matched filter models an isolated source.
The outputs of the -statistic analysis are values of for each frequency bin from which the sideband algorithm then calculates the -statistic Sammut et al. (2014); lal (). The algorithm takes values of the -statistic as input data and values of and as input parameters, and outputs a -statistic for every frequency bin in the search range (as per Eqs. 3 and 4).
The extent of the sideband template, Eq. 1, changes as a function of the search frequency since the number of sidebands in the template scales as . We therefore divide the 500 Hz search band into smaller sub-bands over which we can use a single template. The sub-bands must be narrow enough, so that and hence do not change significantly from the lower to the upper edges of the sub-band, and wide enough to contain the entire sideband pattern for each value of . It is preferable to generate -statistic data files matching these sub-bands, so that the search algorithm can call specific -statistic data files for each template, as opposed to each call being directed to the same large data file. However, the -statistic sub-bands need to be half a sideband width (or Hz) wider on each end than the -statistic sub-bands in order to calculate the -statistic at the outer edges. For a Sco X-1 directed search, single-Hz bands are convenient; for example, even up at Hz, the template width is still less than 0.25 Hz.
The output of the -statistic is compared with a threshold value chosen according to a desired false alarm rate (see Sec. IV.3 below). Any frequency bins returning are designated as candidate events and are investigated to determine whether they can be attributed to non-astrophysical origins, due to noise or detector artifacts, or to an astrophysical signal. The former are vetoed and if no candidates above survive, upper limits are computed (see Sec. V.2 for more information on the veto procedure).
iv.3 Detection Threshold
To define the threshold value for a single trial we first relate it to the false alarm probability , i.e. the probability that noise alone would generate a value greater than this threshold. This is given by
where denotes the cumulative distribution function of a distribution evaluated at 999The -statistic is a distributed variable with degrees of freedom because it is the sum over sidebands of the -statistic, which is in turn distributed with 4 degrees of freedom.
In the case of statistically independent trials, the false alarm probability is given by
This can be solved for the detection threshold in the case of trials, giving
where is the inverse (not the reciprocal) of the function .
The search yields a different -statistic for each frequency bin in the search range. If the -statistic values are uncorrelated, we can equate the number of independent trials with the number of independent frequency bins ( for each Hz band). However, due to the comb structure of the signal and template, frequencies separated by an integer number of frequency-modulated sideband spacings become correlated, since each of these values are constructed from sums of -statistic values containing many common values. The pattern of sidebands separated by Hz spans Hz, meaning there are sideband patterns per unit frequency. Hence, as an approximation, it can be assumed that within a single comb template there are independent -statistic results. The number of statistically independent trials per unit Hz is therefore given by the number of independent results in one sideband multiplied by the number of sidebands per unit frequency, i.e.
This is a reduction by a factor in the number of statistically independent -statistic values as compared to the -statistic.
Using this more realistic value of provides a better analytical prediction of the detection threshold for a given , which we can apply to each frequency band in our search. However, a precise determination of significance (taking into account correlations between different C-statistics among other effects) requires a numerical investigation. To factor this in, we use Monte Carlo (MC) simulations to estimate an approximate performance (or loss) factor, denoted . This value is estimated for a handful of 1 Hz wide frequency bands and then applied across the entire search band. For a specific 1 Hz frequency band, the complete search is repeated 100 times, each with a different realisation of Gaussian noise. The maximum obtained from each run is returned, and this distribution of values allows us to estimate the value corresponding to a multi-trial false alarm probability . For each trial frequency band, we can then estimate as
which can be interpreted as the fractional deviation in using the number of degrees of freedom (the expected mean) as the point of reference. We assume that is approximately independent of frequency (supported by Table 3) and hence use a single value to represent the entire band. We incorporate this by defining an updated threshold
which now accounts for approximations in the analysis pipeline. The MC procedure was performed for 1-Hz frequency bands starting at 55, 255, and 555 Hz. Using the values returned from these bands, we take a value of . Table 3 lists the values of , , and associated with each of these bands.