Model-based Cross-correlation Search for Gravitational Waves from the Low-mass X-Ray Binary Scorpius X-1 in LIGO O3 Data

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 h 0 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.


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 timevarying 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 e in a binary orbit with a companion star of mass ≈0.4M e (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).
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(Abbott et al. , 2019(Abbott et al. , 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 298 lsc-spokesperson@ligo.org 299Deceased, December 2021. 300Deceased, November 2022. 301virgo-spokesperson@ego-gw.it 302 Deceased, March 2022. 303kscboard-chair@icrr.u-tokyo.ac.jpOriginal 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.
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(Suvorova et al. , 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(Abbott et al. , 2022a)).The crosscorrelation 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).
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 f 0 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 followup 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.

Model of Gravitational Waves from Sco X-1
The modeled GW signal from a rotating NS consists of a "plus" . 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 h 0 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 which is equal to h 0 2 for circular polarization and h 8 0 2 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 h 0 eff 2

(
) 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  M is the mass accretion rate onto the NS, which we infer from the X-ray flux F X ; 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): 5.48 10 10 erg cm s 300 Hz 10 km 10 km 1.4 .
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 = r A , 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 F X = 3.9 × 10 −7 erg cm −2 s −1 from Watts et al. (2008) and assuming that the GW frequency f 0 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 orderof-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.

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 T max .The statistic is constructed using short Fourier transforms (SFTs) of length T sft .If the SFTs are labeled by an index K, L, etc., which encodes the detector and time of the SFT, and z K 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 T max or less and W KL is a complex weighting factor constructed using the signal model.Since the choice of frequency bin(s) in the construction of z K and the amplitude and phase of the weighting factor W KL 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 T max can be chosen to "tune" the search: higher T max 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 T sft and T max 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 T max values as in O1 rather than re-optimizing. 305he 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 T max 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, f 0 , we search over the orbital parameters of the system, as summarized in The most likely t asc ¢ is 2019 October 13 15:17:11 UTC (GPS 1255015049).The coordinates t asc ¢ and P ˜are approximately uncorrelated both in the parameter space metric of the search and in the astrophysical prior distribution. 306Note that the included prior probability is an underestimate of the efficiency in covering parameter space, since the "sheared" period coordinate P ˜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.86s (the dashed rectangular boundaries in Figure 1) was less than 0.0625.The search was done over a range of t 1255015049 asc ¢ =  3 185 s ´and with P ˜constrained to lie in an elliptical region centered on 68,023.86s with semiaxes of 3.3 × 0.011 s for P ãnd 3.3 × 185 s for t asc ¢ , as illustrated in Figure 1.For reference, in Figure 2 we show this region in the coordinates t P , asc orb ¢ ( ), 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.
To 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 subbands, 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.
Search templates were placed in the parameter space using an * A n 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 P ˜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 P ) .To account for the nonzero mismatch contribution from the possible P ˜offset, the maximum mismatch of the * A 3 lattice was set to 0.25 P max m -˜.(Note that this actually guarantees coverage at the specified maximum mismatch over a rectangle in t P , asc ¢ ˜rather than just the ellipse shown in Figure 1.) If, on the other hand, the P ˜coordinate was resolved in a particular search job, an * A 4 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 f 0 , 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  f drift | | , changing on a timescale T drift .For O3, in which the duration of the run from start to end is T run = 3.12 × 10 7 and the coherence times used for the initial search are  T 18720 s max , 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.

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 < f 0 < 400 Hz, 6.2 for

Notes.
a The range for the projected semimajor axis a i K P sin 2 1 orb p = ( ) in lt-s was taken from the constraint K 1 ä [40, 90] km s −1 .b This value for the time of ascension t asc ¢ , defined in Equation (7), has been propagated forward by 4125 orbits from the value of t asc 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 P orb .c This is the "sheared" period defined in Equation (8); note that the uncertainty in P ˜has been reduced compared to the marginal uncertainty in P orb by the same factor by which the uncertainty in t asc ¢ has been increased relative to that for t asc , as described in Wagner et al. (2022).The search region in P ˜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. 306The 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%.400 Hz < f 0 < 600 Hz, and 5.8 for 600 Hz < f 0 < 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 < f 0 < 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.
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 T max 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 f a i t P , sin , , 0 asc ¢ ( ˜).For the refinement stage, a grid of 13 × 13 × 13 × 5 points was used, centered on the f a i t , sin , 0 asc ¢ ( )values of the candidate, and covering the full prior range in the initially unresolved P ˜.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  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 t asc ¢ and 3.3 × 0.011 s in P ˜, and the truncating boundaries are at ±3 × 185 s of the most likely t asc ¢ 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 t asc and P orb appearing in Table 1 to t asc ¢ and P ˜.The inner search region, in which we choose a longer T max to do a deeper search, is within ±185 s of the most likely value of t asc ¢ and contains 68.1% of the prior probability, while the full search region, within ±3 × 185 s of the most likely t asc ¢ 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 P ˜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 P ˜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.86s was less than 0.0625.(2017c) is offset in both P orb and t asc ¢ because it used the P orb 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 t asc ¢ induced by the inaccurate P orb value was less for the epoch of the search (2015 rather than 2019).
Figure 3. Illustration of parameter space cells in t asc ¢ and a i sin and example of coherence times T max , 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.
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 T max 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 f a i t P , sin , , 0 asc ¢ ( ˜), 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 T max ) and "level 3" (16× the original T max ) 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 T max 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.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 < f 0 < 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 < f 0 < 400 Hz and 5.8 for 600 Hz < f 0 < 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.Notes.For each range of GW frequencies, this table shows the SFT duration T sft ; the minimum and maximum coherence time T max 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.)a This is the threshold for initiating follow-up, i.e., to produce a level 0 candidate.b This is the number of candidates that would be expected in Gaussian noise, given the number of templates and the follow-up threshold.c This is the actual number of candidates (after clustering) that crossed the S/N threshold and were followed up.d This is the number of candidates remaining after refinement.All of the candidates "missing" at this stage have been removed by the single-detector veto for unknown lines, defined in Section 4. e This is the number of candidates remaining after each has been followed up with a T max equal to 4× the original T max for that candidate.f This is the number of candidates remaining after T max has been increased to 16× its original value.
Note 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 singledetector 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.

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 h 0 2 or h 0 eff 2 ( ) deduced from the highest S/N seen in each band. 307Then, 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 h 0 eff upper limit.As described in Abbott et al. (2017c), to set the adjustment factor for the h 0 limit, we limit attention to injections that were generated with a specified h 0 rather than h 0 eff , of which we recovered 546 of 575, with a resulting adjustment factor of 1.17. 309The upper limits including these adjustment factors are shown in Figure 6.
The 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 r A , at which the disk is truncated by the magnetic field (Pringle & Rees 1972): = ´---- 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 T max ), 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 T max 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. 307We 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). 308These 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 h 0 and 1.21 for h 0 eff .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.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 = 10 8 G and X = 1, gives an Alfvén radius of r A ≈ 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 Figure 6.Upper limits from directed searches in advanced LIGO data.Top: upper limit on h 0 , after marginalizing over NS spin inclination ι, assuming an isotropic prior.The dashed line shows the nominal level balance as a function of GW frequency.Bottom: upper limit on h 0 eff , defined in Equation (2).equivalent the limit on h 0 assuming polarization.(Note that the marginalized upper limit in the top panel is dominated by linear polarization and so is a factor of almost 8 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 h 0 eff levels corresponding to the torque balance h 0 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°± used Zhang al. (2021).Finally, the lighter-shaded band shows the full range of possible h 0 eff 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 of Abbott et (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.
"new" model 2, for which the torque has the form , with the fastness para- . We consider two values for the magnetic field strength at the surface of the star, B = 10 8 and 10 9 G, and consider values of ξ between 0.3 (which sets the upper limit in our plots) and 0.5.
As 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 f 0 ∼ 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., f 0 ≈ 1100 Hz for triaxial emission), and a slower population centered around ν s ≈ 300 Hz (i.e., f 0 ≈ 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 h 0 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 e .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.

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 < f 0 < 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 .(Linear polarization (cos 0 i = ) is at the top, and circular polarization (cos 1 i =  ) 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 = r A ≈ 49 km, which is the value of the Alfvén radius given by Equation (10) for ξ = 1 and B = 10 8 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 model 2, for the range of 0.3 < ξ < 0.5, assuming a magnetic field of 10 8 G (GS8) or 10 9 G (GS9).
consortium.The authors also gratefully acknowledge research  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 e , 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.
max m ˜associated with an offset of 3.3 × 0.011 s.For jobs in which  P max m 0.0625, the initial search was done with P 68023.86 s = ˜and an *A 3 lattice in the other three parameters f (2021), the predictions ofMukherjee 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  f drift | | −T drift model, whose fiducial parameters are taken fromMessenger et al. (2015).

Figure 1 .
Figure1.Search region in terms of parameters t asc ¢ and P ˜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 t asc ¢ and 3.3 × 0.011 s in P ˜, and the truncating boundaries are at ±3 × 185 s of the most likely t asc ¢ 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 t asc and P orb appearing in Table1to t asc ¢ and P ˜.The inner search region, in which we choose a longer T max to do a deeper search, is within ±185 s of the most likely value of t asc ¢ and contains 68.1% of the prior probability, while the full search region, within ±3 × 185 s of the most likely t asc ¢ 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 P ˜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 P ˜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.86s was less than 0.0625.

Figure 2 .
Figure2.The search regions and prior shown in Figure1, expressed in terms of the system parameters t asc ¢ and P orb .For reference, we also show the regions used for the O1 analysis inAbbott et al. (2017c; gray dashed  lines)  and reported for the O2 analysis inZhang 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 forAbbott et al. (2017c) is offset in both P orb and t asc ¢ because it used the P orb estimate inGalloway et al. (2014), while the others used the updated estimate inWang et al. (2018).The analysis ofAbbott 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 t asc ¢ induced by the inaccurate P orb value was less for the epoch of the search (2015 rather than 2019).

Figure 4 .
Figure4.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 Table3); for reference, the thresholds calculated from the numbers of templates in the O1 search ofAbbott et al. (2017c) are shown in orange.Because of the more efficient template placement of algorithm fromWagner 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 < f 0 < 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 < f 0 < 400 Hz and 5.8 for 600 Hz < f 0 < 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 Table3) makes the Gaussian follow-up level an imprecise tool in any event.
(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.

Figure 7 .
Figure7. of upper limits to predictions of torque balance models.The gray band indicates the h 0 upper limit implied by the h 0 eff upper limit in the bottom panel of Figure6, assuming the range of possible inclinations from  0 cos 1 i | |. (Linear polarization (cos 0 i = ) is at the top, and circular polarization (cos 1 i =  ) 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 = r A ≈ 49 km, which is the value of the Alfvén radius given by Equation (10) for ξ = 1 and B = 10 8 G.Note that this is slightly larger than the value of 35 km used in, e.g.,Zhang et al. (2021) andAbbott 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 ofGlampedakis & Suvorov (2021), in model 2, for the range of 0.3 < ξ < 0.5, assuming a magnetic field of 10 8 G (GS8) or 10 9 G (GS9).

Figure 8 .
Figure8.Illustration of how the upper limits shown in Figure6can 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 e , 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.
Wang et al. (2018)ion t asc , 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 inWang et al. (2018)by subtracting P orb /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 P orb , which results in increased uncertainty in t asc and correlations between P orb and t asc ; see Figure2.
d e

Table 2 .
Wagner et al. (2022)jor axis of the orbit is assumed to lie in of[40, 90]km s −1 .The search region in orbital period P orb and time of ascension t asc is constructed using the method ofWagner et al. (2022): the time of ascension t asc is propagated 4125 orbits to define

Table 2
Parameters Used for the Cross-correlation Search

Table 3
Summary of Numbers of Templates and Candidates