{'fraction_non_alphanumeric': 0.03877538864732861, 'fraction_numerical': 0.019030317469771876, 'mean_word_length': 4.324928639391056, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tuning into Scorpius X-1: adapting a continuous gravitational-wave search for a known binary system

We describe how the TwoSpect data analysis method for continuous gravitational waves (GWs) has been tuned for directed sources such as the Low Mass X-ray Binary (LMXB), Scorpius X-1 (Sco X-1). A comparison of five search algorithms generated simulations of the orbital and GW parameters of Sco X-1. Where that comparison focused on relative performance, here the simulations help quantify the sensitivity enhancement and parameter estimation abilities of this directed method, derived from an all-sky search for unknown sources, using doubly Fourier- transformed data. Sensitivity is shown to be enhanced when the source sky location and period are known, because we can run a fully-templated search, bypassing the all-sky hierarchical stage using an incoherent harmonic sum. The GW strain and frequency as well as the projected semi-major axis of the binary system are recovered and uncertainty estimated, for simulated signals that are detected. Upper limits on GW strain are set for undetected signals. Applications to future GW observatory data are discussed. Robust against spin-wandering and computationally tractable despite unknown frequency, this directed search is an important new tool for finding gravitational signals from LMXBs.


Introduction
Continuous gravitational waves (GWs) from neutron stars in binary systems seem likely to be one of the most interesting types detectable by ground-based interferometric observatories. Binary systems constitute 237 of 578 (44%) of known pulsars in the ATNF catalog (v1.53, 2015) [1] with rotational frequency faster than 5 Hz, the approximate lower bound of the frequency range of Advanced LIGO and Advanced Virgo [2,3]. TwoSpect [4,5] is one data analysis method for targeting systems of Low Mass X-ray Binaries (LMXBs), including the prototypical Scorpius X-1 (Sco X-1). This method can apply to continuous-wave (CW) GWs of unknown frequencies from neutron stars in binary systems, with unknown sky locations, orbital periods, or projected semi-major axes. It was used in a prior all-sky search for CWs from unknown sources with data rotates at a frequency ν s , radiating at the quadrupolar GW frequency f = 2ν s , h c ≈ 4 × 10 −27 300 Hz ν s 1/2 F X−ray 10 −8 erg cm −2 s −1 (1) With X-ray flux (3.9 × 10 −7 erg cm −2 s −1 [36]) but unknown frequency for Sco X-1, this limit can be evaluated for f between 50 and 1500 Hz [7], although large systematic uncertainties mean h 0 could be larger or smaller: Angular momentum from accretion would counterbalance that lost to GW emission. If true, f would remain stable for long durations, aside from spin-wandering due to variations in accretion. Computational costs are thus reduced by the restricted range of the frequency-derivative.
The torque-balance levels predicted by Equation 2 are near the estimated thresholds of detection, meriting effort in enhancing the sensitivity of our methods. Directed search improvements skip the computation-saving hierarchical steps of the all-sky search, testing all points in parameter space with templates for enhanced sensitivity to a signal from a known sky location. This paper details the directed, fully-templated method and shows its application to simulated data containing Sco X-1 signals.

Signal model
GW signals are defined by h(t), strain as a function of time. Models of h(t) in a binary search [5] depend on phase evolution, Φ(t), Here solar-system barycentered time is τ (t) and Φ(t = 0) ≡ Φ 0 . We also neglect spinwandering, which could manifest as stochastic variation in Φ 0 . The TwoSpect model does not currently search over amplitude parameters: h 0 is the GW strain amplitude, ι is the inclination angle of the neutron star with respect to the source, ψ is the GW polarization angle. While h 0 can be recovered, its estimation is confounded by ι. Both ι and ψ affect amplitude through detector response, which depends on angle via the plus-and cross antenna functions, F + and F × . Initial GW phase, Φ 0 , further specifies the signal, but is neither explored nor recoverable.
Our search is over unknown Doppler parameters, which drive signal evolution. Sky location (α, δ) is known for Sco X-1, as is orbital period P = 2π/Ω. Initial orbital phase is fixed by time of ascension T asc , to which our method is currently insensitive. Frequency f and projected semi-major axis a sin i (expressed in light-seconds, ls) must be searched over. Because the signal is frequency-modulated by the orbital motion of the source in the circular binary system, the latter parameter manifests through modulation depth ∆f obs = 2πf 0 · (a sin i)/(cP ).

Analysis statistic
TwoSpect has already been described for the all-sky analysis [5,6]. Only a brief summary is given here on the elements common to the directed analysis, with details in Appendix A.
Data containing calibrated GW strain, originally recorded by each observatory as time series, are read into the program from a sequence of short Fourier transforms (SFTs). These SFTs, on the order of minutes long, are shorter than the total observation time T obs , on the order of months. Noise and antenna pattern weights are applied to enhance sensitivity by weighting those SFTs that are more sensitive to a putative source. Each SFT contains K frequency bins. The Earth's motion Doppler shifts the apparent frequency of a source, so these bins must be barycentered: their indices are shifted such that an unmodulated frequency from a given sky location remains in the same bin for all SFTs. These frequency bins represent the instantaneous frequency f of the signal. A second Fourier transform is performed, transforming the power in each frequency bin: the transform is from SFT time to f ′ . Effectively, f ′ represents the orbital frequency of the binary system.
Templates are the M pixel weights w i in the 2-D image-plane (f, f ′ ). Each template corresponds to a signal model from the (f 0 , ∆f obs ) astrophysical parameter space. After the second Fourier transform, the powers Z i of data pixels in the (f, f ′ ) are measured and noise-background λ i estimated. The test statistic, R, is the projection of the data (noise-subtracted powers) onto the templates, normalized by the templates: Estimated p-values for this statistic allow determination of detection probability. Expositions of the detection statistic and signal model can be found in previous methods papers [5,37].

Sensitivity and computational cost
The original design of the all-sky analysis employs hierarchical analysis to control computational costs while still maintaining good sensitivity to a broad parameter space [5]. Each narrow frequency band requires corrections for the Doppler effect caused by Earth's motion. Because this correction depends on sky position, an all-sky search Table 1.
For Sco X-1, distance, eccentricity, X-ray luminosity, sky location, and orbital period are known with good precision and projected semi-major axis is measured to be 1.44 ± 0.18 ls, (Table 1). NS spin frequency, however, is unknown [38]. Other targets, such as XTE J1751-305 [39], have known frequency, reducing computational costs substantially. Sco X-1, XTE J1751-305, and other LMXBs are the principal sources for the new, directed method, described next.
Searching over projected semi-major axis in the range [−3σ a sin i , +3σ a sin i ], around the measured value, and over GW signal frequency f from f min to f max , incurs a predictable computational cost. With a fixed rectangular parameter spacing of 1/(2T coh ) in f and 1/(4T coh ) in ∆f obs , using SFTs of coherence time T coh and analysis bands of width f bw , the search uses N template templates per observatory (Appendix A):

Participation in a Sco X-1 mock data challenge
Several analyses can search for CWs from known neutron stars in binary systems. TwoSpect can also seek unknown systems ‡. The Sco X-1 Mock Data Challenge (MDC) compares five methods [7], including the new, directed version of TwoSpect. While the MDC published the results of the participants, this paper explains how our results were obtained and how similar analyses apply to forthcoming observations.
The MDC simulates 50 open, unblinded and 50 closed, blinded signals, with signal parameters drawn from Sco X-1 astrophysical priors. Each signal is also known as an injection. Three observatories are simulated: H1 (Hanford), L1 (Livingston), and V1 (Virgo). Within 100 frequency bands of 5 Hz each, f is uniformly distributed. Bands cover a range from 50 to 1500 Hz and contain Gaussian noise simulating a detector with noise floor 4 × 10 −24 Hz −1/2 , the expected Advanced LIGO minimum. Projected semi-major axis a sin i is Gaussian-distributed with mean 1.44 ls and standard deviation 0.18 ls. Period P is Gaussian-distributed with mean 68023.70 s and standard deviation 0.04 s. Sky location is (α, δ) = (16:19:55.067, −15 • 38 ′ 25.02 ′′ ). The uncertainty in period is small enough that it is not a free parameter.
Our directed search over f and a sin i, with a fixed period, requires approximately N template = 5.0 × 10 7 templates per observatory, covering 100 bands of 5 Hz each, 500 Hz total. With three observatories, 1.5 × 10 8 templates are needed. Figure 1 shows a wide analysis that is magnified in Figure 2. These analyses require minimal modification from that of the all-sky search: R-statistic and p-value code was identical, aside from looping over all template positions in the parameter grid.

Directed search demonstration and outlier follow-up
The new, directed method is demonstrated through the MDC described in Section 4. This MDC paper [7] presents the performance of each participating pipeline but not the details of our method. Here we highlight methodology not elaborated in the MDC paper, referencing the performance of the directed method for illustration and for contrast with the all-sky method [5].
After writing a loop over templates, post-processing is the main modification to the pipeline: new detection criteria and techniques for parameter and upper limit estimation are required, because the fully-templated output differs significantly from that of hierarchical method. In this section, we detail how these results have been obtained and establish a basis for future applications.

Overview of detection and parameter estimation
A set of extremal p-value outliers in 5 Hz bands is produced for each observatory, subject to a p-value threshold inferred from Gaussian noise. These sets are compared in pairwise coincidence (H1-L1, H1-V1, or L1-V1). Coincidence of outliers allows f or ∆f obs to differ by ≤ 1/T coh between observatories, based on prior experience from the all-sky search [6]. This allowed difference is 2 steps in the f grid or 4 in the ∆f obs grid. Surviving outliers are classified as detections.
For a given detection in one band, the signal parameters are inferred from the values of the template with the highest (single-observatory) p-value. These parameters include f , a sin i, and h 0 . Again, for open signals, the true parameters were known in advance. Total uncertainty in f and a sin i, and non-systematic uncertainty (random) Figure 1. Simulated signal, 100.015 Hz frequency (f ), 1.44 ls projected semi-major axis (modulation depth 0.0133 Hz, ∆f obs ), showing single-template − log 10 p-value. The graph spans 241 templates in ∆f obs , 721 in f , twice as dense in ∆f obs as in f . A signal with h 0 = 4 × 10 −21 is injected into a Gaussian noise amplitude spectral density of 4 × 10 −24 Hz −1/2 , observed for 10 6 seconds. The p-value is extrapolated from the R statistic. For illustration, ∆f obs extends below zero to show that the code is well-behaved and that the algorithm gracefully mirrors results for negative input parameters. Most importantly, − log 10 p is maximal at the true parameters. Signal is partly recovered when template f differs from the true signal by plus or minus ∆f obs , and also when template f and ∆f obs differ equally from the true parameters.
in h 0 , is determined from the standard deviation of the set of parameters of recovered open signals compared to their true parameters, as detailed in section 5.3.
The estimated strain is h rec ∝ R 1/4 . Systematic uncertainty arises from the unknown inclination angle, cos ι, which dominates the total uncertainty in strain. This ambiguity cannot be resolved with the present algorithm and depends partially on the assumed prior distribution of signal amplitudes; the uncertainty is estimated by simulation in Section 5.4. With new algorithm enhancements since the MDC, cos ι can become a searched parameter. Due to the uniform noise floor and low number of injections, a single upper limit value is declared for all 100 bands based on the best estimate of the 95% confidence level of non-detected signals in the open set of injections.
In future applications to real data, the directed analysis can be post-processed using the detection criteria and parameter estimation methods described here. Improved upper limits methods are underway, and a strong candidate signal would likely be followed-up with other analyses, but the core pipeline is the same.

Detection claims
Studying the Gaussian noise in the Sco X-1 MDC open data set, we can set thresholds for detection. In real searches, detection candidates are followed-up, so the same thresholds developed here are used to mark interesting candidates.
To obtain a Gaussian noise sample, we excise injection signals, which are visible in the f vs ∆f obs plane. The excised region depends on injection frequency f inj and modulation depth ∆f inj , on the Earth-orbital Doppler shift (≈ 10 −4 ), and additional bins (at least 10) to avoid spectral leakage. The half-width excised for each injection is δf inj : Hz.
The remaining noise sample is the data set with all intervals In the remaining noise sample, the estimated p-value distribution is not perfectly uniform, due to gaps in the data. Nonetheless, the noise sample provides a distribution of template R-statistics and p-values in the absence of signals. This procedure provides an empirical measure of the estimated p-value that corresponds to an actual false alarm probability of 1% per 5 Hz frequency band. Taken together, these let us establish detection criteria.
If there is any candidate surviving the following criteria in a 5 Hz band, we mark it detected, else not detected: The large discrepancy between the p-value thresholds in the MDC is a historical artifact from a configuration error. The discrepancy is much reduced when this is fixed, as done for future analyses. Our expectation remains that the threshold should be independent of coherence time.
• each candidate must survive at least one double-IFO coincidence test, involving a pairwise comparison of single-IFO candidates to see whether they are within 1/T SFT in both frequency (f ) and modulation depth (∆f obs ).

Parameter estimation and uncertainty for detected signals
Each template is associated with a particular (f, ∆f obs ), so parameters are currently read off from the template with the extremal p-value corresponding to a detection. In the future, accuracy might be improved using interpolation, but the MDC validates that the existing method is highly-accurate. If signals were suspected in real data, this procedure, possibly extended with additional simulations, could generate a parameter space volume for follow-ups to examine. The open set of signals, of which 31 of 50 were detected, are the foundation for understanding parameter estimation uncertainty. The reconstructed h 0 output is h rec = CR 1/4 . The value of C is determined from the mean value of a large number of simulations for circularly-polarized waves over the whole sky and full range of f , P , and ∆f obs .
Then h rec is rescaled twice, first by ρ R for more accurate measurement at the Sco X-1 period and modulation depth, and second by ρ cos ι for unknown cos ι. Thus the final claimed value of h 0 for a signal j is (h 0 ) j : The first scale factor, ρ R = 1.11, corrects the average values of h rec (h 0 -reconstructed) in the open set to match the corresponding h eff (h 0 -effective, given circular polarization weightings). That is,h eff = 1.11 ×h rec , where h eff is defined a priori by h inj (h 0 -injected) and cos ι: In the MDC study [7], 4 of the 31 detected, open signals account for the largest frequency estimation error. This error arose from a misconfiguration that does not affect the other analyses; it was addressed by taking those 4 as one class and the remaining 27 as another, a step that should be unnecessary in future analyses. Then, the uncertainty due to random error for h 0 , f , and a sin i is estimated by the standard deviation between the recovered and true parameters: where N open is the number of open injections, σ f , σ a sin i , and σ h 0 are the uncertainties we state for recovered f rec , (a sin i) rec , and h rec given injected f inj , (a sin i) inj , and h inj . The error between injected and recovered parameters does not show any other clear correlation with p-value or signal frequency, at least in the 31 detected signals. Except for the most marginally detected signals, where noise fluctuations matter, uncertainty in f and a sin i is dominated by the template grid spacing. The σ f and σ a sin i error bars have been used uniformly for claiming uncertainties on the signals in the MDC.
The largest source of uncertainty for h 0 comes from correction for systematic underestimation, multiplying a factor of 1.74 into C. This uncertainty is the ambiguity in cos ι discussed in Section 5.4. Parameter estimation uncertainty for f and a sin i is then just the random error; for h 0 , it is the quadrature sum of random error and cos ι ambiguity.

Ambiguity from cos ι
The largest systematic uncertainty in h 0 comes from the unknown cos ι. The method is optimized for | cos ι| = 1 and computes the R statistic by weighting the SFTs assuming circularly-polarized GWs, which still provides good sensitivity for other polarizations. Recall that h rec must be scaled by ρ R = 1.11 to match h eff . When a source is circularly polarized, the analysis estimates ρ R h rec ≈ h inj . In the case of linear polarization, Equation 9 indicates that h inj will be about 2 3/2 times larger than h eff (and so h rec ). The aim is to find an average conversion factor from h rec to h inj and a robust estimate of the uncertainty.
Since h eff of circularly-polarized signals is greater, those signals are more easily detectable than linearly-polarized signals of equivalent h inj . Therefore the signals that are detectable are biased, near threshold, to being more likely circularly polarized. This "circularizes" the correction factor, depending on the detection efficiency of the pipeline and on the assumed prior distribution of strain amplitudes. Although the effect is minor, estimating its size requires simulation.
The simulation generates 2 million signal amplitudes between h inj = 3 × 10 −26 and h inj = 3 × 10 −24 with a distribution of 1/h 0 , the assumed prior distribution of h 0 values. This simulation code is independent of TwoSpect and should apply to similar directed searches. In this simulation, ρ R is 1 for simplicity, so we can treat h rec = h eff here. We model detection efficiency by assuming no signals are detected below h eff = 1×10 −25 , all are detected above h eff = 3×10 −25 , and the fraction detected is linear in h 0 between those values. Together with a uniform cos ι distribution of [−1, 1], this leads to a trapezoidal distribution of recovered, detected h 0 values with a curved lower (left) edge ( Figure 3).
Part of the domain of the simulation must be excluded. To find the average h inj /h rec for a given h rec , every h rec must correspond to a full sampling of the range of polarizations. A large h inj with linear polarization or small h inj with circular polarization could have the same h rec . No linearly-polarized signals could produce h rec above 1 × 10 −24 , because Equation 9 shows that the largest signal, h inj = 3 × 10 −24 , would be reconstructed a factor of 2 √ 2 smaller, at h rec = 1.06 × 10 −24 (again, ρ R = 1 for the simulation). Above 1 × 10 −24 , polarizations tend to be more circular, thus the average ratio must exclude this region or it will be biased by the limited range of the simulation h inj . Expanding the domain would raise the cutoff, although the resulting ratio would no longer perfectly correspond to the MDC.
With the domain of the simulation determined, we compute that mean ratio of h inj to h rec to be 1.74. Fine-binning h rec , an interval of [1.74 − σ cos ι , 1.74 + σ cos ι ] encloses 68% of corresponding h inj when σ cos ι is 0.37 (found by manual optimization). Therefore our best estimate for the correction factor ρ cos ι , with σ cos ι inferred as the standard deviation, is 1.74 ± 0.37. This factor multiplies ρ R , which is found to be 1.11.
The systematic uncertainty, being the uncertainty σ cos ι in the correction factor, scales with signal strength; the non-systematic (random) is fixed and is also multiplied by the correction factor. The final estimate of the uncertainty in h 0 for signal j is the quadrature sum of the systematic and non-systematic uncertainties: For future data, a similar simulation could be run, with an updated detection efficiency model and prior distribution of strains, to find the uncertainty in h 0 due to cos ι for promising signals.

Accuracy of parameter estimation uncertainty claims
The scheme described above reliably recovers parameters and states uncertainties consistent with the true distribution of errors, as shown in the MDC [7]. Verifying the calibration factors and confidence intervals once more, one can confirm that a conservative fraction of h 0 , f and a sin i are within their 1-σ error bars: 77.4% for σ h 0 , 74.2% for σ f , and 67.7% for σ a sin i .

Upper limits for undetected signals and detection efficiency
Only upper limits have come from CW searches to date. Until GWs are detected from neutron stars, the scientific value of a CW analysis is to constrain the plausible h 0 and inferred ellipticity from those stars. For the MDC-simulated signals, we simplified upper-limit estimation. In a real detector, noise varies with frequency [43]; in this simulation, the noise floor is flat at 4 × 10 −24 Hz −1/2 . Given observational data, we would inject a large number of simulated signals into a number of smaller bands, in order to understand the upper limit as a function of frequency.
To measure detection efficiency, we calculate the h eff for all signals and find the average detection rate for a given h eff . Binomial uncertainty is also calculated and each 1-σ deviation (per 5-signal bin) is graphed in Figure 4, which shows a least-squares sigmoid fit. This detection efficiency curve maps from strain to probability of detection. Next we would like an upper limit function that takes a probability as an input and returns a strain that, with the given probability, is no less than the actual strain.
To characterize the upper limit, we plot the distribution of h rec versus injected h eff in Figure 5. We verify that 95% of non-detected open signals are covered by a naive upper limit of h rec = 2.19 × 10 −25 . This claim does not rely on binning but rather on the sampling density of the injections. This number, when corrected by the ρ R rescaling factor of 1.11 and cos ι correction ρ cos ι of 1.74, yields the upper limit of (1.74) × (1.11) × 2.19 × 10 −25 = 4.23 × 10 −25 . Because of the flat noise floor, this is reported as a single upper limit for any non-detections in the MDC.
Frequency-dependent upper limits are well under development for actual observations. Sigmoid fits to the detection efficiency of a set of injections into real data are generated, one fit per frequency band. Upper limits at a given confidence can then be taken as the h 0 that yields a detection efficiency equal to that confidence. This advancement post-dates the MDC and is planned for future applications.

MDC results
TwoSpect analyses applied to the MDC data set detect more stars than the Radiometer, Sideband, or Polynomial pipelines; only the CrossCorr algorithm found more signals [7]. Each detection includes an estimate of a sin i, which is not produced by Radiometer, Sideband, or Polynomial. The MDC did not model the spin-wandering of the neutron star that is expected in real data, although participants were told to assume its presence, and spin-wandering is planned for future MDCs. TwoSpect is also theoretically highly robust against spin-wandering. This method has already been applied to real data [6], though not using the directed search in a fully-templated mode. This experience validates the program's robustness with respect to non-Gaussian data artifacts. In all, 34   determined for the 16 non-detected, closed signals. Although the distribution of h 0 values in the MDC was astrophysically optimistic, the MDC validated our ability to claim detections and recover orbital and GW parameters accurately.

Future directed CW binary searches
Algorithms such as TwoSpect are designed to find astrophysically-plausible strain from LMXBs. Torque-balance arguments suggest that strain could exceed 10 −25 for Sco X-1 if it rotates at low frequencies.
The previously-published all-sky search in a year of S6 data set an overall upper The unity-slope black line is shown to ascertain whether a further empirical rescaling factor is needed to match the dashed-and-dotted least-squares linear fit (it is: constant 1.11). The zero-slope horizontal dashed line is shown to indicate the ninety-five percent confidence upper limit in the absence of detection.
limit for circular polarization of 2.3 × 10 −24 at 217 Hz [6], set for an H1 amplitude spectral density of 2.0 × 10 −23 Hz −1/2 . For random polarization, a multiplicative factor of 2.2 was applied, for an upper limit of 5.1 × 10 −24 . This corresponds to a sensitivity depth (factor below the noise floor: Appendix A) of 8.7 Hz −1/2 for circular polarization but 4.0 Hz −1/2 for random. Extrapolating to Advanced LIGO design sensitivities of 4 × 10 −24 Hz −1/2 , this implies an upper limit around 5 × 10 −25 for circular polarization and 1 × 10 −24 for random. However, the all-sky paper included an opportunistic search for Sco X-1 on a narrow frequency range (20 to 57.25 Hz), setting a random polarization limit of 2 × 10 −23 at 57 Hz in an L1 amplitude spectral density 1.8 × 10 −22 Hz −1/2 , a depth of 9 Hz −1/2 . This opportunistic search used 1800-s SFTs; longer SFT durations increase theoretical sensitivity. Significant improvement comes from focusing on one sky location; this can be viewed as a reduced trials factor. The directed search demonstrated in this paper achieved an 4.23×10 −25 upper limit in simulated data at design sensitivity, a depth of 9.5 Hz −1/2 . It achieved this depth despite being tested with shorter (360-s and 840-s) SFTs. The directed search is also scalable over a much wider parameter range, like the all-sky method over which it gains twofold in sensitivity. While real data complications may worsen this limit, several simplified and conservative steps were taken. The limit may improve with the enhancements now under developement, when fully tested with injections as in the all-sky search. Even now, the method in this paper is more sensitive for random polarization than the all-sky method is for optimal, circular polarization. Additional improvements to the algorithm, such as coherent SFT summing, have been developed [37] and could further improve this limit in the future, pushing toward the torque-balance strain.
Directed TwoSpect analyses have been demonstrated in this paper. Comprehensively covering the parameter space of Sco X-1 at full sensitivity with the directed search, instead of hierarchically as before, does increase the probability of detection and improve upper limits. When detections do occur, the ability to determine the frequency and projected semi-major axis of the neutron star in the binary system will prove highly informative. Analyses of real data for signals from Sco X-1 and additional neutron stars in binary systems, such as XTE J1751-305, are underway. In the long term, we hope that the discovery of gravitational waves from neutron stars in LMXBs will provide a firm link between our observations and electromagnetic astronomy.

Equation 6b
is the integrated template density over parameter space dimensions. While we are interested astrophysically in a sin i, the observable ∆f obs governs template placement. The equation decomposes into two number densities, Q frequency and Q mod , and corresponding dimension length intervals, L frequency , L mod : where Q frequency is the inverse of the template spacing in frequency, which is 1/(2T coh ), so Q frequency = 2T coh . Index the frequency dimension by j bands. One step j is made per band, width f bw , so the length interval L frequency = f bw . Thus With fixed frequency bands, the number of templates per band does not change. The width in modulation depth, however, depends on a sin i and the frequency f j = (f min + jf bw ). Each template of ∆f obs is indexed by k. Analogous to before, when a sin i is given in light-seconds. Substituting Equation A.3 into Equation A.2, that term depends on neither the frequency nor modulation depth index, and so pulls out in front of the sums. Equation A.5 depends on k; the sum Σ (mod) is evaluated from (a sin i) k = (a sin i − 3σ a sin i ) to (a sin i + 3σ a sin i ) in practice, for an integrated length of 6σ a sin i . Combining these elements, including indexing frequency band steps by j, Writing the limits of the sum of frequency bands, in addition to inserting ones to account for the edges of each sum, yields Equation 6b.

Appendix A.3. Test statistic calculation
The construction of this R-statistic can be described in several steps. The most important points from the original methods paper [5] are reiterated here, with some clarification. Science/observing runs are first parcelled into overlapping short Fourier transforms (SFTs), performed in the detector frame. The SFTs have typical coherence time T coh (referred to as T SFT in newer publications [37]) ranging from 60 s to 1800 s, depending on the hypothesized time-derivative of neutron star frequency [5]. The total number of SFTs with 50%-overlap for an observing time T obs is N, SFT number in the observing run are indexed by n ∈ [0, . . . , N − 1]; SFT frequency bin is indexed by k ∈ [0, . . . , K − 1], where K = T coh f N for a Nyquist frequency f N and only positive frequencies k = T coh f are used. Thus the transformation from time series to SFTs is a map from x(t) (= h(t) + n(t), signal strain plus noise) tox n k . This array of n still depends on detector time t, and the analysis is to be done in Solar System Barycenter (SSB) time t SSB . Travel from the source to SSB introduces an overall phase shift; uncertainty in the distance and proper motion is systematic and the same for gravitational and electromagnetic observations. Detector time is recorded in GPS time, running parallel with Terrestrial Time (TT), and SSB time runs parallel with Barycentric Dynamical Time (TDB). SSB time corrects t by δt R for relativistic effects. Another overall phase shift is caused by Roemer delay ∆t R , the dot product ofn/c from the SSB to the sky location of interest with r from the SSB to the detector [5,12,37]. Barycentering detector-frame data is equivalent to resampling in τ (t) = t SSB (t)+∆t R (t), Each SFT frequency binx n k is Doppler shifted to a frequency bin in the SSB frame corresponding to the sky location, frequency, and time t of the midpoint of the SFT n under investigation: k(f (α, δ, t)) → k(f SSB (τ )). This barycentering procedure corresponds to the time-domain Equation A.8. Henceforth, barycentering is implicit in the k index.
Define the power P n k = 2|x n k | 2 /T SFT . Let P k n be the expected (estimated from a running mean over nearby n) noise-only power in a frequency bin k for SFT n. Also let F 2 n ≡ F 2 +,n + F 2 ×,n for the antenna pattern at the chosen sky location and SFT n -taking this equal-weighted sum of F + and F × polarization components implies an assumption of circular polarization. Then the estimated power in a given binP n k is normalized such that random, white, Gaussian noise will have an expectation value of 1 [5]: Then each row of barycentered frequency bins k is treated as a time series in n. Power for bin k in that time series is Fourier transformed by F f ′ into Z k (f ′ ), where f ′ is the second Fourier transform frequency. During the transform, the background noise power λ(f ′ ) is estimated from the noise in the SFTs, assuming the noise is Gaussian. This second Fourier power Z k (f ′ ) follows a χ 2 distribution with 2 degrees of freedom and mean 1.0, is proportional to h 4 , and is constructed by When re-indexed by sorted template weight, i, Z k (f ′ ) becomes Z i in Equation 5. Extensive discourse on the details of these calculations, as well as the estimation of background and calculation of template weights, is found is the original TwoSpect methods paper [5], and a paper on coherent addition of SFTs rigorously derives the SFT power by including Dirichlet kernel terms. Given the second Fourier power, we calculate the R-statistic its p-value, from which we seek to make a detection.