Precision ephemerides for gravitational-wave searches – III. Revised system parameters of Sco X-1
Neutron stars in low-mass X-ray binaries are considered promising candidate sources of continuous gravitational-waves. These neutron stars are typically rotating many hundreds of times a second. The process of accretion can potentially generate and support non-axisymmetric distortions to the compact object, resulting in persistent emission of gravitational-waves. We present a study of existing optical spectroscopic data for Sco X-1, a prime target for continuous gravitational-wave searches, with the aim of providing revised constraints on key orbital parameters required for a directed search with advanced-LIGO data. From a circular orbit fit to an improved radial velocity curve of the Bowen emission components, we derived an updated orbital period and ephemeris. Centre of symmetry measurements from the Bowen Doppler tomogram yield a centre of the disc component of 90 km , which we interpret as a revised upper limit to the projected orbital velocity of the NS . By implementing Monte Carlo binary parameter calculations, and imposing new limits on and the rotational broadening, we obtained a complete set of dynamical system parameter constraints including a new range for of 40–90 km . Finally, we discussed the implications of the updated orbital parameters for future continuous-waves searches.
keywords:ephemerides – gravitational waves – stars: neutron – techniques: radial velocities – X-rays: binaries – X-rays: individual (Sco X-1)
Low-mass X-ray binaries (LMXBs) are interacting binary systems comprising a low-mass donor star transferring matter on to a compact object (either a neutron star or a black hole) via an accretion disc. They can be divided into two populations based on their long-term behavior; those that accrete persistently at high mass-accretion rates ( erg ), and the transient systems which spend most of their lives in a dim quiescent state, but undergo sporadic outbursts when their X-ray luminosity increases. Neutron stars (NSs) in binaries are considered promising sources of gravitational-wave (GW) emission for ground-based GW detectors. Not only are these NSs typically rotating many hundreds of times a second, the process of accretion can potentially generate and support non-axisymmetric distortions to the compact object (e.g. Bildsten 1998). A quadrupole mass moment will lead to persistent emission of GWs at twice the NS spin frequency.
Sco X-1 was the first LMXB discovered outside the Solar System (Giacconi et al., 1962) and the brightest persistent X-ray source in the sky. The orbital period of 18.9 hr was first reported by Gottlieb, Wright & Liller (1975) from photometry of the optical counterpart V818 Sco (V 12.5) using 1068 plates spanning over 84 yr. Bradshaw, Fomalont & Geldzahler (1999) measured the trigonometric parallax from Very Long Baseline Array (VLBA) radio observations, and obtained a distance of 2.8 0.3 kpc. Further observations revealed the presence of twin radio lobes, which led to an estimate for the inclination of , assuming that the radio jet is perpendicular to the orbital plane (Fomalont, Geldzahler & Bradshaw, 2001). It is believed to harbor a NS based on its X-ray behavior. However, no coherent pulsations or bursts have been observed so far from this source.
Sco X-1 has been a prime target for continuous GW searches for the Laser Interferometer Gravitational-Wave Observatory (LIGO). Due to its relative proximity to Earth and high mass accretion rate, it is predicted to be the strongest GW emitter among the family of known LMXBs within the relevant frequency passband (Abbott et al., 2007). A key difficulty in searching for the GWs emitted by NSs in LMXBs is the lack of precise knowledge about the position and velocity of the NS in the binary orbit (Watts et al., 2008). Searching over possible values of unknown parameters carries a cost in a computationally-bound search, and many methods carry a trade-off between computational cost and sensitivity. The smaller the parameter space region that needs to be searched, the more sensitive the search can be made, given fixed computing resources. Without accurate determinations of the orbital parameters, a directed search may still be carried out, but an observational ‘penalty’ must be paid. Effectively, the signal must be proportionately stronger, compared to a source where the parameters are more precisely constrained, in order to reach the same order of confidence for a detection.
The most important orbital parameters required by a directed search include the orbital period (), the absolute phase of the system (for which we quote the epoch of inferior conjunction of the companion star ), and the eccentricity (e.g., Watts et al. 2008). The projected semi-major axis of the NS (sin) must also be searched over (Messenger et al. 2015; Leaci & Prix 2015; Whelan et al. 2015; Meadors et al. 2016; Sammut et al. 2014), and is related to the amplitude of the projected orbital velocity of the NS, (or sin), as = 2 sin/. For systems in active states, the optical emission is dominated by the reprocessing of hard X-rays in the outer accretion disc, which severely hampers the detection of the much fainter companion. A recent study by Mata Sánchez et al. (2015) found no evidence of companion spectral features for Sco X-1 even in the near-infrared (where the companion star could make a substantial contribution). A novel avenue for obtaining dynamical information on Sco X-1 was opened up thanks to the discovery of extremely narrow, high-excitation emission lines arising from the X-ray illuminated atmosphere of the donor star (Steeghs & Casares 2002, hereafter SC02). These narrow components were strongest in the Bowen region (consisting of a blend of \ionNiii/\ionCiii lines), and are the result of fluorescence of the gas by UV photons from the hot inner disc (McClintock, Canizares & Tarter, 1975). Since the first detection of the donor star of Sco X-1 in emission, it has been shown for several other persistent systems that the use of Bowen fluorescence lines as tracers of counterpart radial velocities provides a unique opportunity to pursue robust radial velocity studies (e.g. Casares et al. 2006; Barnes et al. 2007).
Our project – Precision Ephemerides for Gravitational-wave Searches (PEGS) – aims at providing the most reliable constraints for the key orbital parameters for candidate persistent GW sources. Galloway et al. (2014) (hereafter G14) presented the results of a pilot program for improving the precision of the orbital period and ephemeris of the most promising system, Sco X-1, using the proven Bowen diagnostic of the mass donor. In this follow-up study to G14, we provide revised constraints on key orbital parameters (derived from a re-analysis of existing optical spectroscopic data) in direct support of GW observations of Sco X-1 in the Advanced-LIGO era. We also describe new analysis tools that are relevant to the use of emission lines as dynamical probes of the donor stars in compact binaries in general.
2 Source data and RV determination
This study makes use of all available spectroscopic data previously obtained for Sco X-1 and initially presented in SC02 and G14. Two epochs of data separated by 12 years (1999 – 2011) were obtained using the Intermediate dispersion Spectrograph and Imaging System (ISIS) at the 4.2m William Herschel Telescope (WHT). In addition, Sco X-1 was observed at 44 random epochs between 2011 May 29 – August 23 with the Ultraviolet and Visual Echelle Spectrograph (UVES) on the European Southern Observatory (ESO) Very Large Telescope (VLT; programme 087.D-0278). The superior quality of the combined dataset allowed the application of several techniques for obtaining binary parameter constraints, enabling us to go beyond the initial results provided in G14. Table 1 shows a summary of source data. For details of the observations and data reduction procedures, we refer to SC02 and G14.
We first reanalyzed the Bowen emission components from the irradiated front face of the donor using the well-established radial velocity (RV) fitting method. Following SC02 and G14, we used the molly111http://deneb.astro.warwick.ac.uk/phsaap/software mgfit routine for fitting multiple Gaussians to time-resolved, continuum-normalised spectra and (for each Bowen profile) measuring the common RV of the three strongest narrow lines (\ionNiii 4634 & 4640 and \ionCiii 4647). For the WHT data, we used the same multi-Gaussian model (with 6 free parameters) as was used in SC02 and G14. It is known that during certain phase ranges, especially near phase zero, the narrow lines could be much weaker than at other phases. The weakness of the narrow lines had such a dramatic impact on the RV fitting in G14 (see G14; Fig. 5) that the majority of the WHT measurements within 0.1 of phase 0 (i.e. the inferior conjunction of the donor star) were considered spurious and were thus eliminated from further analyses. To salvage some of the data points in this range, we can go one step further and take into account the spectrum-to-spectrum variation. For each spectrum, the expected RV of the narrow lines was calculated (using a sinusoidal model) and used as the initial fit value of RV to guide the optimizer to a reasonable optimum. The quality of individual Gaussian fits was also carefully assessed by visual inspection, especially for spectra taken near the inferior conjunction of the companion.
In order to exploit the potential of our highest resolution VLT data (see, e.g., G14, Fig. 4), we modified the model for extracting UVES velocities. Instead of fixing the FWHM of all the Gaussians, we fixed the width of the underlying broad component to the value used in G14 but left the (common) width of the narrow components free to vary. This yielded a revised set of UVES RVs that are fully consistent with the values reported in G14 (within the formal errors) as well as the fitted width of the sharp \ionNiii/\ionCiii components. We found that the model with 7 free parameters fitted the UVES spectra better, and that the average statistical uncertainty of the VLT measurements was reduced to 1.1 km (compared to the previously reported precision of 2.2 km ).
|Date||Instrument||Exposure time (s)||Phase coverage||Resolution (Å)||Reference|
|1999 Jun 28–30||ISIS/WHT||300||137||75%||0.84||SC02|
|2011 Jun 16–18||ISIS/WHT||300||157||75%||0.32||G14|
|2011 May 29–Aug 23||UVES/VLT||720||44||34%||0.1||G14|
3 Determination of observables
The number of unreliable RVs that deviate significantly from the predicted values after applying the modified approach is substantially lower than previously estimated. Narrow lines from the donor are still present and can be detectable within the phase ranges 0.1 to 0.1 (albeit being rather weak since we are facing the un-irradiated side of the donor). The derived RVs within those ranges have noticeably larger error bars than data points at other phases but ought to be included in the RV fit (with correspondingly smaller weights). For these reasons, we refrain from screening data (as was introduced in G14 to eliminate a large fraction of spurious RVs), and choose to keep the total of 338 WHT & VLT measurements in the RV analyses of this paper to derive a revised set of Sco X-1 system parameters. The revised and refined RV datasets are provided in the supporting information available online.
3.1 Circular orbit fit
Assuming a negligible (intrinsic) eccentricity of the binary orbit ( = 0), the simplest model for the variation of the projected Bowen RVs is a sinusoid of the form
where K is the velocity semi-amplitude, is the epoch of inferior conjunction of the mass donor, is the orbital period, and is the RV of the system (systemic velocity).
To fit the 4 free parameters (, , K, and ), we started with a least-squares method where each data point was weighted according to the associated error bar given by mgfit whilst giving no differentiation between separate VLT & WHT datasets. The reduced value we obtained ( = 11.3 for 334 degrees of freedom) was greater than that obtained in G14 ( = 8.5 for 238 degrees of freedom) due to the inclusion of all the RV measurements, and much greater than unity. A high value of the -statistic can either be caused by limitation of the sine-wave model, or reflect the fact that the initial measurement uncertainties on the RVs are underestimated. The fit residuals we see are not consistent with those expected from irradiation effects (see Section 4). Therefore, for the purpose of parameter estimation, we considered that the measurement errors had been systematically underestimated by some fractional amount (efac).
Hence we proceeded to introduce additional ‘multiplicative factors’ into the sinusoidal model for adjusting errors on the RV data. To recognize that the data consists of three datasets with different spectral resolution/SNR contributions (1999 WHT, 2011 WHT and 2011 VLT; see Table 1), we invoked 3 multiplicative factors in total to rescale errors on individual multi-Gaussian fits per dataset. The 1999 and 2011 WHT observations (both covering 75% of the orbital cycle) were then refitted separately with the Bayesian Markov Chain Monte Carlo (MCMC) approach (using the python package emcee; Foreman-Mackey et al. 2013). In both cases, the resulting posterior distribution of the multiplicative factor was Gaussian, therefore we took the posterior mean ( = 4.3 and = 2.5) as the best estimate of efac for the WHT datasets. It is also worth noting that the RV semi-amplitude for the 1999 data (77.6 1.3 km ; 68 per cent) is significantly higher than the 2011 velocity (72.0 1.0 km ; 68 per cent). This finding is consistent with the discrepancy between the best-fit velocity semi-amplitude for the combined data ( = 74.9 0.5; G14) and the result of SC02 ( = 77.2 0.4 km ) based on the 1999 observations alone; and might be attributed to variations in the geometry of the line emission over epoch, as discussed in G14.
To determine conservative parameter estimates, we scaled the WHT errors (up) by the efac values derived from the MCMC fit and repeated the least-squares analysis of the combined WHT and VLT data. At the final step, the VLT errors were scaled by a multiplicative factor ( = 2.25). This ensures that the final reduced is equal to 1.0. We also varied the initial parameter value of to obtain the smallest possible cross-term in the covariance matrix between and ((, ) = 3 ). This yielded revised orbital parameters consistent (to within 2.6) with those previously determined (by least-squares fitting to 242 RV measurements) in G14. Fig. 1 shows the improved RV curve (with the best fit sinusoid) consisting of a total of 338 measurements (cf. G14; Fig. 5).
The revised parameter estimates are summarized as follows for direct comparison with G14:
= 2455522.6668 0.0006 HJD(UTC)
= 0.7873132 0.0000005 days
= -113.6 0.2 km
= 74.8 0.4 km .
To compare our new constraint with that of G14, we added 1127 orbital cycles to give 2455522.6682 HJD (UTC), which differs from the estimate in this work by 121 seconds, within the 1 uncertainty of G14’s estimate propagated to 2010. The new measurement improves the constraint on the absolute phase of the system by a factor of two (both by being more precise and by being less correlated with when considering times after 2010) relative to the result of G14. On the other hand, the measurement is at the same precision as that in G14, but differs from that value by 2.6. In Section 3.2.2, we perform tests to show that the new period and ephemeris constraints (derived from refined RV datasets) are indeed more accurate than previous findings.
3.2 Doppler tomography-based method
An alternative route towards orbital parameter estimation is provided by a technique known as Doppler tomography, developed by Marsh & Horne (1988) for recovering accretion structures in mass-transferring binary systems. The idea is that by observing time-dependent line profiles as a function of the binary orbit (preferably covering a complete orbital cycle), sufficient information can be obtained to invert the phase-resolved spectra into an ‘image’ of emission distribution in Doppler/velocity coordinates. In a Doppler tomogram, the Roche-lobe-filling donor star is on the positive -axis, shifted from the origin (, ) = (0, 0) by its projected RV semi-amplitude (); the compact accretor (although often unseen) is below the origin at (0, -); the accretion disc still appears circular in velocity space, but is inside-out (i.e., the outer edge of the disc is nearer the origin). Other major emission sites might be traced at an expected position, including the gas stream leaving the secondary and the stream-disc impact region. For a review of the technique, see Marsh (2001); Steeghs (2004).
Since Doppler tomography uses all observed data simultaneously, it is ideally suited for weak features and provides the only way to detect the donor signature in the cases of low signal-to-noise ratios (e.g. Cornelisse et al. 2008; Wang et al. 2017). The technique has also been applied to the high SNR data of Sco X-1 to reconstruct maps for a range of emission lines present in the optical spectrum (see SC02; Fig. 4). The \ionNiii and \ionCiii maps in particular revealed a compact spot along the positive -axis (consistent with an emission component from the irradiated donor), from which an independent estimate of could be accurately derived through a two-dimensional Gaussian fit.
3.2.1 The Bowen map
To enable a direct comparison to be made between the amplitude obtained from the RV fitting and the tomography-based method, we constructed new Bowen Doppler images (initially separately for the 1999 and 2011 WHT datasets; each covering 75% of the orbital cycle) by assigning the three strongest Bowen emission lines to the same image, and updating the input parameters to Doppler tomography (, , and ) with the values determined in Section 3.1. Both 1999 and 2011 WHT Bowen maps showed a prominent spot at the expected location of the companion, and a much fainter ring-like feature, most likely originating from the accretion disc. We note again that the centroid position of the 1999 donor emission spot is at a higher velocity ( = 2 km and = 76 km ) than that inferred from the 2011 donor spot ( = 72.2 km ). Furthermore, with the aid of the second-generation doppler code222https://github.com/trmrsh/trm-doppler developed by T. Marsh (e.g., Manser et al. 2016), it is possible to construct one Bowen image for the combined 1999 and 2011 WHT data even though these have different spectral resolutions. The combined map and the phase-binned, trailed spectrograms of the Bowen region (from the WHT observations) are displayed in Fig. 2.
3.2.2 Bootstrap Doppler tomography
For robust estimation of from the combined Bowen map we adopted the ‘bootstrap Monte Carlo’ approach recently introduced by Wang et al. (2017) (hereafter, W17). The method was initially developed for statistical uncertainty derivation and, most importantly, robust significance testing for spot features in the low SNR case of XTE J1814338; but can be readily extended to the mid- to high-SNR regime. Following the strategy described in W17, we created 2000 simulated images from 2000 independently generated bootstrapped datasets (Efron & Tibshirani, 1993) in the same manner as for the original map (see W17, Section 5.2.2), using custom wrapper functions. Measurements of various spot properties (e.g., the centroid position and the peak intensity) were then performed for the ensemble of combined images.
Fig. 3 shows the histograms of the most relevant spot parameters as derived from the bootstrap samples. Firstly, we considered the distribution of the peak height in terms of a significance level of the (combined) donor signal. By fitting a Gaussian curve to the histogram plot in the left panel of Fig. 3, we deduce that the mean of the peak intensity distribution is above zero at the >58 level. Moreover, in contrast to the results for the much lower SNR data of XTE J1814338 (which indicated a 4 detection of the mass donor), Fig. 3 (middle panel) shows that the combined Bowen image for Sco X-1 provides a tight constraint on the RV semi-amplitude, = 75.0 0.8 km , for an assumed value of = -113.6 km .
It is known that the choice of the input parameter can have an impact on the resulting reconstruction – and spot features in particular would appear out-of-focus and misplaced if an incorrect was used to obtain the image (e.g. Steeghs 2003; Muñoz-Darias et al. 2009). Thus, we ran the same analysis using a range of values in steps of 2 km , centered on 113.6 km . Within a narrow range of values between 120 and 108 km , a >50 donor detection was achieved, peaking near the expected value from our RV fit. The derived velocity showed a small range of variation (with a maximum of 0.6 km deviation from 75 km ), therefore indicating the presence of a small systematic error on due to the uncertainty in . We note that this constraint on the systemic velocity ( km ) is looser compared to the best-fit value derived from the RVs ( km ) since it is not very sensitive to changes below the spectral resolution of the data. Based on the above results, we obtain
= 75.0 0.8 (stat) 0.6 (sys) km ,
in excellent agreement with the value derived from the traditional RV fitting method. Therefore the novel technique of bootstrap Monte Carlo delivers a robust and more conservative estimate of for the combined observations, which we adopt as our best estimate of the Bowen amplitude. Since Doppler tomography exploits all spectra at once, the technique is also less affected by the problem of extremely weak features near phase zero (see G14; Fig. 2) compared to individual line profile fitting. This makes bootstrap Doppler tomography an attractive tool for deriving binary system parameters across all SNRs. The main disadvantage of the tomography-based method is that and are used as input, instead of being determined directly from the reconstruction. However, the phase shift () between the donor spot and the vertical plane can provide a good indication of the accuracy of the input ephemeris.
With the correct ephemeris, emission from the donor should always be projected on the tomogram at the position = 0 and = (hence = 0), whereas phase errors tend to lead to a rotation of the image. In the right panel of Fig. 3 we show that is consistent with zero (blue solid line) within 2 for our best estimate of . A diagnostic plot of as a function of the assumed (between -120 and -108 km ) gives = 0.003 0.002 (statistical) 0.002 (systematic). It can thus be concluded that the donor spot is indeed centered on the positive -axis using the revised binary period and ephemeris. Table 2 lists the final adopted values of , , , and the cross-term (needed to determine the propagated uncertainty on the future epoch of inferior conjunction; see G14).
A similar test using and of G14 yielded a significant (3) phase shift between the donor spot and the vertical plane, which indicates that the previous results included contributions from systematic errors. Since the revised period and ephemeris were derived from an improved RV curve (including all data points) using a refined analysis approach, and there is no measurable phase shift, the new estimates for orbital parameters presented in this paper can be considered more reliable than the results of G14.
4 Estimation of binary parameters
So far we have exploited the proven Bowen emission line diagnostic that permits a robust and accurate RV study of the companion star in Sco X-1. The main caveat is that the narrow lines are produced on the irradiated inner hemisphere of the secondary, instead of the entire Roche lobe; hence the previously determined Bowen amplitude represents only a lower limit to the true, centre-of-mass velocity of the companion, . Therefore a ‘K-correction’ must be applied to the observed when performing binary parameter calculations.
To precisely constrain the remaining parameters for the purpose of GW searches, and , requires a direct determination of the neutron star orbit, which is particularly challenging due to the lack of NS pulsations. While coherent pulsations of Sco X-1 are not observed, we shall rely on the known set of parameter constraints and apply Kepler’s laws to derive a more conservative estimate compared to that of SC02 obtained using broad disc emission features.
4.1 The K-correction
The deviation between the reprocessed light centre and the centre of mass of the Roche-lobe-filling donor ( = 1) can be modelled by generating RV curves from simulated emission line profiles formed by isotropic irradiation on the illuminated hemisphere. Results of the most extensive study of K-correction modelling to date (Muñoz-Darias et al. 2005; hereafter, MCM05) suggest a weak dependence of on the orbital inclination (); but a strong dependence on the mass ratio between the donor star and the compact object ( = /) and the disc opening angle (), which represents the shielding effect by a flared and opaque accretion disc. By fitting fourth-order polynomial models to the synthetic K-corrections (as functions of ), MCM05 derived parametric approximations to = , for different disc flaring angles between and , in both low ( = ) and high inclination ( = ) cases (see Fig. 4 and Table 1 of their paper).
For Sco X-1, an estimate of the binary inclination is available from modelling of the orientation of radio lobes ( = ; Fomalont et al. 2001). However, the disc shielding parameter is unknown and therefore throughout the rest of our analysis, we consider the full range of possibilities – from the extreme case of neglecting disc shadowing ( = ; maximum displacement) to = where the companion is nearly completely shadowed by the accretion disc (minimum displacement). The key ingredient for improving the K-correction then involves mainly the determination of constraints on .
4.1.1 Center of symmetry search
The ring-like disc emission feature revealed using Doppler tomography (SC02 and this paper) may be used to constrain the NS’s projected RV (). In Doppler coordinates, the projection of an ideal Keplerian disc is always centered on the expected position of the NS (, ) = (0, -), as the disc gas should track the orbital motion of the accreting primary. Thus, although the accretor is rarely seen in emission, a search for the centre of symmetry (CoS) of the disc component in velocity space should, in principle, yield a valid solution for (e.g. van Spaandonk et al. 2010; Longa-Peña et al. 2015).
For a reconstructed Doppler map, one first subtracts a symmetric component around an assumed centre ( = 0, = -), followed by inspecting the resulting residual image for any remaining asymmetries not related to the primary, e.g., emission caused by the stream/disc impact. The process is then repeated for a wide range of values to select an optimal center that minimizes the residuals over a region of the map containing clean disc emission. Since non-disc components and hot spot contamination can be easily masked during the minimization, the method has the advantage of allowing one to remove the distortions by any asymmetries which would normally disrupt the outer disc, and focus on the higher velocity emission arising from the (undisturbed) inner part of the disc (appearing further away from the origin in velocity space).
Under the assumption that the disc is axisymmetric around the NS primary, the technique was adopted in SC02 to obtain an initial estimate for Sco X-1 based on the location of the CoS for the H image ( = 40 km ). Here we focus on applying the method to the newly computed Bowen image – which reveals the best defined accretion disc among the maps reconstructed from principal spectral features (see, e.g., SC02; Fig. 4) – in order to assess the stability and the validity of the CoS measurement.
We cycled the assumed -coordinate of the disc centre between -200 and +20 km in steps of 2 km (while keeping fixed to 0) and measured the mean residuals in the lower half of the reconstruction (to exclude the region near the secondary emission). For the 1999 Bowen data the residuals showed a distinct minimum at = - = -84 km . A similar result was obtained also for the 2011 Bowen image ( = 90 km ), suggesting that a significantly higher value of amplitude than previously estimated (using the same method but a different line feature) needs to be considered. Given the absence of a sharp disc ring in the 1999 H map, compounded by the presence of enhanced emission in the bottom left quadrant, it is most likely that the best fit value still failed to represent the NS’s true projected RV. These significant differences between different lines show that extracting robust constraints from disk lines remains challenging. Hence we decided to adopt the greatest value for based on CoS measurements using maps of a range of emission lines as our most conservative maximum value of . Using a new upper limit on
90 km ,
and apply the K-correction (for = ) to , we obtain an upper limit of 0.76.
4.1.2 Rotational line broadening
Assuming rotation to be the dominant line broadening mechanism, the observed width of the sharp \ionNiii/\ionCiii lines ( sin) to can be used to set a strict lower limit on . For a spectral feature that is produced throughout the Roche lobe, the expected width ( sini) is given by the relation,
(Wade & Horne, 1988). Since the donor is not fully illuminated, this becomes a limit as sini > sin. By taking the weighted average of the set of 12 FWHM (\ionNiii/\ionCiii) values measured from the VLT spectra (see Section 2) near orbital phase 0.5 when the visibility of the irradiated face is maximum, we obtain an estimate for the width of the narrow components of FWHM = 51.9 0.5 km . However, it should be cautioned that this preliminary estimate does not provide a strict lower limit to sini as the measured widths include also the intrinsic broadening effect due to our instrumental resolution (8 km s).
Following Casares et al. (2006), we corrected for the instrumental effect by broadening a Gaussian template of FWHM = 8 km s using a Gray rotational profile (Gray, 1992) without limb-darkening (because fluorescence lines arise in optically thin conditions); and subsequently determined for each of the 12 VLT measurements the amount of rotational broadening required to reproduce the apparent FWHM (\ionNiii/\ionCiii). This leads to a weighted average of sin = 37.8 0.4 km , which combined with the K-correction for = , gives a secure lower limit of 0.22.
The lower limit of can be further improved by applying a correction factor (defined as = 1) to convert the observed rotational broadening to true sini, for the reasons mentioned above. In light of the results of the MCM05 study of , we estimate this correction factor in the low inclination case ( = ) – suitable to be applied to the case of Sco X-1 – by following the same approach for modelling the K-correction described in MCM05. The details of the modelling are described in the appendix. The solution for at = represents the most conservative correction (see appendix), which can be safely adopted to yield an improved lower limit of (at a value slightly higher than 0.22, Section 4.2), and thereby also providing the best lower limit to (see Section 4.2.1).
4.2 Monte-Carlo analysis
The discovery of the Bowen emission lines (made more than a decade ago) arising from the donor star in Sco X-1, combined with the inclination derived by Fomalont et al. (2001), led to some of the first dynamical parameter constraints for this prototypical LMXB (SC02). In this section, armed with the complete set of K-correction models solved by MCM05 and new limits on and the rotational broadening (see Section 4.1), we perform Monte Carlo simulations for deriving the probability density functions (PDF’s) of Sco X-1 system parameters. Of particular note is the distribution (or sin = /(2) ) required by most directed searches for continuous GWs.
As an initial, conservative estimate, we used the binary inclination inferred from the orientation of radio jets and results from the reanalysis of the combined Bowen data (Section 3 only). Synthetic values of , and were selected from a Gauss-normal distribution (with the mean and standard deviation equal to the values listed in Table 2). To perform the K-correction, we started with the most conservative scenario by simulating uniformly-distributed values of and over the entire ranges considered by MCM05 (covering values typical of XRBs) between and ; 0.05 and 0.83. For each pair of randomly generated disc opening angle and mass ratio, we determined a precise value of by interpolating over the grid of MCM05 models. After applying to , we obtained the corresponding amplitude, which was then used (along with , and ) as the input to the mass function equation
for deriving the mass of the compact primary . At the end of each trial, the RV of the primary = and the mass of the donor were also calculated; and the process was repeated times. Assuming a NS accretor, only the outcomes of trials that yielded an value between the minimum ( 0.9 ; e.g., Lattimer 2012) and the maximum stable NS mass (3.2 ) were assembled into the initial, conservative estimate of the PDF’s (Figs 4 and 5; grey).
4.2.1 Radial velocity of the NS
As a next step we imposed an upper limit on determined from center of symmetry analyses of the disc component ( 90 km ), which essentially sets a more realistic upper limit of as discussed in Section 4.1.1. For each Monte Carlo trial, we additionally computed the minimum allowable value by applying the most conservative -correction (derived in the appendix). This is the case of = , dependent only on , and considering the lower limit on 36.6 km s (99.7 per cent confidence) as estimated from the width of the sharp \ionNiii/\ionCiii lines. A conservative and stringent lower limit of could therefore be established using the relationship between and (equation 2). Finally, the output PDF’s that reflect our current best knowledge of Sco X-1 parameters are shown in Figs 4 and 5 (green), with those derived only from the observed distributions of and .
Because of these newly imposed constraints, the CIs for nearly all system parameters (except ) are narrowed down significantly. In particular, the estimate of the centre of the disc component had such a strong impact that the upper limit was reduced from 165 to 90 km ; and thus represents the most valuable determination in the context of refining the parameter space for GW searches. Ruling out high values of (or sin) has a double benefit, since high values of sin increase the required resolution in orbital phase (). We also note that the new sini constraint (combined with a minimum -correction) provides a refined lower limit not only for but also for both and , due to the dependence of on . This yielded the current best hard limits on of 40–90 km , or equivalently,
1.45 ls sin 3.25 ls.
This range of is significantly broader than the previous constraint (40 5 km ) derived solely from centre of symmetry measurements from the H Doppler map (SC02). In the recent search for continuous gravitational-waves from Sco X-1 (e.g., Abbott et al. 2017b), the range for the projected semi-major axis has been expanded in order to cover the full parameter space.
4.2.2 Determination of the component masses
Despite the fact that the RV semi-amplitude of the Bowen emission can be measured with a precision of 1 per cent, the true projected RV of the donor is not well constrained ( = km ; 95 per cent confidence) given the lack of knowledge of and the disc shielding parameter. Using a geometric model for the reprocessing of X-rays in LMXBs, an average disc opening angle of has been proposed for other persistent sources, e.g., GR Mus and V4134 Sgr (de Jong et al., 1996). Similarly, for Sco X-1, an -value at the lowest end of the range between and is perhaps unlikely (given also the clear presence of an accretion disc). However, we chose to consider all values within the most conservative range between 0– in the Monte Carlo analyses to avoid introducing biases into the resulting PDF’s.
Due to the large uncertainties in the final values of (dominated by the unknown K-correction) and (; 95 per cent), it remains difficult to put strong constraints on the component masses. Using the inclination value of – inferred from the orientation of radio jets (assuming that the radio jet is perpendicular to the orbital plane) – we obtain a 95 per cent credible range of the NS mass between 0.9–2.8 and an range of 0.4–1.5 (95 per cent). Interestingly though, the median of our best estimate for the PDF (after applying the and constraints) coincides with the ‘canonical’ value of 1.4 – in support of the presence of a 1.4 neutron star.
The complete set of dynamical system parameter constraints, derived from the latest RV measurements as well as Monte-Carlo simulations, are provided in Table 2.
|(km s)||-113.6 0.2|
|(km s)||75.0 0.8 (statistical) 0.6 (systematic)|
Notes. The range for the projected semi-major axis of the NS sin = /(2) in light-seconds; Adapted from (Fomalont et al., 2001).
4.3 The eccentricity problem
One complication of the Bowen technique is that the established RV curve is not expected to be perfectly sinusoidal owing to the effect of non-uniform flux distribution at the surface of the irradiated donor star. This is a particular handicap for attempts to constrain any true orbital eccentricity,, of the binary. The influence of irradiation effects can be clearly seen in Fig. 6, where we display the simulated residual RV curves obtained by subtracting the sinusoidal fit from the irradiation model (described in the appendix) for the most probable mass ratio of Sco X-1 ( 0.5; as calculated by the Monte Carlo method) and different -values. As noted in the appendix and in previous studies, the amount of deviation of the simulated RVs from a sinusoidal wave form increases with decreasing disc opening angle (see also MCM05; Fig. 2). Closer inspections reveal that for all values the deviations increase at small and large orbital phases ( 0.2 and 0.8), consistent with the distortions as seen, for example, in the observed RV curves of the pre-cataclysmic binary NN Ser (see Parsons et al. 2010; Fig. 8). Based on results from both observations and simulations, we conclude that RV distortions caused by irradiation effects would dominate those produced by any small intrinsic of the binary orbit. Although they are not described by an eccentric model, attempting to fit such RV curves with an eccentric orbit would inevitably produce an apparent eccentricity (E), which varies strongly with disc opening angle333As we mentioned in the appendix, the amount of deviation of RVs from a pure sinusoid increases with increasing mass ratio, therefore E varies also with ..
Worse still, Fig. 6 also shows that the WHT data do not yet offer the RV precision high enough for differentiating between the irradiation and the circular orbit model. Not even the 2011 UVES/VLT data has the precision high enough to firmly differentiate between residual RVs for low and high angles, thus making it currently infeasible to accurately infer the value of the disc opening angle necessary to clearly disentangle E from any true orbital eccentricity.
Due to the inherent limitations of the method, as well as limitations of the current data sets, in this paper we only attempt to place a crude upper bound on by running an MCMC analysis (to fit the combined RV data) and assuming an elliptical orbit model. By leaving the spurious eccentricity () and periastron angle as free parameters, we obtained a 95 per cent credible range of the eccentricity parameter between 0.00086-0.082 (using a uniform prior in log ). However, we stress that the nominal upper bound of 0.082, determined without taking account of any irradiation effects (note that E/ 1), cannot accurately represent the limit to the true orbital , which is expected to be much closer to zero for a Roche-lobe overflow XRB with a high mass-accretion rate (e.g., Sco X-1).
5 Conclusions and Future Perspectives
We presented a reanalysis of the combined spectroscopic data of SC02 and G14 with the aim of providing revised constraints on key orbital parameters of Sco X-1. Here we paid particular attention to the requirements of directed searches for continuous GWs. By exploiting the previously established spectroscopic diagnostic of the mass donor, we derived an improved RV curve of the Bowen blend, arising from the X-ray heated hemisphere of the companion, and demonstrated that both the binary period and ephemeris can be precisely determined using the traditional RV fitting method. Additionally, Doppler tomography of the Bowen emission region clearly reveals a donor spot feature along the positive -axis, providing a viable alternative to reliably measure the RV semi-amplitude (< ). In the absence of coherent X-ray pulsations, a precise measurement of the projected RV semi-amplitude of the NS is still challenging.
The latest centre of symmetry measurements (from the Bowen Doppler tomogram) yielded a CoS of a sharp disc component of 90 km , which we interpret as a revised upper limit to (larger than previously reported by SC02). We explore the implications of these by implementing Monte Carlo binary parameter calculations, combined with the K-correction (to transform to ) and the rotational broadening . By considering a random distribution of mass ratios and disc opening angles, we obtain a looser constraint of (40–90 km , or equivalently, 1.45 ls sin 3.25) than that derived in SC02 (40 5 km ). Finally, we discussed the pitfalls of measuring the eccentricity of the NS orbit; namely, irradiation effects are known to produce an apparent eccentricity, which is hard to disentangle from the true orbital using the current Bowen dataset. By simply allowing for an ellipticity in the fit to the RV curve (without taking account of any irradiation effects), a nominal upper limit for the eccentricity < 0.082 was obtained. However, we stress that the true is expected to be much closer to zero in the case of Sco X-1, and has therefore been ignored in recent directed searches with (advanced) LIGO data (e.g. Abbott et al. 2017a; Abbott et al. 2017b).
In light of our new constraints on orbital parameters, the ranges of search parameters and should be updated with the refined period and ephemeris. More importantly, the range for the projected semi-major axis sin needs to be expanded and updated with the new confidence interval for , which reflects the systematic uncertainty much greater in magnitude than the statistical error estimated in SC02. With the previous values (which spanned an unrealistically narrow range), there was a high probability of missing the GW signal. In the first Advanced LIGO observing run (O1) searches, a preliminary constraint of 10 km 90 km (available at the time the search was constructed) was adopted as the prior assumption. The best marginalized 95% upper limit on the signal amplitude (obtained by Abbott et al. 2017b) reaches 2.3 , which is a factor of 7 stronger than the best upper limits set using data from initial LIGO science runs. In future runs, computing resources will be concentrated on the final refined range of 40 km 90 km , making the search cheaper than in the O1 search, which allows a more sensitive search at the same computing cost.
In the future, we will combine RV measurements provided by our ongoing monitoring campaign of UVES/VLT observations with the 1999 and 2011 measurements to further refine the binary period and ephemeris and ensure a good ephemeris is available throughout the Advanced-LIGO era (see also G14). Superior RV precision derived from future UVES data may also allow us to rule out low values of (see Section 4.2.2) by matching the observed residual RVs to the simulated residuals. This effort would directly improve the K-correction, and thereby improve the constraints on and . Additionally, we will seek to exploit the more complete phase coverage of the UVES/VLT data to compute Doppler tomograms for the VLT spectra, which may provide further insights into the value of . The determination of will remain difficult due to our poor knowledge of both and .
DS and TRM acknowledge support from STFC via grant ST/P000495/1. This project was supported in part by the Monash-Warwick Strategic Funding Initiative. D.K.G. acknowledges the support of the Australian Research Council via the Future Fellowship scheme (project FT0991598). We thank John Whelan and the rest of the LSC continuous-waves search group for useful discussions which helped shape this manuscript. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 087.D-0278.
- Abbott et al. (2007) Abbott B., et al., 2007, Phys. Rev. D, 76, 082001
- Abbott et al. (2017a) Abbott B. P., et al., 2017a, Phys. Rev. D, 95, 122003
- Abbott et al. (2017b) Abbott B. P., et al., 2017b, ApJ, 847, 47
- Barnes et al. (2007) Barnes A. D., Casares J., Cornelisse R., Charles P. A., Steeghs D., Hynes R. I., O’Brien K., 2007, MNRAS, 380, 1182
- Bildsten (1998) Bildsten L., 1998, ApJ, 501, L89
- Bradshaw et al. (1999) Bradshaw C. F., Fomalont E. B., Geldzahler B. J., 1999, ApJ, 512, L121
- Casares et al. (2006) Casares J., Cornelisse R., Steeghs D., Charles P. A., Hynes R. I., O’Brien K., Strohmayer T. E., 2006, MNRAS, 373, 1235
- Cornelisse et al. (2008) Cornelisse R., Casares J., Muñoz-Darias T., Steeghs D., Charles P., Hynes R., O’Brien K., Barnes A., 2008, in Bandyopadhyay R. M., Wachter S., Gelino D., Gelino C. R., eds, American Institute of Physics Conference Series Vol. 1010, A Population Explosion: The Nature & Evolution of X-ray Binaries in Diverse Environments. pp 148–152 (arXiv:0801.3367), doi:10.1063/1.2945024
- Efron & Tibshirani (1993) Efron B., Tibshirani R. J., 1993, An Introduction to the Bootstrap. Chapman & Hall, New York, NY
- Fomalont et al. (2001) Fomalont E. B., Geldzahler B. J., Bradshaw C. F., 2001, ApJ, 558, 283
- Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
- Galloway et al. (2014) Galloway D. K., Premachandra S., Steeghs D., Marsh T., Casares J., Cornelisse R., 2014, ApJ, 781, 14
- Giacconi et al. (1962) Giacconi R., Gursky H., Paolini F. R., Rossi B. B., 1962, Physical Review Letters, 9, 439
- Gottlieb et al. (1975) Gottlieb E. W., Wright E. L., Liller W., 1975, ApJ, 195, L33
- Gray (1992) Gray D. F., 1992, The observation and analysis of stellar photospheres.
- Lattimer (2012) Lattimer J. M., 2012, Annual Review of Nuclear and Particle Science, 62, 485
- Leaci & Prix (2015) Leaci P., Prix R., 2015, Phys. Rev. D, 91, 102003
- Longa-Peña et al. (2015) Longa-Peña P., Steeghs D., Marsh T., 2015, MNRAS, 447, 149
- Manser et al. (2016) Manser C. J., et al., 2016, MNRAS, 455, 4467
- Marsh (2001) Marsh T. R., 2001, in Boffin H. M. J., Steeghs D., Cuypers J., eds, Lecture Notes in Physics, Berlin Springer Verlag Vol. 573, Astrotomography, Indirect Imaging Methods in Observational Astronomy. p. 1 (arXiv:astro-ph/0011020)
- Marsh & Horne (1988) Marsh T. R., Horne K., 1988, MNRAS, 235, 269
- Mata Sánchez et al. (2015) Mata Sánchez D., Muñoz-Darias T., Casares J., Steeghs D., Ramos Almeida C., Acosta Pulido J. A., 2015, MNRAS, 449, L1
- McClintock et al. (1975) McClintock J. E., Canizares C. R., Tarter C. B., 1975, ApJ, 198, 641
- Meadors et al. (2016) Meadors G. D., Goetz E., Riles K., 2016, Classical and Quantum Gravity, 33, 105017
- Messenger et al. (2015) Messenger C., et al., 2015, Phys. Rev. D, 92, 023006
- Muñoz-Darias et al. (2005) Muñoz-Darias T., Casares J., Martínez-Pais I. G., 2005, ApJ, 635, 502
- Muñoz-Darias et al. (2009) Muñoz-Darias T., Casares J., O’Brien K., Steeghs D., Martínez-Pais I. G., Cornelisse R., Charles P. A., 2009, MNRAS, 394, L136
- Parsons et al. (2010) Parsons S. G., Marsh T. R., Copperwheat C. M., Dhillon V. S., Littlefair S. P., Gänsicke B. T., Hickman R., 2010, MNRAS, 402, 2591
- Sammut et al. (2014) Sammut L., Messenger C., Melatos A., Owen B. J., 2014, Phys. Rev. D, 89, 043001
- Southworth et al. (2009) Southworth J., Hickman R. D. G., Marsh T. R., Rebassa-Mansergas A., Gänsicke B. T., Copperwheat C. M., Rodríguez-Gil P., 2009, A&A, 507, 929
- Steeghs (2003) Steeghs D., 2003, MNRAS, 344, 448
- Steeghs (2004) Steeghs D., 2004, Astronomische Nachrichten, 325, 185
- Steeghs & Casares (2002) Steeghs D., Casares J., 2002, ApJ, 568, 273
- Wade & Horne (1988) Wade R. A., Horne K., 1988, ApJ, 324, 411
- Wang et al. (2017) Wang L., Steeghs D., Casares J., Charles P. A., Muñoz-Darias T., Marsh T. R., Hynes R. I., O’Brien K., 2017, MNRAS, 466, 2261
- Watts et al. (2008) Watts A. L., Krishnan B., Bildsten L., Schutz B. F., 2008, MNRAS, 389, 839
- Whelan et al. (2015) Whelan J. T., Sundaresan S., Zhang Y., Peiris P., 2015, Phys. Rev. D, 91, 102005
- de Jong et al. (1996) de Jong J. A., van Paradijs J., Augusteijn T., 1996, A&A, 314, 484
- van Spaandonk et al. (2010) van Spaandonk L., Steeghs D., Marsh T. R., Torres M. A. P., 2010, MNRAS, 401, 1857
Appendix A Irradiation modelling
Using the lprofile code (within the lcurve package) provided by T. Marsh (e.g. Southworth et al. 2009; Parsons et al. 2010), we started by generating phase-resolved synthetic line profiles from an irradiated donor for all possible -values and different mass ratios between the newly determined range of of 0.22–0.76. In order to derive the most appropriate factors, limb darkening effects were not included in the simulations for the optically thin regime examined here. Other key input parameters representing the instrumental setup (e.g., the exposure length and the FWHM spectral resolution) were fixed at values that match our UVES/VLT observations. Next, we fitted a 1D Gaussian profile to the model data sets and measured the widths of the emission components as a function of orbital phase, which always exhibit a single peak at phase 0.5 for all possible combinations of and . Therefore, we take the width measurement at phase 0.5 (after correction for instrumental broadening) as the estimate for the observed that leads to the largest correction factor (i.e. smallest deviation). Since both and true are known, the true sini can be calculated using equation (2) and hence the ratio of can be determined for each irradiation model.
Meanwhile, by establishing synthetic RV curves (through 1D Gaussian fits to the line profiles), it can clearly be seen that the RVs of the donor emission may deviate significantly from a sinusoid at high mass ratios or low disc opening angles (more details will be discussed later in Section 4.3). We reconstructed synthetic Doppler maps to extract (by computing the centroid of the donor emission spot; see Section 3.2), from which we could independently estimate = , to offer a comparison to MCM05.
Fig. 7 presents example results of our simulations of (black) and (blue) versus (for the cases of = and ), with the K-correction of MCM05 produced by a different binary code overplotted (as blue dashed lines) for comparison. One can see that both factors are always significantly below unity, and this can be exploited to derive firmer limits. The K-correction decreases smoothly with both and , and has already been approximated by a set of fourth-order polynomials by MCM05 for from to in steps of . For small disc opening angles ( ), our solutions for coincide almost exactly with the numerical models of MCM05 (as shown in Fig. 7). The results were also in good agreement (within 3.5 per cent) for . The -correction factor, on the other hand, tends to decrease as increases; yet at any given -angle the estimated value of does not vary smoothly with , fluctuating increasingly randomly with increasing .
Supplementary data are available at MNRAS online.