Abstract
We present the results of a model-based search for continuous gravitational waves from the low-mass X-ray binary Scorpius X-1 using LIGO detector data from the third observing run of Advanced LIGO and Advanced Virgo. This is a semicoherent search that uses details of the signal model to coherently combine data separated by less than a specified coherence time, which can be adjusted to balance sensitivity with computing cost. The search covered a range of gravitational-wave frequencies from 25 to 1600 Hz, as well as ranges in orbital speed, frequency, and phase determined from observational constraints. No significant detection candidates were found, and upper limits were set as a function of frequency. The most stringent limits, between 100 and 200 Hz, correspond to an amplitude h0 of about 10−25 when marginalized isotropically over the unknown inclination angle of the neutron star's rotation axis, or less than 4 × 10−26 assuming the optimal orientation. The sensitivity of this search is now probing amplitudes predicted by models of torque balance equilibrium. For the usual conservative model assuming accretion at the surface of the neutron star, our isotropically marginalized upper limits are close to the predicted amplitude from about 70 to 100 Hz; the limits assuming that the neutron star spin is aligned with the most likely orbital angular momentum are below the conservative torque balance predictions from 40 to 200 Hz. Assuming a broader range of accretion models, our direct limits on gravitational-wave amplitude delve into the relevant parameter space over a wide range of frequencies, to 500 Hz or more.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
Rapidly rotating neutron stars (NSs) are primary targets for continuous gravitational-wave (GW) searches with the current network of ground-based detectors, LIGO, Virgo, and KAGRA. In these stars a deformation, or "mountain," sustained by elastic or magnetic strains, may result in a time-varying quadrupole from rotation, leading to the emission of GWs. Similarly, modes of oscillation may also lead to GW emission (see Lasky 2015 for a review).
In particular, NSs in low-mass X-ray binaries (LMXBs) are some of the most promising sources. In these systems magnetically channeled accretion from the companion onto the NS provides a mechanism to create a "mountain" (Ushomirsky et al. 2000; Melatos & Payne 2005; Osborne & Jones 2020; Singh et al. 2020), and the resulting GW torque may provide the solution to an astrophysical conundrum. There appears to be a sharp observed cutoff in the spin frequency (νs ) distribution of NSs in LMXBs at νs ≈ 750 Hz (Chakrabarty et al. 2003; Patruno et al. 2017), well below the theoretical breakup frequency for an NS (Haskell et al. 2018).
Although there are still several uncertainties in the modeling of the spin-up accretion torques (Patruno & Watts 2021; Glampedakis & Suvorov 2021), which may explain this observation (Patruno et al. 2012; Ertan & Alpar 2021), it has been suggested that the spin-down GW torques due to mountains can lead to an equilibrium at high frequencies that naturally explains the observed spins of NSs in LMXBs (Papaloizou & Pringle 1978; Wagoner 1984; Bildsten 1998) and the clustering of systems above νs ≈ 500 Hz and close to the maximum frequency (Patruno et al. 2017; Gittins & Andersson 2019). In such a scenario there is a natural correlation between the observed X-ray flux and the expected strength of the GWs, as a higher accretion rate leads to a stronger spin-up torque and thus requires a stronger GW torque for equilibrium (note that even if this equilibrium holds on average, there is expected to be some slight fluctuation in frequency or "spin wandering"; Bildsten et al. 1997; Mukherjee et al. 2018). Scorpius X-1 (Sco X-1), the most luminous LMXB, which is presumed to consist of an NS of mass ≈1.4M⊙ in a binary orbit with a companion star of mass ≈0.4M⊙ (Steeghs & Casares 2002), is therefore a very promising potential source of GWs. Some of the parameters inferred from electromagnetic observations of the system are summarized in Table 1. Note that the orbital eccentricity of Sco X-1 is believed to be small (Steeghs & Casares 2002; Wang et al. 2018) and is ignored in this search. Inclusion of eccentric orbits would add two search parameters that are determined by the eccentricity and the argument of periapse (Messenger 2011; Leaci & Prix 2015).
Table 1. Observed Parameters of the LMXB Sco X-1
| Parameter | Value |
|---|---|
| Right ascension a | 16h19m55.0850s |
| Decl. a |
|
| Distance (kpc) | 2.8 ± 0.3 |
| Orbital inclination i b | 44(deg) ± 6(deg) |
| K1 (km s−1) c | [40, 90] |
| tasc (GPS s) d | 974,416,624 ± 50 |
| Porb (s) e | 68,023.86 ± 0.043 |
Notes. Uncertainties are 1σ unless otherwise stated. There are uncertainties (relevant to the present search) in the projected velocity amplitude K1 of the NS, the orbital period Porb, and the time tasc at which the NS crosses the ascending node (moving away from the observer), measured in the solar system barycenter.
a The sky position (as quoted in Abbott et al. 2007a, derived from Bradshaw et al. 1999) is determined to the microarcsecond and therefore can be treated as known in the present search. b The inclination i of the orbit to the line of sight, from observation of radio jets in Fomalont et al. (2001), is not necessarily the same as the inclination angle ι of the NS's spin axis, which determines the degree of polarization of the GW in Equation (1). c The projected orbital velocity K1 as estimated by Doppler tomography measurements and Monte Carlo simulations in Wang et al. (2018), which show K1 to be weakly determined beyond the constraint that 40 km s−1 ≲ K1 ≲ 90 km s−1. d The time of ascension tasc, at which the NS crosses the ascending node (moving away from the observer), measured in the solar system barycenter, is derived from the time of inferior conjunction of the companion given in Wang et al. (2018) by subtracting Porb/4. It corresponds to a time of 2010 November 21 23:16:49 UTC and can be propagated to other epochs by adding an integer multiple of Porb, which results in increased uncertainty in tasc and correlations between Porb and tasc; see Figure 2. e The orbital period reported in Wang et al. (2018). Note that this differs from the previous estimate in Galloway et al. (2014) by 2.6σ.References. Bradshaw et al. (1999); Fomalont et al. (2001); Wang et al. (2018).
Download table as: ASCIITypeset image
Given its promise as a source for potentially detectable continuous gravitational waves, Sco X-1 has been the subject of numerous GW searches and search methods to date, starting with a fully coherent search (Jaranowski et al. 1998) of 6 hr from initial LIGO's second science run (Abbott et al. 2007a). Beginning with the fourth science run, results for Sco X-1 have been reported (Abbott et al. 2007b; Abadie et al. 2011) as part of a search for stochastic signals from isolated sky positions, also known as the radiometer search (Ballmer 2006), which has continued in Advanced LIGO's first three observing runs (Abbott et al. 2017a, 2019, 2021a). Sco X-1 has also been included in a search principally designed for unknown binary systems (Goetz & Riles 2011) and subsequently improved to search directly for Sco X-1 (Meadors et al. 2016); these searches were applied to data from initial LIGO's fifth and sixth science runs (Aasi et al. 2014; Meadors et al. 2017). A search method connected to Doppler-modulated sidebands (Messenger & Woan 2007; Sammut et al. 2014) was developed and applied to data from initial LIGO's sixth science run (Aasi et al. 2015a). This was further adopted into the so-called Viterbi search (Suvorova et al. 2016, 2017), which uses a hidden Markov model to track possible spin wandering; the Viterbi search has been applied to data from Advanced LIGO's first three observing runs (Abbott et al. 2017b, 2022a). The cross-correlation method (Dhurandhar et al. 2008; Whelan et al. 2015) used in the present work is an extension of the radiometer search that uses the signal model of GWs from an LMXB such as Sco X-1 to look for correlations between data at different times. It has been applied to data from Advanced LIGO's first and second science runs to set the strongest limits so far on GWs from Sco X-1 (Abbott et al. 2017c; Zhang et al. 2021).
The Advanced LIGO Gravitational-Wave Observatory (Aasi et al. 2015b) has conducted three observing runs, the last two in coordination with Advanced Virgo (Acernese et al. 2015). In these three runs, transient GWs were detected from over 90 coalescences of binary systems of black holes and/or NSs (Abbott et al. 2021b). The LIGO-Virgo O3 observing run (Buikema 2020; Abbott et al. 2020) began on 2019 April 01 15:00:00 UTC (GPS 1238166018), continued until a commissioning break at 2019 October 01 15:00:00 UTC (GPS 1253977218), resumed on 2019 November 01 15:00:00 UTC (GPS 1256655618), and ended on 2020 March 27 17:00:00 UTC (GPS 1269363618). In 2020 April, immediately following the LIGO-Virgo run, the KAGRA detector (Aso et al. 2013; Akutsu et al. 2021) and the GEO 600 detector (Lück et al. 2010; Affeldt et al. 2014; Dooley et al. 2016) conducted joint observations (Abbott et al. 2022b). In this analysis, we use data from the two LIGO detectors, as Virgo and KAGRA data were significantly less sensitive.
We use the calibrated data that are limited to times when a detector was in scientific observing mode (Davis et al. 2021). Due to the presence of transient instrumental glitches that degrade the sensitivity by raising the overall noise spectrum, we apply the "self-gating" procedure (Zweizig & Riles 2021) to remove these glitches when analyzing data below 600 Hz. This reduces the total volume of data included in the analysis below 600 Hz from 243–244 days to 231–240 days for the LIGO Hanford detector and from 250–251 days to 216–248 days for the LIGO Livingston detector. (The ranges are due to differences in the time baseline used in producing Fourier transforms at different frequencies; see Section 3.)
In addition, as in Abbott et al. (2017c), we exclude from our analysis frequencies at which the data are known to be influenced by instrumental disturbances of narrow frequency extent, known as "lines." In practice, this procedure removes data from times at which the signal model has Doppler-shifted the GW signal frequency f0 of the search template into the instrumental line, reducing the sensitivity of the search near known lines.
The remainder of the paper is laid out as follows: In Section 2 we describe the properties of Sco X-1 and the modeled GW signal from it. In Section 3 we describe the specifics of the cross-correlation search as implemented for this analysis. In Section 4 we describe the identification and follow-up of potential signals. Section 5 sets upper limits on the strength of GWs from Sco X-1 from the sensitivity and result of the search on Advanced LIGO data and simulated signals and considers their implications on various torque balance models. Finally, Section 6 contains the conclusions.
2. Model of Gravitational Waves from Sco X-1
The modeled GW signal from a rotating NS consists of a "plus" polarization component ![${h}_{+}(t)={A}_{+}\cos [{\rm{\Phi }}(t)]$](https://content.cld.iop.org/journals/2041-8205/941/2/L30/revision1/apjlaca1b0ieqn2.gif)
and a "cross" polarization component ![${h}_{\times }(t)={A}_{\times }\sin [{\rm{\Phi }}(t)]$](https://content.cld.iop.org/journals/2041-8205/941/2/L30/revision1/apjlaca1b0ieqn3.gif)
.
304
The signal recorded in a particular detector will be a linear combination of h+ and h× determined by the detector's orientation as a function of time. The two polarization amplitudes are
where h0 is an intrinsic amplitude describing the strength of the signal when it reaches the solar system and ι is the inclination of the NS's spin to the line of sight. (For an NS in a binary, the spin inclination ι is not necessarily equal to the inclination i of the binary orbit.) If ι = 0° or 180°, A× = ± A+, and gravitational radiation is circularly polarized. If ι = 90°, A× = 0, and it is linearly polarized. The general case, elliptical polarization, has 
. Many search methods are sensitive to the combination
which is equal to 
for circular polarization and 
for linear polarization (Messenger et al. 2015; this was the convention used in Abbott et al. 2017c but differs by a factor of 2.5 from the definition of 
in Whelan et al. 2015).
In order to understand the astrophysical relevance of the GW strengths we are probing, a useful benchmark is the so-called torque balance level. As already mentioned, it has been suggested that an LMXB, such as Sco X-1, may be in an equilibrium state where the spin-up torque due to accretion is balanced by a spin-down torque due to GW emission (Papaloizou & Pringle 1978; Wagoner 1984; Bildsten 1998). In order to obtain an estimate of the GW amplitude, we start by taking a simple spin-up torque of the form (Pringle & Rees 1972)
where 
is the mass accretion rate onto the NS, which we infer from the X-ray flux FX
; M is the mass of the star; G is the gravitational constant; and r is the lever arm, i.e., the radius at which the accretion torque is applied. By balancing the spin-up torque with the GW torque, i.e., imposing 
, we can obtain the GW amplitude (Watts et al. 2008):
Note that the usual torque balance benchmark assumes that accretion occurs at the surface of the NS, r = R* = 10 km. (If the magnetic field is strong enough to truncate the accretion disk above the surface, the lever arm will instead be the Alfvén radius, r = rA, and the GW amplitude implied by Equation (4) will be larger, as given in, e.g., Zhang et al. (2021) and Abbott et al. (2022a).) For Sco X-1, using the observed X-ray flux FX = 3.9 × 10−7 erg cm−2 s−1 from Watts et al. (2008) and assuming that the GW frequency f0 is twice the spin frequency νs (as would be the case for GWs generated by triaxiality in the NS), the torque balance value is
It is important to note that this amplitude is simply an order-of-magnitude estimate, which we use as a benchmark to understand whether our searches are probing astrophysically significant portions of parameter space. Much of the physics entering the accretion torque is, in fact, highly uncertain and depends on unknown physical parameters, such as the topology of the stellar magnetic field, the disk−field coupling, viscous heating in the disk, efficiency of X-ray emission, or radiation pressure in the disk. All these effects can strongly influence the spin-up torque, leading not only to a large rescaling (of up to an order of magnitude) of the strength of the torque in Equation (3) but in general also to different scalings with the parameters of the system (Patruno et al. 2012; Haskell et al. 2015; Glampedakis & Suvorov 2021). For example, Andersson et al. (2005) have even suggested that, for high accretion luminosities, radiation pressure will lead to a sub-Keplerian disk and a strongly reduced spin-up torque. In this case Sco X-1 would host a slowly rotating NS, which does not emit GWs in our current search band. In light of the various uncertainties, we will retain the standard simplifying assumptions in the derivation of Equation (5) for most of our torque balance comparisons, but keep in mind that the torque balance level is uncertain and model dependent and should not be interpreted too strictly.
3. Setup of Cross-correlation Search
The cross-correlation (CrossCorr) search method (Dhurandhar et al. 2008; Whelan et al. 2015) has been used to search for GWs from Sco X-1 in LIGO data from the first two observing runs of Advanced LIGO and Advanced Virgo (Abbott et al. 2017c; Zhang et al. 2021). It uses the signal model described in Section 2 to construct an appropriately weighted statistic ρ including correlations between data separated by up to a coherence time 
. The statistic is constructed using short Fourier transforms (SFTs) of length Tsft. If the SFTs are labeled by an index K, L, etc., which encodes the detector and time of the SFT, and zK
is an appropriately normalized combination of the Fourier data at the frequency of interest, we can write the statistic ρ as
where 
is the set of all pairs of SFTs whose start times differ by 
or less and WKL
is a complex weighting factor constructed using the signal model. Since the choice of frequency bin(s) in the construction of zK
and the amplitude and phase of the weighting factor WKL
for each SFT pair depend on the unknown parameters of the signal, we must conduct the search at a set of points in parameter space, each of which defines a "search template."
The maximum separation 
can be chosen to "tune" the search: higher 
values produce a more sensitive search but can significantly increase computing cost, both due to the increased number of correlation terms in the statistic and especially due to the increased density of search templates needed in the parameter space. As detailed in Abbott et al. (2017c), the Tsft and 
values were chosen as a function of GW frequency and orbital parameters in order to optimize the search at a given computing cost. For the present search, we used the same 
values as in O1 rather than re-optimizing.
305
The one exception is for the GW frequency range 400–600 Hz, for which the achievable sensitivity is closer in O3 than it was in O1 to the signal strength nominally expected from the torque balance model Equation (5); for those frequencies, we used double the 
of the O1 search. Note that, even with the same coherence times, the O3 search would require more computing resources than the O1 search, due to the increased observing time. However, by using a more efficient template lattice and convenient coordinate choices as described in Wagner et al. (2022), we are able to offset the increase in observing time and maintain a manageable computing time.
In addition to the GW signal frequency, f0, we search over the orbital parameters of the system, as summarized in Table 2. The projected semimajor axis of the orbit is assumed to lie in the range ![$a\sin i\in [1.44,3.25]$](https://content.cld.iop.org/journals/2041-8205/941/2/L30/revision1/apjlaca1b0ieqn18.gif)
lt-s, corresponding to a range in projected orbital velocity of [40, 90] km s−1. The search region in orbital period Porb and time of ascension tasc is constructed using the method of Wagner et al. (2022): the time of ascension tasc is propagated 4125 orbits to define
and the "sheared" orbital period is defined as
The most likely 
is 2019 October 13 15:17:11 UTC (GPS 1255015049). The coordinates 
and 
are approximately uncorrelated both in the parameter space metric of the search and in the astrophysical prior distribution.
306
Note that the included prior probability is an underestimate of the efficiency in covering parameter space, since the "sheared" period coordinate 
was unresolved in many search jobs, i.e., the mismatch associated with an offset of 3.3 × 0.011 s from the most likely value 68023.86 s (the dashed rectangular boundaries in Figure 1) was less than 0.0625.
Figure 1. Search region in terms of parameters 
and 
defined in Equations (7) and (8), respectively. The lattice is constructed to completely cover, with maximum mismatch 0.25, the solid black truncated ellipse. The solid black ellipse has semiaxes 3.3 × 185 s in 
and 3.3 × 0.011 s in 
, and the truncating boundaries are at ±3 × 185 s of the most likely 
value. For comparison, the thin colored ellipses show curves of constant prior probability corresponding to 1σ, 2σ, and 3σ (containing 39.3%, 86.5%, and 98.9% of the prior probability, respectively), including the effects of changing coordinates from tasc and Porb appearing in Table 1 to 
and 
. The inner search region, in which we choose a longer 
to do a deeper search, is within ±185 s of the most likely value of 
and contains 68.1% of the prior probability, while the full search region, within ±3 × 185 s of the most likely 
value, contains 99.2% of the prior probability. The slight misalignment of the prior and search ellipses is due to a software bug, which led to a definition of 
that differed slightly from the optimal one, as described in Section 3. The dashed rectangular boundaries show the region effectively covered by the majority of search jobs for which the "sheared" period coordinate 
was unresolved in many search jobs, i.e., the mismatch associated with an offset of 3.3 × 0.011 s from the most likely value 68023.86 s was less than 0.0625.
Download figure:
Standard image High-resolution imageTable 2. Parameters Used for the Cross-correlation Search
| Parameter | Range |
|---|---|
| f0 (Hz) | [25, 1600] |
 (lt-s)
a
| [1.44, 3.25] |
 (GPS s)
b
| 1,255,015,049 ± 3 × 185 |
 (s)
c
| 68,023.86 ± 3.3 × 0.011 |
Notes.
a The range for the projected semimajor axis
in lt-s was taken from the constraint K1 ∈ [40, 90] km s−1.
b
This value for the time of ascension 
, defined in Equation (7), has been propagated forward by 4125 orbits from the value of tasc in Table 1 and corresponds to a time of 2019 October 13 15:17:11 UTC, near the middle of the O3 run. The increase in uncertainty is due to the uncertainty in Porb.
c
This is the "sheared" period defined in Equation (8); note that the uncertainty in 
has been reduced compared to the marginal uncertainty in Porb by the same factor by which the uncertainty in 
has been increased relative to that for tasc, as described in Wagner et al. (2022). The search region in 
is given by the elliptical boundary shown in Figure 1, but it is defined as "unresolved" if only one template is needed to cover 68,023.86 ± 3.3 × 0.011 at a maximum mismatch of 0.0625.Download table as: ASCIITypeset image
The search was done over a range of 
and with 
constrained to lie in an elliptical region centered on 68,023.86 s with semiaxes of 3.3 × 0.011 s for 
and 3.3 × 185 s for 
, as illustrated in Figure 1. For reference, in Figure 2 we show this region in the coordinates 
, along with the search regions used in Abbott et al. (2017c) and Zhang et al. (2021), propagated forward in time to the epoch considered in this analysis.
Figure 2. The search regions and prior uncertainties shown in Figure 1, expressed in terms of the system parameters 
and Porb. For reference, we also show the regions used for the O1 analysis in Abbott et al. (2017c; gray dashed lines) and reported for the O2 analysis in Zhang et al. (2021; orange dotted lines), propagated to the epoch of this search (which transforms the rectangular search regions into parallelograms). Note that the search region for Abbott et al. (2017c) is offset in both Porb and 
because it used the Porb estimate in Galloway et al. (2014), while the others used the updated estimate in Wang et al. (2018). The analysis of Abbott et al. (2017c) is still believed to have covered the plausible signal space because of the underresolution of the period direction and the fact that the offset in 
induced by the inaccurate Porb value was less for the epoch of the search (2015 rather than 2019).
Download figure:
Standard image High-resolution imageTo perform the analysis, the parameter space was divided into small jobs that could be run in parallel. Each 5 Hz GW frequency band was divided into between 100 and 4000 sub-bands, and the orbital parameter space was divided into 9 or 20 cells, as illustrated in Figure 3. These subdivisions of parameter space were chosen so that each analysis job would run in a reasonable amount of time (≲10 hr), allowing the analysis to be done quickly via distributed computing.
Figure 3. Illustration of parameter space cells in 
and 
and example of coherence times 
, in seconds, chosen as a function of the orbital parameters of the NS. Increasing coherence time improves the sensitivity but increases the computational cost of the search. The values used are the same as in Abbott et al. (2017c), except between 400 and 600 Hz, where they have been doubled; see Table 3 for details.
Download figure:
Standard image High-resolution imageSearch templates were placed in the parameter space using an 
lattice with a maximum mismatch of 0.25 as described in Wagner et al. (2022) and Wette (2014). Because the use of the sheared coordinate 
reduces the 1σ prior uncertainty from 0.043 to 0.011 s, the period was unresolved for most search jobs. We computed the maximum mismatch 
associated with an offset of 3.3 × 0.011 s. For jobs in which 
, the initial search was done with 
and an 
lattice in the other three parameters 
. To account for the nonzero mismatch contribution from the possible 
offset, the maximum mismatch of the 
lattice was set to 
. (Note that this actually guarantees coverage at the specified maximum mismatch over a rectangle in 
rather than just the ellipse shown in Figure 1.) If, on the other hand, the 
coordinate was resolved in a particular search job, an 
lattice in all of the coordinates with maximum mismatch of 0.25 was used.
As noted in Section 1, even in approximate equilibrium, Sco X-1 may undergo stochastic variation of the GW frequency f0, also known as "spin wandering." As in Abbott et al. (2017c), we can apply the estimates in Whelan et al. (2015) to set bounds on the loss of signal-to-noise ratio (S/N) under a simplistic model in which the GW frequency undergoes a net spin-up or spin-down of magnitude 
, changing on a timescale Tdrift. For O3, in which the duration of the run from start to end is Trun = 3.12 × 107 and the coherence times used for the initial search are 
, the expected loss of S/N is
which indicates that spin wandering is not likely to be an important effect for this search. In addition, as noted in Zhang et al. (2021), the predictions of Mukherjee et al. (2018) based on the time variation of the X-ray flux from Sco X-1 imply considerably less spin wandering than the naive 
−Tdrift model, whose fiducial parameters are taken from Messenger et al. (2015).
4. Candidates, Outliers, and Follow-up
The detection statistic ρ is normalized to have zero mean and unit variance in Gaussian noise. As in Abbott et al. (2017c), we make a naive estimate of the expected background by assuming that each search template represents an independent Gaussian random number, and we use this value to set the threshold at an approximate level of one expected false-alarm level per 50 Hz. As shown in Figure 4, the threshold for follow-up for this search was set at 6.3 for 25 Hz < f0 < 400 Hz, 6.2 for 400 Hz < f0 < 600 Hz, and 5.8 for 600 Hz < f0 < 1600 Hz. Due in part to our more efficient template placement, we were able to use a lower follow-up threshold than in Abbott et al. (2017c), except for 400 Hz < f0 < 600 Hz, where we used the same threshold, despite having a search with twice the coherence time. Table 3 shows the resulting expected numbers of false alarms in each GW frequency range.
Figure 4. Selection of follow-up threshold as a function of GW frequency. If the data contained no signal and only Gaussian noise, each template in the parameter space would have some chance of producing a statistic value exceeding a given threshold. Within each 5 Hz frequency band, the total number of templates was computed and used to find the threshold at which the expected number of Gaussian outliers above that value would be 0.1. This is shown with short blue lines for the templates in the present search (see Table 3); for reference, the thresholds calculated from the numbers of templates in the O1 search of Abbott et al. (2017c) are shown in orange. Because of the more efficient template placement of algorithm from Wagner et al. (2022), the O3 search has fewer templates, and therefore a lower implied threshold, than the O1 search, which used the same coherence times. The exception is for 400 Hz < f0 < 600 Hz, where the same threshold of 6.2 was used for the O1 and O3 searches, and the latter used twice the coherence time (and therefore a denser parameter space lattice) as the former. The present search uses a threshold of 6.3 for 25 Hz < f0 < 400 Hz and 5.8 for 600 Hz < f0 < 1600 Hz (black dashed line), which are lower than in the O1 search (magenta dashed line). Note that the large number of non-Gaussian outliers (see Table 3) makes the Gaussian follow-up level an imprecise tool in any event.
Download figure:
Standard image High-resolution imageTable 3. Summary of Numbers of Templates and Candidates
| f0 (Hz) | Tsft |
 (s) | ρ | Number of | Expected Gauss | Follow-up Level | |||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| Min | Max | (s) | Min | Max | Thresh a | Templates | False Alarms b | 0 c | 1 d | 2 e | 3 f |
| 25 | 50 | 1440 | 10080 | 18720 | 6.3 | 5.68 × 109 | 0.8 | 63 | 31 | 10 | 1 |
| 50 | 100 | 1020 | 8160 | 14280 | 6.3 | 2.88 × 1010 | 4.3 | 131 | 114 | 40 | 2 |
| 100 | 150 | 840 | 6720 | 10920 | 6.3 | 6.13 × 1010 | 9.1 | 169 | 166 | 73 | 3 |
| 150 | 200 | 720 | 5040 | 8640 | 6.3 | 6.69 × 1010 | 10.0 | 171 | 170 | 68 | 2 |
| 200 | 300 | 600 | 2400 | 4800 | 6.3 | 4.54 × 1010 | 6.7 | 66 | 66 | 14 | 4 |
| 300 | 400 | 510 | 1530 | 3060 | 6.3 | 2.09 × 1010 | 3.1 | 19 | 19 | 1 | 1 |
| 400 | 600 | 360 | 720 | 2160 | 6.2 | 3.46 × 1010 | 9.8 | 343 | 226 | 20 | 9 |
| 600 | 800 | 360 | 360 | 360 | 5.8 | 1.80 × 109 | 6.0 | 15 | 15 | 0 | 0 |
| 800 | 1200 | 300 | 300 | 300 | 5.8 | 4.36 × 109 | 14.5 | 226 | 70 | 2 | 0 |
| 1200 | 1600 | 240 | 240 | 240 | 5.8 | 4.36 × 109 | 14.5 | 346 | 55 | 2 | 0 |
Notes. For each range of GW frequencies, this table shows the SFT duration Tsft; the minimum and maximum coherence time 
used for the search, across the different orbital parameter space cells (see Figure 3); the threshold in S/N ρ used for follow-up; the total number of templates; and the number of candidates at various stages of the process. (See Section 4 for detailed description of the follow-up procedure.)
equal to 4× the original 
for that candidate. (True signals should approximately double their S/N; any candidates whose S/N goes down have been dropped.) All of the signals present at this stage are shown in Figure 5, which also shows the behavior of the search on simulated signals injected in software.
f
This is the number of candidates remaining after 
has been increased to 16× its original value.Download table as: ASCIITypeset image
For candidates exceeding the follow-up threshold, we applied the procedure detailed in Abbott et al. (2017c). Here we highlight the basic steps, as well as details that were changed for this analysis:
- 1.Candidates were "clustered" together in GW frequency, with all templates within 0.01 Hz of a peak in S/N above the threshold being represented by the parameters of the peak. These are known as the "level 0" results.
- 2.A "refinement" search was performed on each level 0 candidate, with the same
as the original search and a resolution ∼3× as fine as the original lattice. This and later stages of follow-up were run on a rectangular grid in the "sheared" parameters . For the refinement stage, a grid of 13 × 13 × 13 × 5 points was used, centered on the values of the candidate, and covering the full prior range in the initially unresolved . To deal with the effects of unknown narrowband features ("lines") present in a single detector, we computed a detection statistic using only data from the LIGO Livingston Observatory (LLO) detector and another using only LIGO Hanford Observatory (LHO) data. If either of these exceeded the detection statistic constructed from all the data, we vetoed the candidate as a likely instrumental artifact. Candidates from the refinement search that survive this veto are known as the "level 1" results. - 3.Two successive rounds of follow-up were performed, starting with the level 1 candidates. At each stage, the coherence time was quadrupled from the previous stage, and the density of templates in each direction was increased by a factor of 3. The grid used was 13 × 13 × 13 × 13 in , centered on the peak of the previous level's results. Candidates that increased their S/N from the previous level were known as the "level 2" (4× the original ) and "level 3" (16× the original ) results.
A total of 22 candidates survive level 3 of follow-up. Two checks were done to determine whether any of them represent convincing detection candidates: one using the results of the cross-correlation search, and one using an independent pipeline.
For the first check, in Figure 5 we plot the ratio by which the S/N increases from level 1 to level 2 and from level 2 to level 3. We also plot the corresponding ratios for all of the candidates surviving level 2 (the 16× original 
follow-up is not available for candidates that fail level 2), as well as for the simulated signal injections described in Section 5. We would naively expect a real signal to double its S/N between level 1 and level 2, and again between level 2 and level 3, but none of the candidates from the search come close to this. As in Abbott et al. (2017c), none of the candidates double their S/N from level 1 to level 3, let alone in a single follow-up stage. On the other hand, all but one of the injections, while not doubling their S/N with each stage of follow-up, increase their S/N noticeably more than any of the candidates from the search.
Figure 5. Ratios of follow-up statistics for search candidates and simulated signals. This plot shows all of the candidates that survived level 2 of follow-up (see Section 4 and Table 3), both from the main search and from the analysis of the simulated signal injections described in Section 5. It shows the ratios of the S/N ρ after follow-up level 1 (at the original coherence time 
), level 2 (at 4× the original coherence time), and level 3 (at 16× the original coherence time). (The boxes labeled "single detector" are outliers or injections at GW frequencies where only one detector's data were included in the analysis because of known instrumental artifacts in the other detector.) The green dashed lines are at constant values of ρlevel 3/ρlevel 1 equal to 2 and 4. There are no points with ρlevel 2/ρlevel 1 < 1 because those candidates do not survive level 2 follow-up and are therefore not subjected to level 3 follow-up. From the construction of the statistic in Whelan et al. (2015), the naive expectation is that the S/N will roughly double each time 
is quadrupled. Empirically, the follow-ups of injections do not show exactly that relationship, but all but one (which was injected at a frequency contaminated by instrumental artifacts) show significant increases in S/N that are not seen in any of the follow-ups of search candidates. We thus conclude that no convincing detection candidates are present.
Download figure:
Standard image High-resolution imageNote that, at some GW signal frequencies, known instrumental lines lead us to omit all the data from one of the two LIGO detectors from the search. This produces a single-detector search for which the unknown-line veto is not applicable. Of the 22 search outliers that survived our veto process, 6 were in this category, as well as 4 of the injections, including the one injection that increased its S/N negligibly under the follow-up procedure. Also note that of the 821 injected signals (out of 918) that produced ρ values above their respective thresholds, 817 survived all the levels of follow-up. (There were three vetoed at level 1, one at level 2, and zero at level 3, all because of the single-detector unknown-line veto.) We thus conclude that our follow-up procedure is relatively robust and that there are no convincing detection candidates from the search.
An additional, complementary, multistage MCMC follow-up using the method described in Tenorio et al. (2021) using the PyFstat package (Ashton & Prix 2018; Keitel et al. 2021; Ashton et al. 2022) was applied to the 22 outliers using the same configuration as in Abbott et al. (2022c). This method places templates adaptively to compute the semicoherent 
-statistic (Jaranowski et al. 1998; Cutler & Schutz 2005) around the candidate of interest using a diminishing number of coherent segments (660, 330, 92, 24, 4, and 1). The coherence times of the corresponding segments range from half a day to the full observing run. A Bayes factor is computed using the 
-statistic values from subsequent coherence stages corresponding to the loudest template. The signal hypothesis assesses the consistency of these values, whereas the noise hypothesis states the (in)consistency of the final value with the background distribution, taking the final-stage factor into account.
The resulting Bayes factor values are significantly lower than what would be expected for a signal detectable by this search, as confirmed by an analogous follow-up of a similar number of injected signals.
5. Upper Limits and Implications
Since our search produced no convincing detection candidates, we set upper limits on the strength of GWs from Sco X-1 as a function of GW frequency, using the method described in detail in Abbott et al. (2017c). First, naive Bayesian upper limits were set within each 0.05 Hz band, using the 95th percentile of the posterior on 
or 
deduced from the highest S/N seen in each band.
307
Then, a series of simulated signals were added to the data at a variety of amplitudes, and a Bayesian logistic regression analysis was performed to estimate the factor by which to multiply the naive 95% Bayesian upper limit on amplitude in order to reach the threshold of 95% signal recovery. (As in Abbott et al. (2017c), "recovery" was defined as an increase in the maximum S/N seen in a band, over the value with no injection present.) We performed a total of 918 injections
308
between 25 and 500 Hz, of which 863 were recovered, with a resulting adjustment factor of 1.19 to the amplitude of the 
upper limit. As described in Abbott et al. (2017c), to set the adjustment factor for the h0 limit, we limit attention to injections that were generated with a specified h0 rather than 
, of which we recovered 546 of 575, with a resulting adjustment factor of 1.17.
309
The upper limits including these adjustment factors are shown in Figure 6.
Figure 6. Upper limits from directed searches in advanced LIGO data. Top: upper limit on h0, after marginalizing over NS spin inclination ι, assuming an isotropic prior. The dashed line shows the nominal expected level assuming torque balance (Equation (5)) as a function of GW frequency. Bottom: upper limit on 
, defined in Equation (2). This is equivalent to the upper limit on h0 assuming circular polarization. (Note that the marginalized upper limit in the top panel is dominated by linear polarization and so is a factor of almost 
higher). The blue dotted–dashed line (labeled as "tb w/ι = 44°") corresponds to the assumption that the NS spin is aligned to the most likely orbital angular momentum and ι ≈ i ≈ 44° (see Table 1). The blue diagonal bands show 
levels corresponding to the torque balance h0 in the top panel. The darker-shaded band corresponds (5th to 95th percentiles) to a Gaussian distribution with mean and standard deviation corresponding to ι = 44° ± 6°, as used in Zhang et al. (2021). Finally, the lighter-shaded band shows the full range of possible 
values corresponding to torque balance, with circular polarization at the top and linear polarization on the bottom. For comparison with the "CrossCorr O3" results presented in this paper, we show in the top panel the isotropic marginalized limits from the previous cross-correlation searches in Abbott et al. (2017c) ("CrossCorr O1") and Zhang et al. (2021) ("CrossCorr O2"). In the bottom panel we include the limits assuming circular polarization from other searches of O3 data: "Radiometer O3" is the narrowband radiometer analysis of Abbott et al. (2021a), which used data from Advanced LIGO's first three observing runs, and "Viterbi O3" is the analysis of Abbott et al. (2022a) using a hidden Markov model.
Download figure:
Standard image High-resolution imageThe upper limits placed by this search improve on those from previous observing runs and are now probing a theoretically significant portion of parameter space. This is usually quantified in terms of the torque balance amplitude, i.e., the GW amplitude that would be required for GW spin-down torques to balance the accretion torque. As discussed in Section 1, this limit is a useful benchmark but is highly uncertain, reflecting the high level of uncertainty in the theoretical modeling of accretion torques. We illustrate this in Figure 7 by comparing the upper limits from our search not only to the standard limit of Equation (5), obtained by assuming the lever arm to be the radius of the NS, r = R* = 10 km, in Equation (4), but also to an example of the range of torque balance amplitudes allowed by current models for accretion onto a magnetized star. First, we consider a model of the same form as in Equations (3) and (5), but we do not assume accretion to occur on the surface, but rather take the torque arm to be the Alfvén radius rA, at which the disk is truncated by the magnetic field (Pringle & Rees 1972):
where 0.1 ≲ X ≲ 1 is a phenomenological parameter that encodes the uncertainty in the truncation radius of the disk. Using the mass accretion rate of Sco X-1 inferred from X-ray observations (Watts et al. 2008), along with B = 108 G and X = 1, gives an Alfvén radius of rA ≈ 49 km, which is used to generate the curve in Figure 7. We also consider one of the parameterized models of Glampedakis & Suvorov (2021), which encompass a wide range of physics and can successfully fit spin-up episodes in a number of observed LMXBs. In particular, for illustrative purposes, we use their "new" model 2, for which the torque has the form 
, with the fastness parameter 
and the corotation radius 
. We consider two values for the magnetic field strength at the surface of the star, B = 108 and 109 G, and consider values of ξ between 0.3 (which sets the upper limit in our plots) and 0.5.
Figure 7. Comparison of upper limits to predictions of torque balance models. The gray band indicates the h0 upper limit implied by the 
upper limit in the bottom panel of Figure 6, assuming the range of possible inclinations from 
. (Linear polarization (
) is at the top, and circular polarization (
) is at the bottom.) The dashed blue line is the usual conservative torque balance estimate assuming accretion at the surface of the NS, r = R* = 10 km. The dotted–dashed red line is the same model assuming a lever arm of r = rA ≈ 49 km, which is the value of the Alfvén radius given by Equation (10) for ξ = 1 and B = 108 G. Note that this is slightly larger than the value of 35 km used in, e.g., Zhang et al. (2021) and Abbott et al. (2022a), since we use the inferred mass accretion rate of Sco X-1. The colored bands show the range of predictions for the models of Glampedakis & Suvorov (2021), in particular model 2, for the range of 0.3 < ξ < 0.5, assuming a magnetic field of 108 G (GS8) or 109 G (GS9).
Download figure:
Standard image High-resolution imageAs can be clearly seen, there is a wide portion of parameter space allowed by theoretical models. Furthermore, at higher frequencies the theory is intrinsically uncertain; it is in fact unclear whether one should expect the source to be rapidly rotating, as some of the models of Glampedakis & Suvorov (2021) predict spin equilibrium due to accretion alone below GW frequencies of roughly 1 kHz, without any need for gravitational radiation–in fact, it is also possible that the NS may not be in the frequency range we are searching. This "fuzziness" in the torque balance limit (which may be even further enhanced by additional effects not considered in the models of Glampedakis & Suvorov (2021), such as viscous heating in the disk or the efficiency of X-ray emission) makes it therefore impossible to draw firm conclusions on the equation of state (EOS) of the NS, or magnetic field strength, directly from our upper limits, without committing to a particular accretion model. It is therefore important to remember, when analyzing our results, that torque balance may not be active in this system, or that even if it is, the GW torques may only be active for a fraction of the time, with a low duty cycle for GW emission (Haskell et al. 2015).
Nevertheless, it is clear from Figure 7 that in the range we are most sensitive to, approximately between 30 and 400 Hz, our upper limits are probing below even the more stringent limit set by accretion onto the surface, thus searching a physically significant portion of parameter space. Furthermore, our results are probing the parameter space predicted by the models of Glampedakis & Suvorov (2021) up to f0 ∼ 500 Hz. Searching at these higher GW frequencies is important; although the frequency of Sco X-1 is not known, the distribution of spin frequencies of the observed AMXPs appears to be bimodal (Patruno et al. 2017), with a "fast" population of pulsars, for which Patruno et al. (2017) make the hypothesis that GW emission may play a role, centered around νs ≈ 550 Hz (i.e., f0 ≈ 1100 Hz for triaxial emission), and a slower population centered around νs ≈ 300 Hz (i.e., f0 ≈ 600 Hz).
To give an illustration of how our results can constrain possible torque balance models, consider in more detail the simplest accretion model, i.e., that obtained by setting r = R* in Equation (3). In this case we may ask whether, by comparing the upper limits from our searches to the theoretical value for the h0 in Equation (3), it is possible to put constraints on the physical parameters of the star for which the torque balance scenario is still viable. To answer this question, we consider two parameters: the mass of the star and the inclination angle ι, for two EOSs taken from the CompOSE database (Typel et al. 2015; Oertel et al. 2017; Typel et al. 2022) both for a softer EOS, GR15 (Gulminelli & Raduta 2015), and for a stiffer EOS, GPPVA (Grill et al. 2014). The results can be seen in Figure 8 for the mass of the NSs. We see that, for the range of GW frequencies in which our search is most sensitive, we can exclude the torque balance scenario for higher-mass NSs, especially in the case of a stiffer EOS. The GW amplitude is, however, clearly very strongly affected by the inclination angle, so in Figure 8 we also plot our constraints in terms of the inclination angle ι, holding the stellar mass fixed at M = 1.4M⊙. In this case we also see that we can rule out torque balance models with nearly circular polarization (small ι) over a wide range of frequencies for both choices of EOS.
Figure 8. Illustration of how the upper limits shown in Figure 6 can constrain models for Sco X-1 that include torque balance due to gravitational waves. In the top row, we show constraints on the NS mass, assuming torque balance due to GW at the specified frequency, in the simple model where r = R*. We consider two EOSs, a softer model, GR15 (Gulminelli & Raduta 2015), and a stiffer model, GPPVA (Grill et al. 2014). The largest exclusion region is for an NS inclination angle ι = 0° (or 180°), where the GWs would be circularly polarized. Assuming the worst-case scenario of linear polarization ι = 90° gives constraints that hold for any inclination, and the most likely value of ι = 44° (aligned with the binary orbit) gives an intermediate case. We see that for both EOSs (and especially for the stiffer EOS) the torque balance scenario can be excluded for higher-mass NSs for the GW frequency range in which our searches are more sensitive. Since the inclination angle plays a strong role in the constraints, we present in the bottom row, for the same EOSs but for a fixed mass of M = 1.4M⊙, the limits that can be set on ι. We see that, for both EOSs, our observations can exclude nearly circular polarized (small ι) GW emission at the torque balance level for a wide range of frequencies.
Download figure:
Standard image High-resolution image6. Conclusions and Outlook
We have presented the results of the most sensitive search to date for GWs from Sco X-1, using LIGO detector data from the third observing run of Advanced LIGO and Advanced Virgo. We have set upper limits across a range of GW signal frequencies 25 Hz < f0 < 1600 Hz, corresponding to NS spin frequencies of 12.5 Hz < νs < 800 Hz. The sensitivity of our search is now probing possible models of torque balance equilibrium over a range of GW frequencies spanning hundreds of hertz and, for the first time, approaches the standard conservative torque balance prediction even under pessimistic assumptions about NS inclination angle. We expect to see a further improvement in sensitivity from upcoming LIGO-Virgo-KAGRA observing runs (Abbott et al. 2020).
This material is based on work supported by NSF's LIGO Laboratory, which is a major facility fully funded by the National Science Foundation. The authors also gratefully acknowledge the support of 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 GEO 600 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 Netherlands Organization for Scientific Research (NWO) 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, the Department of Science and Technology, India, the Science & Engineering Research Board (SERB), India, the Ministry of Human Resource Development, India, the Spanish Agencia Estatal de Investigación (AEI), the Spanish Ministerio de Ciencia e Innovación and Ministerio de Universidades, the Conselleria de Fons Europeus, Universitat i Cultura and the Direcció General de Política Universitaria i Recerca del Govern de les Illes Balears, the Conselleria d'Innovació Universitats, Ciència i Societat Digital de la Generalitat Valenciana and the CERCA Programme Generalitat de Catalunya, Spain, the National Science Centre of Poland and the European Union—European Regional Development Fund, Foundation for Polish Science (FNP), the Swiss National Science Foundation (SNSF), the Russian Foundation for Basic Research, the Russian Science Foundation, the European Commission, the European Social Funds (ESF), the European Regional Development Funds (ERDF), the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the French Lyon Institute of Origins (LIO), the Belgian Fonds de la Recherche Scientifique (FRS-FNRS), Actions de Recherche Concertées (ARC) and Fonds Wetenschappelijk Onderzoek—Vlaanderen (FWO), Belgium, the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, the Natural Science and Engineering Research Council Canada, Canadian Foundation for Innovation (CFI), the Brazilian Ministry of Science, Technology, and Innovations, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Leverhulme Trust, the Research Corporation, the Ministry of Science and Technology (MOST), Taiwan, the United States Department of Energy, and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, INFN, and CNRS for provision of computational resources.
This work was supported by MEXT, JSPS Leading-edge Research Infrastructure Program, JSPS Grant-in-Aid for Specially Promoted Research 26000005, JSPS Grant-in-Aid for Scientific Research on Innovative Areas 2905: JP17H06358, JP17H06361, and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grant-in-Aid for Scientific Research (S) 17H06133 and 20H05639, JSPS Grant-in-Aid for Transformative Research Areas (A) 20A203: JP20H05854, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, National Research Foundation (NRF), Computing Infrastructure Project of KISTI-GSDC, Korea Astronomy and Space Science Institute (KASI), and Ministry of Science and ICT (MSIT) in Korea, Academia Sinica (AS), AS Grid Center (ASGC) and the Ministry of Science and Technology (MoST) in Taiwan under grants including AS-CDA-105-M06, Advanced Technology Center (ATC) of NAOJ, and Mechanical Engineering Center of KEK.
This paper has been assigned LIGO Document No. LIGO-P2100110-v13.
Software: LALSuite (LIGO Scientific Collaboration 2018), LatticeTiling (Wette 2014), PyFstat (Ashton & Prix 2018; Keitel et al. 2021; Ashton et al. 2022) numpy (Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020).
Footnotes
- 304
The preferred (spin-2) polarization basis is constructed from orthonormal unit vectors in the plane of the sky: one along the projection of the NS spin axis, and the other along the intersection of the NS equatorial plane with the plane of the sky. This basis is rotated relative to a fiducial north-and-east-on-the-sky basis by an (unknown) polarization angle ψ, equivalent to the position angle of the NS's polarization axis (Jaranowski et al. 1998; Prix & Whelan 2007).
- 305
- 306
The correlations would have been even smaller had the optimal coefficient 2.25 × 10−4 been used in Equation (8), but a software bug led to the use of the coefficient 2.42 × 10−4. The impact of this error is negligible, however, reducing the prior probability covered by the search region from 99.4% to 99.2%.
- 307
We used a simple extreme value likelihood assuming independent Gaussian distributions for the detection statistics from the templates in the initial bank. Future work may leverage more sophisticated methods of estimating this distribution, such as those of Tenorio et al. (2022).
- 308
These are the same injections that were used to validate the follow-up procedure, as described in Section 4.
- 309
For comparison, the adjustment factors in Abbott et al. (2017c) were 1.44 for h0 and 1.21 for
. The fact that the factors derived from the current search are comparable is evidence that the template bank modifications of Wagner et al. (2022) do not significantly reduce the sensitivity of the search.
