Search methods for continuous gravitational-wave signals from unknown sources in the advanced-detector era

Continuous gravitational waves are long-lasting forms of gravitational radiation produced by persistent quadrupolar variations of matter. Standard expected sources for ground-based interferometric detectors are neutron stars presenting non-axisymmetries such as crustal deformations, r-modes or free precession. More exotic sources could include decaying ultralight boson clouds around spinning black holes. A rich suite of data-analysis methods spanning a wide bracket of thresholds between sensitivity and computational efficiency has been developed during the last decades to search for these signals. In this work, we review the current state of searches for continuous gravitational waves using ground-based interferometer data, focusing on searches for unknown sources. These searches typically consist of a main stage followed by several post-processing steps to rule out outliers produced by detector noise. So far, no continuous gravitational wave signal has been confidently detected, although tighter upper limits are placed as detectors and search methods are further developed.


Introduction
The search for continuous gravitational-wave signals (CWs), long-duration forms of gravitational radiation, is one of the endeavours of gravitational-wave astronomy. These signals are produced by long-standing quadrupolar variations, such as rapidly-spinning neutron star (NSs) sustaining a crustal deformation, undergoing an r-mode instability or in free precession [1][2][3][4], as well as more exotic sources such as the annihilation of ultra-light boson clouds around spinning black holes [5][6][7][8] or compact dark matter objects (CDOs) in the Solar System [9].
Detecting a CW signal could shed some light on NS physics, as well as open a new channel to test general relativity (extra polarizations, Lorentz violations) or detect dark matter [10][11][12][13][14]. No confident CW detection has been reported to date. The current product of CW searches are source-agnostic upper limits on the nominal CW amplitude h 0 . These results can then be mapped into different astrophysical scenarios, such as the ellipticity of nearby NSs, the mass of ultralight bosons around black holes [15,16], or the nearby population of planetary-mass primordial black hole binaries [17,18].
The present document reviews search methods and pipelines employed to look for CW signals in the observing runs performed by the second generation of ground-based interferometric detectors (advanced detectors). Reviews on the physical mechanisms of CW emission by NSs can be found in [1][2][3][4]. Basic data analysis techniques are discussed in [19]. Finally, [20] discusses the main results of previous CW searches up to 2017.
The standard CW signal model consists of a quasi-monochromatic source emitting gravitational waves at a certain frequency f 0 . For the case of a NS sustaining a certain ellipticity, f 0 corresponds to twice the spinning frequency of the star. Further time derivatives of the frequency arise due to different physical mechanisms affecting the source, such as arXiv:2111.12575v2 [gr-qc] 7 Dec 2021 energy emission as gravitational or electromagnetic radiation. In the case of sources in binary systems, the orbital motion induces a Doppler modulation.
Regardless of the specific intrinsic frequency modulation of the source, CW signals as seen from Earth are Doppler-modulated due to the detector motion around the Solar System barycenter (SSB). For a source with a given intrinsic frequency evolutionf , the detector-frame frequency is given by where the phase-evolution parameters λ include the CW frequency f 0 and the sky position of the source n, as well as any other parameter describing the intrinsic frequency evolution of the signal. v(t)/c refers to the detector velocity expressed as a fraction of the speed of light, and contains both the daily and yearly motion of Earth around the SSB, of orders O(10 −6 ) and O(10 −4 ), respectively. The yearly modulations can be clearly seen in the left panel of Fig.  1; the daily modulation is contained within a frequency bin. The amplitude of a CW signal is described using four parameters, namely the initial CW phase φ 0 , the sperical angles describing the orientation of the source {ψ, cos ι}, and the nominal gravitational wave amplitude h 0 . The response of a ground-based detector to a passing CW is better described in terms of the so called JKS representation [21] h(t; λ, A) = where the four time-independent A µ depend on the four amplitude parameters A = {φ 0 , ψ, cos ι, h 0 } and the antenna-pattern response of the detector is contained in the four quadratures h µ (t; λ), which only depend on the phase-evolution parameters. The basic effect of the antenna-pattern response on a CW signal is a daily amplitude modulation, clearly visible in Fig. 1.
As opposed to the short signals produced by compact binary coalescences (CBCs), with typical durations between a few minutes and less than a second for current detectors, CW signals are expected to last for years, spanning several observing runs of the current and future generation of ground-based interferometric detectors [22][23][24][25][26]. This difference in duration is crucial in terms of detecting and estimating the parameters of a CW signal.
Modelled searches are usually performed using matched filtering [27], comparing the datastream to a set of templates in order to find a high correlation. Due to the typical duration of a CW signal (spanning the entire observing run), the required number of templates to perform a blind search is prohibitively high even for current computing standards [28][29][30].
The standard strategy, in a broad sense, is to reduce the effective length of the datastream by performing matched filtering over shorter segments; the segment-wise results can then be combined into a final statistic. These kind of schemes are usually referred to as semicoherent searches [19]: the segment-wise analysis is typically referred to as coherent, as it compares the phase evolution of a signal with the datastream throughout a coherence time T coh . The resulting coherent filters are then combined incoherently (ignoring phase information) into a final detection statistic. This incoherent combination allows to recover part of the sensitivity lost due to the split of the initial datastream; the final sensitivity, however, is lower than that of a fully coherent search. Figure 2 illustrates the principle of operation of semicoherent searches. By comparing shorter streams of data, a looser constraint is imposed when comparing phaseevolution templates. As a result, the characteristic width of detection-statistic peaks widens, reducing the required number of templates to ensure a good covering of the parameter space. This strategy is at the core of multi-stage approaches such as those discussed in Sec. 4.
Parameter estimation, on the other hand, is positively affected by the long signal durations. Typical frequency resolutions are under a mHz for initial stages, achieving nHz  Effect of different values of T coh on one of the standard detection statistics for CW searches, the F -statistic (see Sec. 2.1). Longer coherence times impose a greater penalization to deviations with respect to the signal model; as a result, peaks tend to become narrower around the true signal parameters. resolution at fully-coherent follow-ups using a year of data. Sky localization is also significantly improved. If we think of CBC sky localization, neglecting contributions from anntenna-pattern amplitude modulations and higher GW modes, the problem is basically that of determining the direction a GW pulse came from, which can be solved by means of measuring the pulse from N detectors and finding the overlapping sky-positions [31]. The case of CW signals is much simpler, as they are not just single pulses, but continuously arrive at the detector as it moves around the SSB. Due to this movement, a single interferometric detector recieving a CW signal at different positions with respect to the SSB is essentially equivalent to an arbitrary large set of different detectors reciving the same pulse from a source. Hence, for CW sources, sharp sky localization can be achieved using a single detector by simply extending the duration of an observing run.
CW searches require a very fine parameter-space resolution in order not to miss a signal, increasing the computing cost of a matched-filtered search up to unaffordable figures [28]. As a result, searches for CW signals from unknown sources tend to follow a hierarchical approach [32][33][34][35][36][37][38]: wide parameter-space regions are analized using a less-constraining statistic so that a coarser template bank can be used, see Fig. 2. Interesting regions are then typically small enough to follow up using a more sensitive statistic. A quantitative description of this strategy can be found in [39].
The structure of this work goes as follow: in Sec. 2 we review the main search methods and pipelines employed to search for CW signals from unknown sources during the era of the advanced detectors. Section 3 discusses different post-processing stages employed by said searches; these include both specific prescriptions to select interesting candidates (e.g. clusterings) and consistency arguments to assess the overlap of a specific set of templates with instrumental disturbances. Section 4 reviews follow-up strategies to further analyze interesting candiadtes. So far, no search has claimed a confident CW detection, reporting instead constraints on the maximum detectable amplitude achieved. Different approaches to construct said constraints are listed in Sec. 5. Finally, a summary of the reviewed methods is presented in Sec. 6.

F -statistic searches
The F -statistic is a standard detection statistic for CW signals. Initially derived as a maximum-likelihood estimator with respect to amplitude parameters [21,70], it was later rederived in a Bayesian context as a Bayes factor, gauging the presence (or lack) of a signal in a Gaussian noise data stream, in which amplitude parameters are marginalized using a rather unphysical set of amplitude priors [71][72][73][74][75][76]. This detection statistic can be extended for more generic types of sources, such as binary white-dwarf systems [77] or the inspiral phase of binary black-hole coalescenses [78].
Semicoherent searches balance sensitivity and computing cost by choosing a suitable number of coherent segments whose combination into a semicoherent quantity can be performed in an efficient manner. The basic "stack-slide" procedure used for many F -statistic semicoherent searches [32][33][34]79] is to set up a template bank in each coherent segment to compute the segment-wise coherent detection statistic. Then, in the semicoherent stage, the semicoherent detection statistic is computed on a finer template bank by combining results from the segment-wise coherent detection statistics. This template bank refinement can be seen as a consequence of using multiple coherent segments: the higher the number of coherent templates to be combined, the higher the resulting number of distinct semicoherent segments. Discusion on the specifics of template bank refinements are deferred to Secs. 2.1.1 and 2.1.2.
We discuss three different implementations of semicoherent F -statistic searches, namely the Global Correlation Transform hierarchical search (GCT), Weave and Time-domain F -statistic. They mainly diverge in the manner of constructing semicoherent quantities, taking different trade-offs in terms of computing cost, memory requirement and robustness to non-Gaussianities. We note, however, that fully-coherent searches for targets at a specific sky position, such as supernova remnants [59,66], are still performed nowadays.

GCT hierarchical search
Detection statistics (and the F -statistic in particular) present a set of characteristic correlations across the CW parameter-space due to some level of degeneracy present in CW signals [80]. Understanding the structure of said correlations offers a simple method to reuse coherent F -statistic values to construct a semicoherent statistic, reducing the overall computing cost of the search. The GCT introduces a set of coordinates defined by the intersection of parameter-space correlation surfaces to identify nearby parameter-space points where the F -statistic achieves a high value due to the presence of a signal [81][82][83].
Due to its interesting trade-off between sensitivity and computational cost [84], semicoherent GCT searches have been used by the Einstein@Home project to perform deep searches throughout different observing runs [40,49,61,62,69]. Einstein@Home searches distribute the computational load of a search across a volunteer-computing network using BOINC [85] to analyze a high number of parameter-space candidates.
At the core of the GCT search there are combinations of F -statistic values computed at different coherent segments containing data from one or multiple detectors. These are combined into the semicoherent F -statistic using a finer template bank with refinement only on the spindown parameters [82]. This allows for the easy implementation of Bayesian extensions over the semicoherent F -statistic, such as the "line-robust" B S/GL or B S/GLtL statistics [86][87][88][89], which combine F -statistics from different detectors to suppress singledetector artifacts.
No general computing model is available for the GCT search, requiring extensive softwareinjection campaigns in order to numerically tune the sensitivity of a search to the available computing resources [36]. As discussed in Sec. IV B of [90], this is due to a core assumption in [32,82,83] neglecting refinements in the sky parameter-space. Further developments in the field, discussed in Sec. 2.1.2, proposed a new strategy to solve this problem. For the case of directed searches, however, optimal setup strategies are available [91].

Weave
Local parameter-space structure can be understood in terms of the mismatch µ [92,93]. Given a signal parameterized by a set of parameters λ 0 , the parameter-space mismatch µ quantifies the fractional loss of (squared) signal-to-noise ratio ρ 2 produced by an arbitrary parameter-space off-set ∆λ with respect to the true signal parameters: In a neighbourhood of λ 0 , the mismatch can be expressed in terms of a quadratic form as where the symmetric 2-rank tensor g plays the role of a Riemannian metric in the parameterspace.
It was quickly realized that, for the case of flat parameter spaces, the metric g offers a complete description of the parameter-space, allowing for the easy construction of optimal setups [29,94]; alas, the standard CW parameter-space (specifically, the sky-position subspace) presents a non-trivial structure resulting in a curved parameter-space.
Weave [90] represents the coming together of a number of search strategies. It implements, for the first time, a semicoherent search with a well-understood computational model based on a flat parameter-space metric. The setup is capable of constructing template banks at the suitable resolution to achieve an optimal computing cost [30,90,95,96]. Setting up Weave requires only a list of time stamps delimiting semicoherent segments, a coherent mismatch µ, with which single-segment template banks will be constructed, and a semicoherent mismatchμ, to setup the semicoherent template bank. The identification of a template from the semicoherent bank to its corresponding ensemble of coherent templates is handled by the Weave code using the coherent parameter-space metric to identify the nearest neighbor in each segment. In a sense, Weave retains the general characteristics of the GCT search, mainly being an engine to combine coherent F -statistics, but, as opposed to it, Weave requires no extensive numerical calibration to be deployed using an optimal setup. The success of Weave as an all-sky search, however, is related to the characteristics of the F -statistic's parameter-space structure. As discussed in depth in [96,97], the optimal setup for a realistic computing budget tends to yield semicoherent mismatch values well beyond the validity of the metric approximation (i. e. µ 1). Nevertheless, empirical studies [96] show that in this regime the F -statistic actually falls off more slowly than predicted by the metric, meaning the resulting template bank will contain more templates than strictly required. An alternative approach circumventing the empirical characterization of the µ 1 regime is discussed in [94,98,99].
Despite achieving a better sensitivity than the GCT search at a fixed computing cost, the increased memory requirements of Weave make it so far unsuitable for its deployment on Einstein@Home [84]. Nonetheless, its sensitivity and setup flexibility have already been proved in the implementation of a novel CW search strategy [51].

Time-domain F -statistic
Both GCT and Weave searches rely on a common implementation of the F -statistic [100], publicly available under LALSuite [101]. The Time-domain F -statistic pipeline [102] uses a different implementation, based on the F -statistic's time-dependent behavior throughout the observing run. Instead of computing semicoherent quantities, it focuses on significant parameter-space points at coherent-segment level, looking for coincident candidates across different time segments and detectors. Using this coincidence criterion (which we review in more depth in Sec. 3) the pipeline automatically becomes robust to strong instrumental features, since they tend to overlap with different parameter-space regions as an observing run progresses.
This particular implementation of the F -statistic uses its own template bank setup in order to optimize the number of fast Fourier transform (FFT) computations [103,104], which normally takes a significant part of the overall computing cost of the search. Further improvements at parameter-estimation level are optimization algorithms to resolve the characteristic frequency multi-modality of the F -statistic [105] and the inclusion of machine-learning algorithms to filter out non-astrophysical candidates [106].

Fourier-transform-based searches
The efficiency of F -statistic searches stems from the marginalization with respect to amplitude parameters, removing four parameter-space dimensions from the search [29]. Semicoherent pipelines, moreover, assume amplitude parameters to be independent across different coherent segments. This approach makes it difficult to search for amplitude parameters, such as specific CW polarizations.
An alternative family of methods use Fourier transforms of short data segments (Short Fourier Transforms, SFTs) as the basic unit of operation. The duration of an SFT is typically such that CW signals are contained within a single frequency bin. For an isolated NS, this length is about 30 minutes [107], although the exact value depends on the considered frequency range. We note that SFTs themselves can, in some situations, corresponds to coherent segments of a semicoherent search; in such cases, the effective coherent length is proportional to the SFT length.
In order to construct a detection statistic S, Fourier transform amplitudesx are combined using a set of weights (effectively a kernel) K taking into account source polarization, antennapattern amplitude modulations, Doppler modulations and relative phase deviations across different detectors: where parameters other than time dependency have been kept implicit for the sake of simplicity. Frequency modulations are assumed to be contained within the specific set of Fourier amplitudes being combined, although K can be configured to increase robustness against different kinds of spectral leakage by including neighboring frequency bins into the kernel [39,108]. We review two families of searches stemming from Eq. (5), depending on whether their primary target is to increase the robustness against deviations from the intended CW model (PowerFlux & Falcon) or to improve sensitivity by increasing the effective amount of data used (Cross-Correlation).

PowerFlux & Falcon
The PowerFlux pipeline [109][110][111][112] estimates power from a CW source depending on its sky position. The basic implementation [41,42,54] uses a diagonal kernel K to combine Fourier power from each SFT. In this case, the role of K is to diminish the contribution of unfavourable frequency bins (due to high noise floors or a low antenna-pattern response) to the total weighted Fourier power. Since the final statistic ignores relative phase shifts between SFTs, this corresponds to a semicoherent search for a specific CW polarization.
Loosely coherent methods [39,113,114] exploit the flexibility of Eq. (5) to set up a kernel K to account for unmodeled phase shifts, be it due to parameter-space mismatches or unaccounted physics such as small binary orbital modulations [115]. To do so, phases are allowed to drift at most by a specific amount δ across contiguous data segments, obtaining as a result [39] where data is assumed to come from a single detector for the sake of simplicity. (The framework presented in [39] is flexible enough to treat data from multiple detectors.) This simple kernel illustrates the principle of operation of loosely coherent searches: if δ = π, K behaves like a delta function and phases among contiguous segments are uncorrelated, effectively performing a semicoherent search; if δ = 0, K = 1 and phases are correlated throughout the full data stream, performing a fully-coherent search. Tuning δ to intermediate value allows to trade sensitivity and computing cost: low δ values involve simpler kernels (less non-zero terms), easing the computing cost of Eq. (5) but imposing tighter constraints with respect to the chosen CW model. A discussion on physically relevant δ values can be found in [39], whereas an efficient implementation of low-δ kernels can be found in [113,114]. Due to its high computing cost, the initial implementation of loosely coherent methods on PowerFlux was mainly used on directed searches, such as spot-light surveys [55] or follow-up stages (see Sec. 4). An efficient implementation, Falcon, was later developed to perform all-sky searches throughout O1 [43,44] and O2 [46][47][48], setting competitive constraints on the nearby population of Galactic neutron stars.

Cross-Correlation
The Cross-Correlation search [116][117][118], implemented in the LALSuite library [101], closes the sensitivity gap between semicoherent and coherent searches by increasing the effective amount of data using the correlation of Fourier amplitudes at different times.
Although initially designed to look for stochastic GW backgrounds, CW signals are a perfect target to be looked for using Cross-Correlation, as they are long-lasting and deterministic. This allows to use not only the cross-correlation of different instruments during a certain period of time, but also the cross-correlation of data segments at different periods of time. In the context of CW searches, Cross-Correlation searches have been employed to look for specific high-priority targets, such as CW emission from the LMXB system Scorpius X-1 [58,67]. Complementarily, searches for stochastic GW backgrounds, such as the Radiometer search [119][120][121], are able also to constrain CW amplitudes from specific sky locations, such as those of SNRs, LMXBs such as Scorpius X-1 or the Galactic Center [122][123][124].
To do so, the kernel K is constructed following the expected correlation of a signal across different stretches of data at different times and detectors. Specifically, only amplitudes within a certain time range T max are combined together. As opposed to PowerFlux and Falcon, however, polarization angles are averaged out using uniform priors. The T max parameter plays a similar role to that of δ (see Sec. 2.2.1) in terms of trading computing cost and sensitivity: longer T max increases the number of cross-correlations to perform, but also imposes a tighter constraint with respect to the specified CW model. Latest developments on this pipeline include the use of re-sampling techniques to accelerate its evaluation [125] and the use of shear parameter-space transform to optimize the setup of template banks [118].

Hough-transform semicoherent searches
Loud instrumental features in the data tend to saturate detection statistics due to their strong resemblance to CW signals along short periods of time. This problem pushed forward the development of methods capable of suppressing narrow-band features in the data, such as the application of the Hough transform to the search for CW signals.
The basic idea is to limit the contribution of each semicoherent segment to a bounded quantity, which contributes by a limited amount to the overall statistic. This is done by binarizing a power-like quantity (whitened Fourier power formatted as SFTs, as discussed in Sec. 2.2, or the F -statistic) into ones and zeroes using a predetermined threshold [107]. Such a binarization, however, comes at the cost of ignoring noise-floor variations and antennapattern amplitude modulations, meaning highly-contaminated frequency bands around lowsensitivity sky positions will contribute the same amount as clean frequency bands at the most favored sky positions. This problem is usually solved introducing a set of weights into the detection statistic [126,127].
The Hough transform [128] can be used to identify shapes in binarized images. Given a parameter-space parametrizing a family of curves, each set of parameters is assigned a score, the number count, proportional to the number of pixels in the image consistent with the corresponding curve. In CW searches, the parameter space is typically the set of phaseevolution parameters, according to which the data spectrogram is traversed adding (weighted) ones and zeroes depending on whether each frequency bin contains excess power or not.
Hough-transform based searches have been widely used during the latest all-sky searches [41,42,45,50,52,53]. We discuss the two main implementations of the Hough transform for CW searches, SkyHough and FrequencyHough.

SkyHough
The SkyHough pipeline was the initial implementation of the Hough transform to the search for CW signals [107,129]. Due to the Earth's movement around the SSB, a CW signal arriving at the detector at a certain time t with a given frequency f 0 has a very specific set of sky positions from which it could have originated. Said sky positions take the shape of thick annuli ("circles in the sky") [80,107,130] which, for consistency, can be described themselves as one/zero regions in the sky patch. The binarization of the data spectrogram, thus, gets mapped into the selection and summation of sky patches containing different selected annuli.
The "circles in the sky" are weakly affected by local changes in f 0 ; hence, once computed for a specific frequency, these structures can be saved into a look-up table (LUT) to analyze several frequency bins. The use of LUTs, on the other hand, affects the maximum sensitivity of SkyHough since only an approximated frequency-evolution track is being used. Further developments on the pipeline include the re-analysis of interesting candidates using the exact frequency-evolution track and using more sensitive statistics in order to improve parameter estimation [131,132]. Also, the combination of LUTs can be easily accelerated using GPU parallelization [133].
The use of LUTs depends on a basic set of CW parameters, namely frequency and sky position, meaning it can be arbitrarily extended to look for different types of CW signals as long as the source's intrinsic and extrinsic frequency evolutions are uncoupled from the Earth's Doppler modulation. This includes transient-like CW signals, such as those produced by newborn NSs [134], or NSs in binary systems [133].

FrequencyHough
FrequencyHough is another pipeline based on the Hough transform [135,136]. As opposed to SkyHough, binarized spectrograms are mapped onto the frequency and spindown subspace ( f 0 , f 1 ). In this case, sky position is fixed for all the analyzed templates. This presents two main advantages with respect to SkyHough. First, it allows to increase frequency resolution independently of other parameters, resulting in a better parameter recovery. Second, all candidates analyzed at once are related to the same sky position. As a result, sky regions where templates tend to overlap with instrumental features can be dealt with in a simpler manner [137,138]. This different mapping, however, limits the generalization possibilities of FrequencyHough to linear relations between parameters.
This search is usually combined with different data formats: wide parameter-space searches are usually run using the so called Short Fourier Data Base (SFDB) [139], which includes time-domain cleaning of raw data to reduce the effect of transient noise. For searches at a specific parameter-space region at hand, such as supernova remnants, the Galactic center or a specific outlier from a wider search, Band Sampled Data (BSD) [140,141] is used to efficiently apply heterodyning filters, reducing the computing cost of analyzing such a local parameter-space region.
Several generalizations of these methods have been developed in order to look for other kinds of source. Power-law frequency evolution was covered using a generalization of FrequencyHough [142], allowing to search for binary neutron star merger remnants [143]. A study on the suitability of FrequencyHough to probe planitary-mass primordial black holes was F -statistic [65,68] Dual-harmonic Viterbi [152] F -statistic [68] Transient Viterbi [158] Norm. Fourier power [143,158] SOAP [154] Line-aware statistic -Viterbi 3.0 [153] B-statistic - presented in [17]. A method to conduct all-sky searches for CWs from the evaporation of boson clouds around spinning black holes was developed in [144], using the BSD framework [140,141]. A related method to place constraints on dark-photon dark matter interacting directly with GW interferometers was developed in [145] and applied in [146].

Viterbi searches
CW signals could contain stochastic contributions (spin wandering) whose behaviour would not be well represented by the standard model introduced in Eq. (1). This could affect objects in both isolated and binary systems [147,148], and should be taken into account in order to describe the underlying physics.
A simple approach is to describe the frequency-evolution model itself as a stochastic process, namely a Markov Chain (MC), and infer the most likely instantaneous frequency of a signal from the data. This is usually developed under the framework of Hidden Markov Models (HMM) [149][150][151][152][153], but an equivalent description can be done using Bayesian probability [154]. As discussed in [12], the search for ultralight boson-cloud evaporation around spinning black holes could also benefit from this kind of approach. A specific example of such a search, directed toward Cygnus X-1 using Advanced LIGO O2 data, was presented in [16].
Given a set of measurements at discrete times x = {x j , j = 0, . . . , N − 1}, we want to infer the instantaneous frequency of a signal in the data f = { f j , j = 0, . . . , N − 1}. This problem can be readily expressed as an inference one on f by means of Bayes theorem [155] where the sampling distribution of different measurements has been (conservatively) factored assuming logical independence. The prior probability distribution on the instantaneous frequency P( f ) is usually specified in terms of the initial frequency P( f 0 ) and the transition probabilities of the MC P( f j | f j−1 ) P( f 0 ) is generally taken as a uniform distribution over the searched frequency band; the choice of transition probabilities (thus prior probabilities) P( f j | f j−1 ) and sampling distribution P(x j | f j ) is dependent upon the method and source of interest, as summarized in Table  1. Regardless of the astrophysical scope, searches using a MC evolution model [Eq. (7)] obtain the most likely (maximum-posterior) frequency-evolution path f * using the Viterbi algorithm [159] f * = arg max f P( f |x) . (9) Two basic choices exist for transition probabilities [152], depending on whether the dominant frequency-drift time-scale is given by the source's spin wandering or secular spindown. The former allows transitions to any neighbouring frequency bin at each time step. The latter uses a so-called biased HMM: given a frequency at bin j, the following frequency bin must be equal or lower, f j+1 ≤ f j ; the specific number of bins is usually between two and three [65,156].
Since no template bank is involved, Viterbi searches are able to benefit from a wide variety of detection statistics at a small tuning cost. Specifically, binary modulations can be easily folded in using the J -statistic, which improves over the C-statistic [160] by combining frequency side-bands using complex weights. In order to track phase information, Viterbi 3.0 uses an efficient implementation of the B-statistic [71,72] first proposed in [113].
Finally, the SOAP pipeline [154] uses Bayesian spectral analysis to avoid relying on the F -statistic, looking for CW signals displaying a sinusoidal behaviour during short periods of time. To do so, the sampling distribution is taken proportional to the data's Schuster periodogram [161], adding an extra hypothesis to increase the robustness against strong monochromatic features in the data. This pipeline assumes the CW signal frequency is contained in a single frequency bin; as a result, the maximum recoverable spindown corresponds to that of a NS with a small ellipticity, in an analogous manner to the latest Falcon searches [46][47][48].

Machine learning
The use of machine-learning (ML) techniques in the search for CW signals follows one of two trends: either to classify and summarize the results of a search's main stage, as discussed in Sec. 3.2, or to substitute the search step itself, acting as a detection statistic. A recent review on ML applications to GW data analysis in general can be found in [162].
A first approach is to train a classifier directly over Fourier-transformed raw time-series data to distinguish the presence of a signal within background noise. The specific data format is dependant upon the signals being looked for: for CW signals, which last for long periods of time over narrow frequency bands. In [163,164], a convolutional neural network (CNN) is trained on the real and imaginary parts of short-time Fourier transforms. Slight variations on this proposal were employed to look for postmerger signals. These signals, compared to CW signals, last a shorter period of time over a wider frequency band [165]. In order to reflect the non-trivial time dependency of the signal, Ref. [166] took Fourier transforms using shorter time durations, finally using the data spectrogram as input for, again a CNN.
The second approach takes a less radical point of view. Instead of starting from raw data, machine learning algorithms are applied on the output of a search pipeline to construct a new detection statistic [167,168]. This approach could be beneficial, as the output of search pipelines typically enhances signal features across the output parameter space. The specific format with which the output data is better represented, however, remains a point of discussion..
Further applications of ML to post-process the output of a search [106,169,170] will be covered in Sec. 3.2.

Post-processing strategies
The main stage of a wide parameter-space search usually returns the loudest templates in terms of a specific detection statistic. A good portion of these templates are correlated, either because their corresponding time-frequency evolution tracks sweep over similar data or because of the presence of non-Gaussianities in the data, producing broad parameter-space artifacts [88,171,172]. The idea behind post-processing stages is to reduce the number of candidates into an affordable quantity to be followed up.
Complementarily, veto strategies can be applied in order to reduce further the number of candidates resulting from a search. A veto is a simple method to assess the consistency of a CW candidate with respect to a signal hypothesis. The outcome is usually a boolean answer; as opposed to a post-processing or follow-up stage, the primary objective is to quickly reject an inconsistent CW candidate at a low computing cost.
The following subsections summarize common post-processing strategies employed in contemporary CW searches. After a brief overview of conicidence and clustering steps, we review four families of vetoes. Other techniques employed in previous searches can be found in [137].

Coincidences
In a network of gravitational-wave detectors with a comparable level of sensitivity, a CW signal is expected to produce significant candidates in the analysis of every detector's data. Imposing a coincidence criterion, that is, focusing on common parameter-space regions highlighted in every single dataset, reduces the false-alarm probability of the search, as noise fluctuations are less likely to be coincident across different detectors than CW signals [102]. This approach obeys a robustness versus sensitivity trade-off, as combining data from different detectors into a single analysis increases the sensitivity of a search without increasing the required number of templates to be evaluated, keeping the computing budget under control [93].
This approach has been widely employed in FrequencyHough and SkyHough searches. In both cases, a parameter-space distance, based on an Euclidean ansatz, is used to identify closeby candidates in each detector's results. Other searches, such as Time-domain F -statistic, PowerFlux or Falcon, impose coincidence criteria based on the overlap of enhanced parameter-space regions across both datasets, rather than using a parameter-space distance.
In particular, as discussed in Sec. 2, Time-domain F -statistic is the most prominent user of this strategy [102]. As opposed to other semicoherent searches, it does not combine coherent segments into a semicoherent statistic, but looks for coincident candidates across multiple segments (including different detectors), reducing the overall false-alarm probability of the search.

Parameter-space clustering
Another option is to group together nearby templates according to some notion of distance. This process, usually referred to as clustering in the data analysis literature, has been extensively employed by PowerFlux, Einstein@Home, FrequencyHough, and SkyHough using different implementations, both in terms of clustering strategy and parameter-space distance [40][41][42]45,49,50,52]. Clustering is typically implemented using an unsupervised approach [130,173]: the parameter-space and the clustering algorithm itself act as prior information to construct meaningful groupings of the resulting template bank. Supervised ML approaches [106,169,170], aimed at identifying specific parameter-space structures, have also been proposed.
Unsupervised approaches try to unveil structure from a data set. Prior information is encoded both in terms of the distance used to compare nearby candidates and the linkage criteria. The resulting clusters are sieved through selection criteria which can take into account parameters like the maximum significance in the cluster or its number of elements. We focus our exposition on the choice of parameter-space distance; a complete description of the clustering algorithms themselves is given in the corresponding references. Initial implementations, used by FrequencyHough and SkyHough [41,42,50] assumed a Euclidean parameter-space distance on the CW signal parameters λ = { f 0 , f 1 , . . . , n, .
pairing candidates within a certain distance threshold to form the final clusters. The exact parameters λ (i) and parameter-space resolutions δλ (i) are search-dependent. A more informative approach, still based on a Euclidean ansatz, was proposed in [173]. In this case, clusters are classified according to topographic parameters by projecting detection statistics over planes, namely over the frequency-spindown plane and the ecliptic plane. Instead of using the basic CW parameters, distance was computed after projecting the sky position of the candidate onto the ecliptic plane; thus allowing a greater variance around the ecliptic, where sky localization tends to become more uncertain [80,107].
These distance measures are effective for local analyses, but quickly become unrealiable whenever more involved parameter-space structures such as correlated parameters with periodic boundaries come into play, as is the case for signals from sources in binary systems [174]. A parameter-space distance for a generic, quasi-monochromatic CW signal was proposed in [130] using the instantaneous detector-frame frequency associated to a CW template f (t; λ). Concretely, the distance between two templates λ and λ * can be defined as the average mismatch between their corresponding detector-frame frequency tracks throughout an observing run This prescription is consistent with the F -statistic's parameter-space correlations and can be simplified into a discrete sum for a faster implementation [52,53,130]. Candidate post-processing, and clustering in specific, is also a suitable step in which machine learning strategies are able to deliver an improvement of sensitivity. As opposed to raw data, on which CW signals are typically a subdominant contribution, the structures produced by different features in the data on the parameter space of a search are suitable to be classified using a supervised approach, as long as a clear classification of the features at hand is available. To date, two approaches have been proposed.
The first one [106], framed within the Time-domain F -statistic pipeline, uses a convolutional network to classify different representations of the main search's output into three possible classes, namely noisy bands containing Gaussian-like noise, narrow spectral artifacts, and CW signals. A similar (albeit more manual) approach was reported in [40].
The second approach is developed as an alternative to the clustering algorithm employed in Einstein@Home searches [169,170]. In this case, a neural network is trained to recognize signal-induced patterns on a certain projection of the parameter space in order to identify typical structures associated to CW signals. The training set is produced by manually identifying software-injected signals using an image editing tool.

Detector-consistency vetoes
The first family of vetoes tests the consistency of a CW candidate across the network of GW detectors. The basic implementation is formulated as follows: given a CW candidate with parameters p and N detectors, a detection statistic S(p) is computed using data from each detector alone {S 1 (p), . . . , S N (p)} and combining the datasets from all detectors at once S M (p). Then, a function of single-detector statistics F{S 1 (p), . . . , S N (p)} is compared to the multi-detector statistics in order to decide whether the candidate behaves consistently with a CW signal or not. The decision boundary, usually expressed in terms of the difference F{S 1 (p), . . . , S N (p)} − S M (p), can be calibrated by means of a software-injection campaign. As discussed in [86][87][88][89], the B S/GL and B S/GLtL detection statistics are a Bayesian approach to this implementation of the detector-consistency veto. In their implementation [40,49,175], however, the information is processed during search-time, discarding inconsistent candidates at an earlier stage and potentially improving the detection of marginally significant signals.
An example of detector-consistency vetoes can be found in the SkyHough contribution to [41], where F was taken to be the expected multi-detector statistic computed from the single-detector statistics including the varying noise floors but ignoring the antenna pattern modulations. Another example is found in the Weave search [51], where F is simply the maximum detection statistic over all the involved detectors. A slight variation was employed in the BinarySkyHough analysis of early O3 data [52,53], where the detection statistic of one of the LIGO detectors was compared against the other one. This was motivated due to the asymmetric behaviour of said detectors during the third observing run, with H1 more affected by noise disturbances than L1.

χ 2 vetoes
The second family of vetoes considers the behaviour of a putative CW signal within a particular dataset, namely, whether the detection statistic accumulates throughout the observing run in a way which is more consistent with an instrumental artifact than an astrophysical signal. This is the idea behind the χ 2 veto, initially proposed in [176] for CBC signals and later implemented for CW signals in [177]. In this case, the dataset is partitioned into p segments over which a detection statistic is computed {S 1 , . . . , S p }. These segment-wise statistics are then compared to the expected segment-wise statistics under the presence of a signal in Gaussian noise, usually characterized by a mean and standard deviation µ and σ, and combined into a chi-squared discriminant Under the assumption of Gaussian noise, this discriminant follows a χ 2 distribution with p − 1 degrees of freedom; real-data applications, however, must calibrate this test using a suitable injection campaign [41]. In [52,53], extreme deviations of the segment-wise detection statistic were used to identify stretches of data in which the CW template showed a high degree of overlap with an instrumental feature; this approach is equivalent to using a χ 2 discriminant in the limit of a very strong sample S j − µ j σ j .

Vetoing narrow spectral features
The third family of vetoes relies on detector characterization to identify frequency bands in which an instrumental feature is present. CW searches integrate long periods of time looking for quasi-monochromatic signals concentrated around a fraction of a Hertz. Quite often, those narrow bands are populated by narrow spectral features (lines) due to instrumental or environmental disturbances (defective power supplies, blinking LEDs, wind blowing, local fauna interacting with the detector. . . ) which, under very general conditions, are able to mimic the effect produced by a CW signal in the detector, usually producing a high number of candidates in a search [171,178,179]. Catalogs listing narrow spectral features and their cause (if known) for the latest runs of the advanced detectors are publicly available [180][181][182][183][184][185][186]. This kind of features are generally not a problem for CBC searches, as those signals sweep wide frequency bands in a relatively short time duration, although noise subtraction techniques are applied in severe cases [187,188].
A common approach, usually referred to as the (known) line veto (see e.g. [41,42,45,52]), is to check whether the frequency evolution of a CW candidate overlaps with any frequency band containing such instrumental features, in which case the candidate is discarded. This requires a high degree of manual intervention, as line catalogs must be created and properly understood, and incurs the risk of removing a genuine signal candidate due to an unfortunate line crossing at a potentially insignificant period of the run.
The DM-off veto was proposed in [189] as a hypothesis-test version of the line veto: it compares the significance of a CW candidate, which includes Doppler modulations due to the Earth's movement around the SSB, versus the significance obtained after analyzing its surrounding frequency band using an unmodulated template bank (i.e. without Doppler Modulation, hence DM-off). This veto was applied with great success in [40]. As happens with the detector-consistency veto and the B S/GL statistic [88], the DM-off veto can be refactored as a detection statistic in a Bayesian framework to construct a proper line-robust statistic: instead of testing against a single-detector artifact, the hypothetical B S/GMU statistic would test against an ensemble of monochromatic and unmodulated signals in any number of detectors [178].
Alternatively, as performed in [40,49,190], frequency bins in the data containing lines can be replaced by Gaussian noise drawn from the distribution of neighboring bins to suppress the presence of candidates, preventing any candidate of instrumental origin from polluting the search results. This process is generally performed on SFT data before starting a search.
Short-duration loud instrumental glitches also present a problem to CW searches, as they tend to degrade the noise floor across a wide frequency band. This sort of artifacts, however, are typically dealt with before starting a search using a cleaning procedure such as gating [139,179,191,192].

Null-hypothesis vetoes
The fourth family of vetoes are essentially a reformulation of the standard null-hypothesis test, in which a CW candidate is deemed as uninteresting if it is consistent enough with respect to the background noise distribution. A simple proposal, usually referred to as offsourcing [64,193], is to evaluate a CW candidate on nearly independent noise realizations by shifting its sky position away. This is based on the fact that detector artifacts tend to imprint wider parameter-space regions with significant templates than CW signals. Off-sourced time-frequency tracks are able to break signal-induced correlations while still being affected by instrumental disturbances. The resulting distribution is a proxy of the noise hypothesis' sampling distribution and can be used to construct significance arguments about the CW candidate of interest.
This veto was adapted by [51] to evaluate the final surviving candidate of a search. In their use case, however, they considered the distribution of the loudest candidate from a CW search, that is, considering the number of trials performed by evaluating a template bank on a data stream. Such a distribution has been previously studied in the CW literature [194,195], but it was not until recently that a method applicable to a generic CW search based on extreme value theory was proposed [172].
Steps towards a fully Bayesian treatment of loudest-candidate null-hypothesis vetoes were taken in [38], which proposed a Bayes factor to evaluate the loudest candidate of a CW search, B * S/G , whose noise-hypothesis component was constructed fitting a Gumbel distribution to the loudest outliers of off-sourced template banks. This approach was used to develop a complete hierarchical follow-up framework, discussed in Sec. 4.2.

Follow-up
Follow-up stages are contextualized within a hierarchical search, as discussed in the introduction. They improve the parameter estimation of a CW candidate by imposing tighter constrains on its expected behaviour. This leads to the factual use of simple follow-up stages as signal-consistency veto strategies. Base search stages construct less-sensitive detection statistics by effectually using less-constraining CW signal models. For example, a semicoherent search, in which phase information is contained in discrete, non-overlapping segments, is insensitive to arbitrary phase jumps between coherence segments. In this sense, the sensitivity loss with respect to a fully-coherent search is due to the increased trials factor of this looser family of signal models [43].
There are two ways in which follow-ups may be performed. Following the notation established in [33], we refer to them as fresh data mode (FDM) and recycling data mode (RDM). As their names suggest, FDM looks for a CW candidate in a new dataset containing brand new information; RDM, on the other hand, re-analyzes the same dataset using a different method. The typical example of FDM is to look for a CW candidate obtained in a certain observing run using data from subsequent observing runs [40,54,62]. Examples of RMD include, for example, multi-stage semicoherent or loosely-coherent searches aiming towards a fully-coherent search in a restricted parameter space region [36][37][38][39].
Most CW searches are conducted focusing on a single dataset, usually the latest available observing run of the advanced detectors. Hence, follow-up strategies operate on RDM (even though sometimes FDM assumptions are used for simplicity, as they turn out more conservative and simpler to implement [195,196]). Nevertheless, this sort of strategies could be detrimental towards detecting a CW signal, as standard RDM operations leading to effectively longer coherence time may lose candidates affected by some form of unmodelled behaviour, such as glitches [197,198] or accretion-induced spin wandering [148]. Strict FDM must be used in order for a search to follow up a CW candidate with a consistent signal model. Looking into previous observing runs, on the other hand, runs into a sensitivity problem, as marginal candidates may end up completely lost due to the lower quality of the detectors. Examples of searches in which a FDM follow-up looked into a posterior observing run include [40,55].

Single-stage follow-up
The simplest follow-up strategy calibrates a threshold on a different (more sensitive) detection statistic according to some prescription and compares it to the score returned by reanalyzing the CW candidate. As opposed to a veto, this approach points towards evaluating the consistency (or discrepancy) of the CW candidate with respect to a certain population of signals.
The standard FrequencyHough follow-up as employed in [41,45] belongs to this category: baseline Fourier transform length is increased, imposing tighter constraints on the signal model and discarding short-duration candidates. Later searches also make use of the BSD framework [140,141]. Similar strategies, in this case using the fully-coherent F -statistic, were proposed in [199,200].
In order to overcome the curse of dimensionality, BinarySkyHough searches [50,52,53] employ an MCMC-based follow-up implemented in PyFstat [37,38,201]. The multi-detector F -statistic [21,70] is used to allow for arbitrarily long coherence times. In this sense, following up a CW candidate is equivalent to sampling the posterior probability distribution of the phase-evolution parameters λ given a data stream [37] P(λ|x) ∝ e F (λ;x) · P(λ) , where the prior probability distribution P(λ) represents the parameter-space region of interest identified by a search. The result of this grid-less approach is the F -statistic evaluated at the loudest candidate of the parameter space at a negligible mismatch. This approach generalizes that of [199,200], which specializes in the single-stage fully-coherent follow-up of CW candidates. As discussed in Sec. 4.2, the use of MCMC methods simplifies the setup of a generic multi-stage follow-up, as no calibration of parameter-space grids are required. The onus in this case is on the search pipeline to deliver a small-enough prior support for the MCMC to converge. This can always be achieved by starting from a shorter coherence time [37], at the expense of increasing the number of follow-up stages.
A similar strategy is used by the Time-domain F -statistic follow-up, focusing on the optimization aspect of the procedure. In this case, a max-finding algorithm such as [202] is employed to travel around the parameter space. The result, as with PyFstat, are the most favored signal parameters, corresponding to the ones reporting the loudest F -statistic value.

Multi-stage follow-up
As discussed during Sec. 1, multi-stage follow-ups are the natural continuation to a wide parameter-space search after identifying interesting parameter-space regions. Given a data stream, each subsequent follow-up stage operates in RDM, gradually increasing the coherence time with respect to previous stages and, as a consequence, imposing a tighter version of the selected signal model. The effect on a CW candidate is twofold: first, non-astrophysical CW-like artifacts tend to get rejected as coherence time increases (though, as earlier discussed, this could have detrimental effects on more complex CW signals too); second, increasing coherence time results in a refinement in parameter-space resolution, improving the parameter estimation of the candidate at hand. This last effect must be considered carefully, as it also implies an increased number of templates to analyze, quickly becoming unaffordable if parameter-space regions are not gradually narrowed down.
We start by discussing the follow-up strategy introduced in [36], which was applied as a generic follow-up to multiple CW searches with minor modifications [40][41][42]49]. This example is paradigmatic in the sense that it fully exposes the two main challenges of a multistage setup. A similar approach, albeit at much smaller scale, was employed to follow-up SkyHough results in [45]. First, one must set up a proper set of parameter-space grids such that CW signals are not lost due to a bad parameter-space covering. If an analytical model of the follow-up method at hand is not available, as is the case in [36], one must resort to an extensive software injection campaign to construct a suitable setup. Second, a criterion must be set up to select/reject CW candidates after performing different stages. If the detection statistic's behaviour across different stages is well understood (see e.g. [196]), an analytical criterion can be derived from first principles; otherwise, a software injection campaign must be used to calibrate a rejection criterion.
Latest developments on the follow-up of CW candidates [37,38] are able to simplify the setup for the F -statistic, although further work is required for its application to generic detection statistics. As discussed in Sec. 4.1, the use of MCMC samplers simplifies the setup due to their lack of grids: if a CW signal is within the prior support, parameter-space samplers can get arbitrarily close to the injection parameters given enough time to wander around the prior volume. This argument can be posed quantitatively in terms of the so-called coherence time ladder [37], which make use of the parameter-space metric to increase the parameter-space resolution in a controlled manner. Since the follow-up is typically a local analysis, rough estimates of the number of templates using typical parameter-space resolutions are usually a valid approximation [198].
Proper comparison of detection statistics from different stages in a Bayesian framework is currently restricted to the F -statistic [196], as the sampling distribution under a signal hypothesis must be known given a (squared) signal-to-noise ratio ρ 2 . This result was used in [38] to propose B * S/N , a (meta) Bayes factor (as the F -statistic itself is a Bayes factor) evaluating the result of a multi-stage follow-up, pushing forward the development of a fully-Bayesian follow-up of CW candidates. In this case, the probability under the noise hypothesis was derived from a combination of off-sourcing and extreme value theory [172]. The use of detection statistics (Bayes factors) as data proxies to construct Bayesian arguments is also discussed in [195,203].
An alternative family of follow-up methods were developed under the name of loose coherence [39,113,114], already introduced in Sec. 2.2.1. In this case, instead of following the ad hoc recipe of increasing coherence time until a fully-coherent search is achieved, phase information across neighbouring time segments is gradually correlated in a controlled manner by combining (complex) Fourier amplitudes. As opposed to semicoherent methods, which allow for arbitrary phase jumps at the border of a segment, loosely coherent methods allow for phase shifts within pre-specified ranges, depending on the required robustness of the method. Under this framework, semicoherent methods are rediscovered imposing delta-correlation between phases at consecutive time segments. This follow-up approach is fully integrated within the PowerFlux and Falcon searches [41][42][43][44]46,47,55]

Upper bounds on h 0
Since no CW detection has been reported to date, the main data product of CW searches are bounds on the nominal gravitational-wave amplitude h 0 produced by a population of sources consistent with the target signal model. As discussed in Sec. 1, astrophysical information can be extracted from this quantity by taking different assumptions, such as the maximum allowed ellipticity from a galactic neutron star at a certain distance from the detector.
Two basic approaches are pursued at production level to derive said upper bounds, depending on whether the aim is for a strict frequentist upper limit or population-based sensitivity estimations. The calibration and establishment of upper bounds of any kind usually involves an extensive software injection campaign in real data with a non-negligible computing cost; due to this, a common approach lies in between both extrema, quoting a proper estimation at a definite set of representative frequency bands and interpolating the results across the rest of the spectrum.
Population-based sensitivity estimations are based on estimating the false-dismissal probability of a search given a certain setup (be it a threshold at a fixed false-alarm probability or a more intricate procedure). The p% detection probability amplitude h p% 0 corresponds then to the amplitude h 0 associated to a false dismissal of (1 − p)% after properly marginalizing with respect to other amplitude parameters using a set of priors reflecting the studied source distribution [194]. In practise [40][41][42]45,49,50,52], detection probabilities are usually estimated numerically by means of an injection recovery campaign.
Strict frequentist upper limits, on the other hand, return a conservative estimate of the upper bound, in the sense that false-dismissal probability is at most (1 − p)%. An example of this kind is the universal statistic procedure [204], which is able to construct strict frequentist upper limits regardless of the underlying noise distribution. This procedure has been extensively combined with the PowerFlux and Falcon pipelines to efficiently produce robust upper limits under different GW polarization assumptions [41][42][43][44]46,47,55].
To date, all wide parameter-space searches have made use of one of these two upper bounds to report on their results. These upper bounds, however, describe the probability of detecting a signal given an ensemble of equivalent noise realizations, rather than the range where a signal could be found given the data stream at hand. Work towards reporting the latter, Bayesian upper bounds, for wide parameter-space searches was developed in [195,203].

Summary
We reviewed the methods employed by current wide parameter-space searches for continuous gravitational waves from unknown sources conducted on advanced-detector data. The most widespread approach consists of a hierarchical setup in which parameter-space regions are analyzed using more sensitive (and consequently more expensive) methods as they are gradually narrowed-down. Detecting a CW signal requires both an instrumental and computational effort to confidently unveil such a weak signal using the current generation of gravitational-wave detectors. The use of multiple methods taking different tradeoffs in sensitivity and robustness against instrumental artifacts provides an ideal environment to pursue new strategies towards CW detection and parameter estimation. For the sake of completeness, Table 2