Einstein@Home all-sky search for periodic gravitational waves in LIGO S5 data

This paper presents results of an all-sky searches for periodic gravitational waves in the frequency range [50, 1190] Hz and with frequency derivative ranges of [-2 x 10^-9, 1.1 x 10^-10] Hz/s for the fifth LIGO science run (S5). The novelty of the search lies in the use of a non-coherent technique based on the Hough-transform to combine the information from coherent searches on timescales of about one day. Because these searches are very computationally intensive, they have been deployed on the Einstein@Home distributed computing project infrastructure. The search presented here is about a factor 3 more sensitive than the previous Einstein@Home search in early S5 LIGO data. The post-processing has left us with eight surviving candidates. We show that deeper follow-up studies rule each of them out. Hence, since no statistically significant gravitational wave signals have been detected, we report upper limits on the intrinsic gravitational wave amplitude h0. For example, in the 0.5 Hz-wide band at 152.5 Hz, we can exclude the presence of signals with h0 greater than 7.6 x 10^-25 with a 90% confidence level.

This paper presents results of an all-sky search for periodic gravitational waves in the frequency range [50, 1 190] Hz and with frequency derivative range of ∼ [−20, 1.1] × 10 −10 Hz s −1 for the fifth LIGO science run (S5). The search uses a non-coherent Hough-transform method to combine the information from coherent searches on timescales of about one day. Because these searches are very computationally intensive, they have been carried out with the Einstein@Home volunteer distributed computing project. Post-processing identifies eight candidate signals; deeper follow-up studies rule them out. Hence, since no gravitational wave signals have been found, we report upper limits on the intrinsic gravitational wave strain amplitude h0. For example, in the 0.5 Hz-wide band at 152.5 Hz, we can exclude the presence of signals with h0 greater than 7.6 × 10 −25 at a 90% confidence level. This search is about a factor 3 more sensitive than the previous Einstein@Home search of early S5 LIGO data.

I. INTRODUCTION
A promising class of sources for detectable gravitational wave signals are rapidly rotating neutron stars with non-axisymmetric deformations [1][2][3][4][5]. Such objects are expected to emit long-lived continuous-wave (CW) signals. In the rest frame of the neutron star, these waves have a constant amplitude and are quasi-monochromatic with a slowly decreasing intrinsic frequency. They are received at Earth-based detectors with a Doppler modulation due to the relative motion between the source and the detector. Consequently the observed phase evolution depends on the intrinsic signal frequency, the first frequency time-derivative (also called spindown), and the source sky position; these parameters shall collectively be called the phase evolution parameters. While using higher order frequency derivatives could potentially improve the astrophysical detection efficiency in a part of the parameter space (see Sec. III), we shall not consider them in this paper for computational reasons. Finally, the received signal has a time-dependent amplitude modulation due to the (time-dependent) relative geometry of the wave and the detector.
The previous two decades have seen the construction and operation of several kilometer-scale laser interferometric gravitational wave detectors [6][7][8][9][10][11]. The detectors and the data analysis tools have steadily improved over this period. These have made it possible to search for various gravitational wave signals with ever-improving sensitivity. In this paper we focus on data from the fifth science run (S5) of the LIGO (Laser Interferometer Gravitational wave Observatory) detectors, collected between the GPS times of 815 155 213 s (Fri Nov 04 16:00:00 UTC 2005) and 875 145 614 s (Sun Sep 30 00:00:00 UTC 2007). The LIGO network [6] consists of three detectors: one at Livingston, Louisiana, USA, with an arm length of 4 km (L) and two in the same vacuum envelope at Hanford, Washington, USA, with arm lengths of 4 km (H) and 2 km, respectively. Only data from H and L detectors are used in this paper. The Virgo and GEO 600 detectors also collected data during the same time interval but are not used in this analysis, which is optimized for two detectors with similar sensitivities.
A coherent strategy for extracting faint CW signals buried in noisy data using standard maximum-likelihood techniques in the presence of "nuisance parameters" was derived in [12]. The resulting detection statistic is the so-called F-statistic, which has since been generalized to the case of multiple detectors [13,14]. The F-statistic has also been shown to arise as a special case in a more general Bayesian framework [15]. Using the F-statistic means that we need to search explicitly only over the phase evolution parameters.
Coherent wide-parameter-space searches utilizing the F-statistic have been carried out since the second LIGO science run [16,17]. The amplitude sensitivity of this type of search improves as the square root of the time baseline. However, the template bank spacing decreases dramatically with the time baseline, and the computational requirements of the search increase rapidly. Even with a coherent time baseline of just a few days, a wide-frequency-band all-sky search is computationally extremely challenging. It becomes completely unfeasible if one considers instead time baselines on the order of months.
As is often the case with computationally bound problems, hierarchical approaches have been proposed [18][19][20]. In these strategies, the entire data set is split into shorter segments. Each segment is analyzed coherently, and afterwards the information from the different segments is combined incoherently (which means that the phase information is lost). The amplitude sensitivity grows at best with the fourth root of the number of segments. Such methods have been used in previous wide-parameter-space searches published by the LIGO and Virgo collaborations [21][22][23][24][25][26].
A subset of these searches [21,22,25,26] used segments sufficiently short (1 800 s) that the signal remains within a single Fourier frequency bin in each segment. In these cases, a simple Fourier transform suffices for each segment. Three different variants of such methods have been developed that combine the results from the different short segments incoherently: the "stack-slide", the "Hough-transform" and the "PowerFlux" schemes. The stack-slide procedure [18,20] averages the normalized power from the Fourier transform of 30-minute segments (Short time baseline Fourier Transform, SFT for short) of the calibrated detector strain data. The Pow-erFlux scheme [22,25] can be seen as a variation of the stack-slide method, where the power is weighted before summing. The weights are chosen according to the detector noise and antenna pattern to maximize the signal-tonoise ratio (SNR). The Hough-transform method [19,27] sums weighted binary counts, depending upon whether the normalized power in an SFT bin exceeds a certain threshold.
As the segment duration is increased, it becomes necessary to account for signal modulations within each segment by computing the F-statistic over a grid in the space of phase evolution parameters. This results in a significant increase in the computational requirements of the search. The distributed volunteer computing project Ein-stein@Home [28] has been created to address this need. Two previous papers [23,24] report on results of such CW searches from the fourth LIGO science run and from the first two months of S5, respectively. The method used was based on the computation of the coherent F-statistic on data segments from either the H or L detectors separately, and only parameter space points with values of 2F larger than 25 were returned back to the Einstein@Home server for further inspection. The threshold value of 25 limited the ultimate sensitivity of that search: if a signal was not loud enough to surpass that threshold on at least some of the segments it would not be detected. The threshold value was set by bandwidth constraints on the size of the results file uploaded back to the server by the host, i.e. on the maximum number of significant points that could be returned. These results were subsequently combined by a coincidence scheme, performed offline in the post-processing phase. In contrast, in the search presented here, the combination of the results from the coherent searches takes place directly on the host machines using a Hough-transform scheme. This makes it possible to use a much lower threshold on 2F, equal to 5.2, that defines the parameter space points to be passed on to the Hough-transform. Moreover, in this search, data from the H and L detectors are coherently combined [13,14]. Finally, more data was searched in this analysis compared to any previous Einstein@Home search. The Ein-stein@Home runs presented here refer to searches based on the first (S5R3) and second year of S5 LIGO data. This latter search was run on Einstein@Home in two separate steps, called S5R5 and S5R6. Since the S5R6 run used the same data as S5R5, but extended the search region above 1 kHz, in this paper we simply refer to these two runs as S5R5.
The paper is structured as follows. In Sec. II, we briefly review the Hough-transform method. Section III describes the Einstein@Home distributed search used to analyze the data set. Section IV gives a detailed description of the S5R5 post-processing, which is based on the pioneer S5R3 post-processing (described in Appendix B). Upper limit computations from the more-sensitive S5R5 data are provided in Sec. V. The study of some hardwareinjected signals is presented in Sec. VI. In Sec. VII we make some concluding remarks.

II. THE DATA ANALYSIS METHOD
A. The waveform model Let us begin by briefly describing the standard signal model for CW signals. In the rest frame of the neutron star, the gravitational wave signal is elliptically polarized with constant amplitudes A +,× for the two polarizations h +,× (t). Thus, we can find a frame such that The two amplitudes are related to an overall amplitude h 0 and the inclination angle ι between the line of sight to the neutron star and its rotation axis The value of h 0 is model dependent. For emission due to non-axisymmetric distortions, the amplitude h 0 depends on the ellipticity ε of the star defined as Here I zz is the principal moment of inertia of the star, and I xx and I yy are the moments of inertia about the other axes. The amplitude is given by where f is the frequency of the emitted GW signal (which is also twice the rotational frequency of the star), G is Newton's constant, c is the speed of light and d is the distance to the star. The distribution of ε for neutron stars is uncertain and model dependent since the breaking strain for a neutron star crust is highly uncertain (see e.g. [2,[29][30][31] for further discussion). Energy loss from the emission of gravitational and/or electromagnetic waves, as well as possible local acceleration of the source, causes the signal frequency arriving at the solar system to evolve. To first order, it can be expressed asf where τ is the arrival time of a wavefront at the solar system barycenter (SSB), f 0 is the frequency at a fiducial reference time τ 0 , andḟ denotes the first time derivative of the frequency. The astrophysical implications of ignoring higher order derivatives in this Taylor expansion will be discussed later. The phase of the signal, φ(t), follows directly from the frequency evolution with an initial phase φ 0 at the reference time.
As the detector on the Earth moves relative to the SSB, the arrival time of a wavefront at the detector, t, differs from the SSB time τ : Here r(t) is the position vector of the detector in the SSB frame, and n is the unit vector pointing to the neutron star; ∆ E and ∆ S are respectively the relativistic Einstein and Shapiro time delays [32]. In standard equatorial coordinates with right ascension α and declination δ, the components of the unit vector n are given by (cos α cos δ, sin α cos δ, sin δ).
Ignoring the relativistic corrections, the instantaneous frequency f (t) of a CW signal, as observed at time t by a detector on Earth, is described by the well-known Doppler shift equation: where v(τ ) is the detector velocity with respect to the SSB frame; v(τ ) is the sum of two components, from the yearly Earth motion around the Sun ( v y ) and from the rotation of Earth around its axis ( v d ).
Finally, the received signal at the detector is where F +,× are the detector beam pattern functions which depend on the sky position n and the relative polarization angle ψ of the wave-frame [12,33]. There are thus altogether eight signal parameters, which include the four phase evolution parameters (f 0 ,ḟ , α, δ), and four other parameters (h 0 , ι, ψ, φ 0 ).

B. The Hough-transform algorithm
For completeness, we summarize the Hough detection statistic in this section. Further details of the method are given in [19] and previous searches with this method applied to short coherent times are reported in [21,22].
The Hough-transform is a well-known technique used mainly in digital image processing for robust extraction of patterns. Such a procedure is employed here for identifying points in the time-frequency plane that match the pattern expected from a signal. The time-frequency data in our case is the F-statistic computed as a function of signal frequency, for each of the data segments of duration T seg , over a grid of points in the space of (α, δ,ḟ ). The grid in (α, δ,ḟ ) space used for this F-statistic computation is called the coarse grid because its resolution is determined by the coherent time baseline T seg , and makes no reference to the full observation time or the number of segments N seg . The result of this computation is thus a collection of F-statistic values F i α,δ,ḟ (k) where the integers k and i label a frequency bin (with spacing δf as defined below) and a data segment respectively.
The frequency and frequency derivative spacings for the coarse grid are based on choosing the maximum allowed fractional loss in the F-statistic when the signal and template points are slightly mismatched. This leads naturally to the notion of a metric in parameter space [34,35] and this has been studied for the CW case in [14]. The grid spacings in f,ḟ are given respectively by [17,36] δf = √ 12m πT seg (9) and where m represents the nominal single dimension mismatch value. For all the Einstein@Home runs discussed here m has been taken equal to 0.3. The span T seg of each segment has been set equal to 25 hours for all the runs. The frequency resolution, given by Eq. (9), is δf ∼ 6.7 µHz for all the Einstein@Home runs described in this paper. As we shall see shortly, for technical reasons it turned out to be necessary to use a finer spacing forḟ than given by Eq. (10). In combining these N seg different F-statistic vectors, it is necessary to use a finer grid in (α, δ,ḟ ) centered around each coarse grid point. Our implementation of the Hough-transform algorithm assumes that the fine grid is a Cartesian product of a rectangular sky-grid and one dimensional grids in f andḟ . Moreover the fine skygrid is assumed to be aligned with the α and δ directions in the sky. In order to completely cover the sky with the different fine sky-grids, it is thus simplest to choose a rectangular coarse sky-grid aligned with the (α, δ) directions. We choose a coarse grid such that the spacing in δ is a constant and the spacing in α is proportional to (cos δ) −1 . This ensures that each cell of the coarse grid covers a fixed solid angle.
Since the coherent integration time is very close to a sidereal day, it is reasonable to assume that the Hough sky-patch size dθ (which in our case is the same as the coarse sky-grid resolution) is determined by v d , the Earth's rotation speed at the equator. At frequencyf we have [19]: In practice, this estimate was verified by Monte-Carlo studies with signal injections, and we used where the factor R = √ 3 increases the size of the skypatch dθ for the S5R5 run. The Monte-Carlo studies showed that the grid spacing for S5R5 corresponded approximately to a mismatch of m ≈ 0.3 so that, for any other value of m, the spacing would be approximately m/0.3Rdθ.
We also need to set the resolution of the refined skygrid used by the Hough algorithm. Since the full observation time is of the order of a year, the relevant scale here is set by the speed v y of Earth as it orbits the Sun. Following [19], the resolution for the fine sky-grid at a frequencyf is given by The parameter ℘ scales the resolution compared to the conservative estimate c δf /(f v y ) and in practice, again based on Monte-Carlo studies, we used ℘ = 0.5. One can see that the increase in the number of sky position points from the coarse F-statistic grid to the fine Hough grid is Finally, let us turn to the coarse and fine grids foṙ f . Ideally, we should refine the coarseḟ grid spacing of Eq. (10) by a factor N seg [19]. However, in our implementation of the Hough transform, using this refinement turns out to increase the maximum memory footprint of the different searches and would make it unsuitable for Einstein@Home. As a compromise, it was decided not to use any refinement inḟ and instead to use a finer resolution for the coarse grid. Based on Monte-Carlo analyses, an acceptable value for theḟ spacing turns out to be This value was used for S5R5. However, for S5R3, this was incorrectly set to leading to a corresponding loss in sensitivity for S5R3. The flow chart of the search algorithm used for this search is depicted in Fig. 1. The input data set is composed of 30-minute baseline SFTs. This set is partitioned in subsets such that no more than 25 hours of data are spanned by each segment and such that there is overall (including data from both detectors) at least 40 hours of data in each segment. Let T obs be the observation time spanned by the N seg segments constructed in this way 1 . The multi-detector F-statistic is computed for each segment at each point of the search parameter space (f ,ḟ , α, δ). The next step consists of selecting parameter space points for which the F-statistic is above the fixed threshold. For every set of (f ,ḟ , α, δ), we assign a value n i = 1 or 0 in the i th segment depending on whether the corresponding F-statistic is above the threshold 5.2 or not; this threshold turns out to be optimal [19]. The values n i (f ) are called a "peakgram", which is the input to the Hough-transform.
The final statistic used by the Hough search is a weighted sum of binary counts n i , giving the so-called Hough number count n c , expressed by [22] The weight w i for a frequency f and for a particular skylocation is determined from the average antenna response and average detector noise over the duration of the i th segment. Since the input data in each segment consists of SFTs, we perform the averaging over each SFT. Let N i be the number of SFTs in the i th segment. Let S i,γ h be the single-sided power spectral density (PSD) of the γ th SFT in the i th segment, and averaged over a narrow frequency band containing the search frequency. The w i are given by where F +(i,γ) and F ×(i,γ) are the detector antenna pattern functions for the γ th SFT in the i th segment, and High level schematic of the pipeline used in the searches. For each data segment, the multi-detector Fstatistic is computed and frequency bins are selected setting a threshold on the F-statistic. Such selected frequency bins are then used to create the Hough map. The output is a set of candidates in the parameter space. The color-filled boxes indicate the steps performed on the volunteer computers. The acronym "WU" is defined in Sec. III A and refers to the independent computing tasks into which we partition the computational work of the search.
It is easy to see that, in the hypothetical case when the data is exactly stationary, so all the S i,γ h are identically equal to each other, then w i,γ = 1. More realistically, the w i,γ are approximately unity for stationary data and the use of the harmonic mean in Eq. (19) ensures that the w i,γ do not deviate too far from unity in the presence of non-stationary noise.
The weight normalization is which ensures that the Hough number count n c lies within the range [0, N seg ]. In [22] it is shown that the weights w i , first derived in [37], maximize the sensitivity, averaged over the orientation of the source (see Appendix A for a further discussion of the weights and some technical problems that were encountered in the search). From the Hough number count n c we define the significance (or critical ratio), CR: that measures the significance of the measured n c value as the deviation from the expected valuen c = N seg p in absence of any signal, in units of the expected noise fluctuations σ; p is the probability that a parameter space pixel is selected in the absence of a signal. In case of unity weighting, the standard deviation is simply that of the binomial distribution: σ = N seg p(1 − p). When the weights are used, the standard deviation is given by where || w|| 2 = Nseg−1 i=0 w 2 i [22]. The CR is the detection statistic returned by the hierarchical searches presented here.

III. THE EINSTEIN@HOME DISTRIBUTED SEARCH
The Einstein@Home project is built upon the BOINC (Berkeley Open Infrastructure for Network Computing) architecture [38-40], a system that exploits the idle time on volunteer computers to solve scientific problems that require large amounts of computer power. During the S5R5 run, Einstein@Home had approximately 225 000 registered volunteers and approximately 750 000 registered host machines that contributed a total of approximately 25 000 CPU (Central Processing Unit) years.

A. Details of the S5R5 search
The computational load is partitioned in independent computing tasks, called "workunits" (WUs), each of which is analyzed by a volunteer machine. In particular, 7 369 434 and 10 945 955 WUs have been generated for the S5R3 and S5R5 runs, respectively. The The analyzed data consists of 5 550 and 5 010 SFTs from the LIGO H and L interferometers, respectively. The number of data segments used for the S5R5 run is 121, each spanning no more than T seg = 25 hours and with at least 40 hours of data, as already said. Similar details for the S5R3 run can be found in Appendix B.
The total search frequency range of the S5R5 search is [50, 1 190] Hz, with a frequency resolution δf ∼ 6.7µHz and spindown resolution δḟ ∼ 0.12 nHz s −1 . Each WU analyzes a constant frequency band B 20 mHz, the full spindown interval, ranging roughly from −2 nHz s −1 to 0.11 nHz s −1 , and a region of the sky, as we shall see in Sec. III C.
The original data contained instrumental artifacts in narrow frequency bands that were known before the launch of the Einstein@Home run. Those bands were identified and the corresponding frequency bins in the SFTs were replaced with white Gaussian noise at the same level as the neighboring frequencies. Table I shows which bands were treated in this manner and what instrumental artifact they harbored. These control bands are useful to compare and contrast the results obtained on real data against pure theoretical noise. Measurements and studies after the Einstein@Home run refined the frequencies and widths of these artifacts and identified additional ones; the final lists of artifacts are given in Appendix C, and were used to discard candidates (as we shall see in Sec. IV B). The "cleaning" process affected ∼ 27 Hz of search bandwidth, in addition to bands that were eliminated later in post-processing.
The output data files from each WU are stored as ZIPcompressed ASCII text files containing the 10 000 most  significant candidates ranked according to the significance (as defined in Eq. (21)) over the parameter space searched by that WU. The decision to keep the top 10 000 candidates was based on the maximum upload volume from the hosts to the Einstein@Home servers. All in all, on the order of 10 11 candidates were returned to the Ein-stein@Home server from each run, corresponding roughly to 2.3 TB of data.
The files contain nine quantities for each candidate. The first four are the values (on the coarse grids) of the frequency f , sky-position (α, δ), and spindownḟ . The fifth is the significance of the candidate as defined in Eq. (21). The remaining four quantities are connected with the refined sky-grid centered on each coarse skygrid point (recall that refinement is performed only on the sky). In particular, the location of the most significant point on the fine sky-grid, and the mean and standard deviation of the Hough number count values on all points of the fine grid are returned.

B. Validation of returned candidates
In order to eliminate potential errors, due to defective hardware and/or software or to fraud, BOINC is configured so that each WU is processed redundantly by computers owned by at least two different volunteers. An automated validation process checks the consistency of the results, ruling out those that are inconsistent, in which case new WUs are generated to run again independently. The first step of the validation is to check that the file syntax is correct and that the first four values, i.e. (f, α, δ,ḟ ), are within the appropriate ranges. Next, for the pair of result files from each WU, the validator checks that the values of (f, α, δ,ḟ ) agree to within floating point accuracy (the frequency is in fact checked to double precision). Finally, the significance values CR 1 and CR 2 from the two result files are compared and are validated if In S5R5, about 0.045% of results that were processed by the validator were marked as "invalid", including both syntax errors in individual files and errors in comparisons of different results files. Excluding the syntax errors in individual files, the error rate arising from comparisons of pairs of distinct result files (most likely due to differences in floating point arithmetic on different computational platforms) was ∼ 0.015%. Is it possible that two invalid results could agree with each other and thus end up being marked as valid? While it is difficult to exclude this scenario with complete certainty, an upper limit for the probability of this happening is (0.015/100) 2 ≈ 2.2 × 10 −8 (only the 0.015% error rate due to comparisons of distinct result files is relevant here). As mentioned earlier, there were a total of ∼ 1.1 × 10 7 WUs. It is therefore unlikely that even a single pair of result files would be incorrect and still pass validation.
The threshold of 0.12 on the value of ∆ defined above turns out to be much looser than necessary. The differences in the actual observed values of ∆ from a pair of matching result files are usually much smaller. The observed standard deviation of ∆ turns out to be ∼ 0.012, i.e. an order of magnitude smaller than the threshold. In addition, the standard deviation of the difference CR 1 − CR 2 is measured to be ∼ 0.15, which corresponds to a standard deviation of ∼ 0.7 in terms of the number count. As we shall see later (see e.g. Fig. 3), the loudest events in every 0.5 Hz band have an average loudest number count 70. Thus, these differences correspond to a 1% effect in the number count at the 1-σ level, and we expect this to have a negligible effect on our analysis.

C. Workunit design
The design of the WUs must satisfy certain requirements. The first is that the WUs must be balanced, i.e. each WU must cover the same number of parameter space points so that they can be completed in roughly the same amount of time by a typical host machine. Second, the amount of data that must be downloaded by each host machine and the maximum memory footprint of each job must be within appropriate limits. Finally, one needs to choose the computational time for each WU on a typical host machine and the total time that the project should run, given the total computational power that is available. To meet these requirements, we need to understand how to split up the parameter space, and to measure the CPU core time spent by the search code on each part of the analysis.
We start with the basic parameters of the search, namely the total observation time T obs , the coherent time baseline T seg and the number of segments N seg , the resolution for the coarse and fine sky-grids, dθ F , dθ H , given by Eqs. (12) and (13), respectively and the frequency and spindown resolutions, δf and δḟ , given by Eqs. (9) and (15), respectively. Recall that the limit on the maximum memory footprint of each job already forced us to forego any refinement inḟ .
Unlike the previous Einstein@Home search [24], where each WU searched the whole sky, here we choose each WU to cover a fixed frequency bandwidth B, the entire spindown search range, and a limited area of the sky. Let f b be the highest frequency in the b th search band. We want the computation time for every WU to be approximately the same, hence every WU must search the same number of coarse sky-grid points. Since the resolution in the sky, dθ F , is inversely proportional to the frequency of the signal that we are searching for, WUs at higher frequencies will be searching smaller portions of the sky. Let N b be the number of coarse grid sky points over the whole sky for the b th band: which is shown in Fig. 2 for the S5R5 run. At frequency f b , we then partition the sky in P b parts, each containing N sky points: sky . In practice, to limit the number of sky-grid files needing to be downloaded by the host machines, the sky-grids are constant over 10-Hz frequency bands and are determined based on the highest frequency in each band.
We choose a constant range ofḟ values over the entire frequency band. To see what this implies for potential sources, we define the spindown age for a system emitting at a frequency f and spindownḟ to be τ = f /|ḟ | 2 . It is clear that at a given search frequency, the minimum spindown age is determined by the maximum value of |ḟ | included in the search. Our choice of a constant range ofḟ means that the minimum spindown age is frequency dependent: We choose the minimum age at 50 Hz to be 800 years, and a fixed range inḟ : ∼ [−20, 1.1] × 10 −10 Hz s −1 .
The total number of WUs is simply the sum of the total number of sky-partitions at each frequency f b : In the third step we have replaced the sum over frequencies with an integral between a minimum (f min ) and maximum (f max ) frequency. As a consistency check, we shall see later that B 20 mHz, which is sufficiently small that this is a good approximation.
The total computing time for one WU can be expressed as where N f , Nḟ , and N sky represent the number of coarse search frequency bins, spindown values, and sky-gridpoints respectively, τ F and τ H are the times needed to compute the F-statistics and Hough number count for one point of the coarse grid; N sb 3 represents additional "sideband" bins needed to compute the Houghtransform. The need for these sidebands can be understood by thinking about computing the Hough number 2 The reader is warned that the spindown age as defined here can be rather different from the true chronological age of the star. For a neutron star with a present day frequency f and spindowṅ f due purely to the emission of GWs, the age would be f /4|ḟ |; more generally this depends on the physical mechanisms that are responsible for the spindown [41]. 3 The bins N sb can be calculated by ∆f sb /δf , where the average Hough "sidebands" ∆f sb can be estimated from the Hough count for a frequency near the edge of a given search band: that number count will involve summing the peakgrams along a curved track that can extend a small distance to either side of the target frequency. The time τ F needed to compute the F-statistic for one point of the coarse grid can be expressed as where N SF T is the total number of SFTs used in the search; τ 1 F is the time needed to compute the F-statistic per single SFT. The time τ H in Eq. (28), needed to compute the Hough number count corresponding to one point of the coarse grid, can be written as where τ 1 H is the time to sum the Hough number count per single data segment and per single point of the fine grid. Here N ref sky is an overall refinement factor, i.e. the number of grid points analyzed by the Hough algorithm for each coarse grid point. In our case, since the only refinement is over the sky, N ref sky is given by Eq. (14). The computational time is thus determined by the two timing constants τ 1 F and τ 1 H . For our implementation of the algorithm, these constants were measured to be These numbers are of course only average values for a typical host CPU core available at the time of the Ein-stein@Home runs. The presence of N sb leads to an overhead for the computation. We want to control this overhead and keep it below some acceptable level. Thus, we define the overhead to be the ratio between the time spent in computing the F-statistic for the "sidebands" and the total computational time: The bandwidth B needs to be sufficiently large so that is sufficiently small, but a too large value of B can lead to excessively high download volumes for the Ein-stein@Home clients. From the above equation, we see that the frequency bandwidth B can be determined by fixing : For all the Einstein@Home runs presented here we choose = 5%. For the S5R5 run, this leads to B 20 mHz. The total run-time τ p of the project is where N CPU represents the number of volunteer CPU cores. Given a certain number of CPU cores and having fixed τ p , the maximum search frequency f max can be derived from the above equations to be For the S5R3 and S5R5 runs, the nominal project duration was chosen to be 6 months. With the above choices, the search frequency ranges for S5R3 and S5R5 turn out to be respectively [50, 1 200] Hz and [50, 1 190] Hz.

D. Accuracy of spindown model
Let us briefly discuss the second order spindownf which, as mentioned previously, is not a part of our search. For our frequency resolution δf , given by Eq. (9), and the full observation time T obs ,f would have to be at leastf in order for the signal to move by a single frequency bin over the full observation time. On the other hand, for a minimum spindown age τ min at a frequency f , a useful estimate for the range off that we should search is f /τ 2 min =ḟ 2 max /f . Thus,f is potentially more important at lower frequencies and higher spindown values. Using the maximum value of |ḟ | in the search and the minimum value of the search frequency gives us a value off that might be potentially of astrophysical interest: Comparing with the minimum value off obtained above, we see that there is potentially a region in (f,ḟ ) space where we could improve our astrophysical detection efficiency by includingf ; for our chosen range ofḟ there is no effect off above ∼308 Hz. It is important to note that the calculation of Eq. (36) is too conservative because it does not include any correlations between the phase evolution parameters. On the other hand, there is considerable uncertainty in the value off ast . If the neutron star has a braking index n ≡ ff /ḟ 2 , thenf ast increases by a factor n. If a star is spinning down purely due to gravitational wave emission, then n = 5. On the other hand, for the Vela and Crab pulsars, observed values of n are ∼ 1.4 and ∼ 2.5 respectively [42,43]. The actual impact on our astrophysical reach is thus hard to quantify. Let us consider as an example the extreme case when the spindown is entirely due to gravitational radiation so that the braking index is n = 5. For this case, a conservative estimate of the part of the spindown range |ḟ | cons that is included in our search is: Here we have used Eqs. (36) and (37)  IV.

S5R5 POST-PROCESSING
As said earlier, roughly 10 11 candidates from the S5R5 run were returned to the Einstein@Home server. They were then transferred to the 6720-CPU-core Atlas Computing Cluster [44] at the Albert Einstein Institute in Hannover, and post-processed. The goal is to filter the set of 10 11 candidates, excluding false candidate events. The post-processing strategy consists of the following steps: • selection of 100 most significant candidates in 0.5 Hz frequency bands; • removal of known instrumental noise artifacts; • removal of unknown data artifacts through the Fstatistic consistency veto; • follow-up of the most significant candidates with S5R3 data; • fully-coherent follow-up of the surviving candidates.
The items outlined above are described in the next subsections.

A. Selecting the top candidates in frequency bands
As is commonly done in CW searches, the results are examined separately in fixed-size search frequency bands; here we choose to perform the analysis in 0.5 Hz bands. As described earlier, in designing the WUs, we have previously been led to break up the frequency range in ∼ 20 mHz bands and the sky has been partitioned as well. This was however done for purely technical reasons to make the search on Einstein@Home feasible. The choice of frequency bands for the post-processing is based on different requirements. First, we would like the detector to have roughly constant sensitivity within each frequency band. Furthermore, as we shall see, the search does not result in a convincing detection candidate, and upper limits will be set over each of these frequency bands. Having a large number of very narrow bands would make the calculation of the upper limit very computationally intensive. The choice of 0.5 Hz is a compromise between these two requirements. This choice is in fact comparable to previous CW searches and will make comparisons straightforward. Finally, we note that all other things being equal, having a larger frequency band will in principle also lead to a decrease in sensitivity simply because of having a larger number of templates. This is however a relatively minor effect in the present case.
For each of the 0.5 Hz-wide frequency bands, we select the 100 most significant candidates for further analysis, leading to a set of 228 000 loudest S5R5 candidates. This choice was dictated by the available computational and human resources for the post-processing. As will be illustrated in the following, at the end of the automated post-processing procedure there will remain of order 10 candidates which survive all selection criteria. This is about the number that we can afford to follow-up manually with further investigations. As our follow-up procedures are further automated and optimized, it will become possible to consider lower thresholds and to inspect a correspondingly larger number of candidates.
The number count of such candidates is plotted in Fig. 3 as a function of the frequency. For most bands, it is consistent with expectations. The number count generally increases with increasing frequency because the number of sky points searched over increases (see Eq. (11)) and the maximum expected value of a random variable, over repeated trials, grows with the number of independent trials. Let us ignore the effect of the weights, and take the Hough number count n c (defined in Eq. (17)) to be an integer random variable following a binomial distribution. The cumulative probability for obtaining n c or lower is The binomial parameters are N seg = 121 and p seg = 0.267 (consistent with a 2F threshold of 5.2). The probability p max that the maximum number count is n c over a set of N trials independent trials is We compute Eq. (40) as a function of N trials , which we take to be the number of templates searched to cover 0.5 Hz bands, i.e. 7.5 × 10 4 frequency values × 18 spindown values × 8 444 Hough pixels per F sky point × N b , the total number of F sky points shown in Fig. 2. The number of Hough pixels has been computed using Eq. (14) with m = 0.3 and ℘ = 0.5. Figure 3 shows the expected value of the maximum (central black curve), computed using the probability function given by Eq. (40), superimposed on our measurements (red circles) and confirming that there is broad agreement, in most bands, between our results and the expectations for Gaussian noise.

B. Removing known data artifacts
As a first step of the post-processing pipeline, we eliminated from the list of top candidates any candidate whose frequency was too close to that of either a known artifact or to the cleaned noise bands described in section III A. Specifically, we discarded those candidates whose detection statistic could have been constructed with contributions either from: • data polluted in either of the two instruments by spurious disturbances; details of such detector disturbances are given in Appendix C and, in particular, a list of known spectral disturbances for the H and L instruments are listed in Tables VI and VII.
• from fake noise that had been inserted by the cleaning process and hence could not host a CW signal (see Table I in Sec. III A).
After this veto, about 25% of the candidates were eliminated from the original set of 228 000 loudest S5R5 candidates. More precisely, a total of 172 038 S5R5 candidates survived this veto. The bandwidth removed due to the lines listed in Appendix C amounted to ∼243 Hz; an additional 27 Hz was removed due to the cleaned noise bands.

C. The F-statistic consistency veto
We have thus far considered only known instrumental disturbances for vetoing candidates. However, we expect there to be more such disturbances present in the data that have not yet been explicitly identified. The idea is to discriminate between disturbances in a single detector and signals, which should produce consistent values of the F-statistic in both detectors [45]. We refer to this method as the F-statistic consistency veto.
For each of the 172 038 S5R5 surviving candidates, the single-detector and multi-detector F-statistic was computed for each of the 121 data segments and then averaged over the segments. We refer to these averaged 2F values as 2F H and 2F L for the H and L detectors respectively, and 2F HL for the coherent combination of the data from the two detectors. Candidates were discarded if either 2F H or 2F L were greater than 2F HL . Using this veto, a small fraction (4.1%) of candidates was eliminated, leaving 164 971 surviving S5R5 candidates.
The impact of this veto is limited in this case due to the prior removal of the bulk of instrumental artifacts. However, the F-statistic consistency veto represents an efficient method to remove disturbances that clearly stand out of the noise in the absence of independent instrumental evidence. Figure 4 shows the average 2F-values

D. Distribution of candidates
We have now applied all of our vetoes that try to remove instrumental artifacts. While there will of course remain other low amplitude instrumental spectral lines and hardware signal injections (described later in Sec. VI), we now need to deal with the possibility that, say, even Gaussian noise can mimic a signal in some cases. All remaining candidates will need to undergo detailed individual inspection and we will only be able to afford this for a few candidates. As our follow-up techniques become more refined, optimized and automated, we will be able to improve this part of the pipeline and dig deeper into the noise. Figure 5 shows the histogram of 2F HL for the 164 971 surviving candidates up to 2F HL -values of 9. The distribution actually extends up to 2F HL ∼ 542.8 but we show only the low 2F HL -values distribution in order to explain our next choice of threshold. There are 166 candidates with 2F HL > 9 and they are all clustered at the two frequencies of ∼108.9 Hz and ∼575.2 Hz, corresponding to two simulated signals injected in the data stream, as discussed in Sec. VI. The region 2F HL < 6.5 contains well over 99% of the candidates and, as seen in Figure 5, below 6.5 the density of candidates increases very sharply. We will take 6.5 as a threshold for the next step in our follow-up procedure.
The number of candidates expected to survive the 6.5 cut on 2F HL in fixed 0.5 Hz bands increases with frequency because the number of sky locations searched scales with the square of the searched frequency. In or- Values of 2F averaged over 121 data segments for the single-detector case, 2FH (green dots), 2FL (blue dots) and the multi-detector case, 2FHL (red dots), against those for the combined multi-detector statistic. The top (bottom) plot shows such values for 164 971 (7 067) surviving (vetoed) S5R5 candidates such that 2FH and (or) 2FL is less (greater) than 2FHL . der to compute the false alarm probability corresponding to this threshold in different frequency bands, note that in the absence of a signal, the value of 2F HL in the i th segment, 2F (i) HL , follows a χ 2 distribution with 4 degrees of freedom. Furthermore, since it is clear that 2F HL × 121 is a χ 2 random variable with (4 × 121) degrees of freedom. The false alarm probability corresponding to a threshold at (6.5×121) for such a random variable is ∼10 −16 . This corresponds to expected false alarm rates of about 0.1%, 0.6%, 2.6% and 10% for searches in Gaussian noise over over 0.5 Hz bands at 100 Hz, 250 Hz, 500 Hz and 1 000 Hz, respectively, considering the number of independent trials given by the number of searched templates in the respective bands. Disregarding the non-Gaussian line features evident in Fig. 6, the ratio of the number of candidates observed above the 6.5 threshold at lower frequencies (say below 800 Hz) to that at higher frequencies (say above 800 Hz) is not inconsistent 4 with the ratios of the false alarm rates computed above. We note that the false alarm probability given above overestimates the number of expected candidates above the 6.5 threshold because that threshold is not the only cut applied to the data. The previous cuts, discussed in the preceding sections, lower the actual false alarm probability of the surviving candidates. There are 184 remaining S5R5 candidates, for which 2F HL > 6.5 and they are shown in Fig. 7 as a function of the frequency. They are clustered at twelve frequencies and only the most significant candidate from every cluster has been followed up.

E. Following up candidates with S5R3 data
For the next step in our post-processing pipeline, recall that our underlying signal model of Eq. (1) assumes that the signal is long-lived; thus its amplitude is constant in time and its intrinsic frequency evolves smoothly according to Eq. (5) with a constant spindown. This is an idealized model: although pulsars are the most stable clocks in the Universe, neutron stars are nonetheless known to glitch, to be perturbed by external agents, and in some cases to be affected by significant timing noise. Furthermore, for sufficiently long observation times, the spin-down evolution model that we use here, including only the first spindown order, may not be adequate to describe the actual signal model (see also the discussion in Section III D). However, since the data set used in the S5R3 run ends just about a week before the S5R5 data set, it is reasonable to assume that any putative signals should be present in both data sets. Moreover, the average noise floor level in the detectors turns out to be approximately stable between S5R3 and S5R5. Thus, we might expect that a detectable signal should be visible in both searches. The next follow-up step for each candidate then consists of a hierarchical search carried out on the same WU as done for S5R5 (i.e. over the same parameter space), but using the S5R3 data set. The closest 5 candidate to the original one from such a search was identified and the value of its detection statistic was compared with what one would expect if the S5R5 candidate were due to a signal. In particular, the expected number counts in S5R3 and S5R5 should be related according to where 84 is the number of data segments used in S5R3, as we shall see in Appendix B. Possible reasons for this not to be a good approximation would be if the detector noise floor were to vary significantly between S5R3 and S5R5, or the fact that the relative geometry between the detector and source varies in time. We have already remarked that the noise floor is, on the average, stable between S5R3 and S5R5. Furthermore, albeit in each segment 5 The distance used to judge closeness between candidates is a Euclidean distance expressed in bins in the four dimensions (f,ḟ , α, δ). were discarded as not being consistent with a CW signal, where σ was computed using Eq. (22). As shown in Table II, two candidates were discarded by this followup test, at ∼80.9 Hz and ∼108.9 Hz. However, the second of these, as well as the candidates at ∼52.8 Hz and ∼575.2 Hz, represent three simulated signals injected only part of the time during S5, as discussed in Sec. VI.
As we can see from Fig. 7, the bulk of the candidates arise from the strong hardware-injected pulsar 2 and 3 signals (see Sec. VI). In particular, 84 and 88 candidates are clustered near the frequencies of these two injected signals respectively.

F. Fully-coherent follow-up
Excluding the hardware-injected simulated signals, the post-processing up to this point has left us with 8 surviving candidates. These have been significant enough to pass our thresholds and have not been clearly identified as instrumental artifacts, or eliminated by inconsistency between the H and L detectors, or by the follow-up with the S5R3 data set. We therefore need to consider other more sensitive methods. If these candidates are real signals, then their SNRs and significance should increase if the parameter space grids are made finer, or as the

coherent integration time becomes larger.
We use a three-step procedure consisting of a gridbased semi-coherent Hough search, followed by a semicoherent and a final fully-coherent F-statistic search, using the Mesh Adaptive Direct Search (MADS) algorithm for constrained optimization. The reference implementation of the MADS algorithm is publicly available through the NOMAD library [46,47]. Hence, in the following, we refer to such searches simply as NOMAD searches. Contrary to the traditional grid-based methods, a mesh adaptive search constructs the trial points as the search evolves aiming to find the maximum of the statistic.
The three steps of the follow-up procedure are the following: 1. Re-run the Hough search around a given candidate, but with a finer grid to reduce the mismatch with a putative signal. The search region includes 5 frequency bins on either side around the candidate and 16 neighboring coarse sky-grid points. The fine Hough sky-grid is refined by a factor of 2 in each direction by using ℘ = 1 (see Eq. 13) instead of 0.5 as in the original search. Furthermore, we refine the coarseḟ grid spacing of Eq. (10) by a factor N seg = 121 [19].
2. The loudest candidate from the first step is used as a starting point for the semi-coherent F-statistic NOMAD optimization. The detection statistic in this step is the sum of the F-statistic values from each segment. This search has been performed in a fixed parameter space box around the starting point. The dimensions of the box are ∆f = 10 −4 Hz, ∆α = 0.10 rad, ∆δ = 0.24 rad and ∆ḟ = 10 −10 Hz s −1 . The loudest candidate found in this semicoherent F-statistic NOMAD search is passed on to the next step.
3. In the third step, the loudest candidate from the previous step is used as a starting point for the fully-coherent F-statistic NOMAD search. This search spans the entire duration of the S5R5 data set and has been carried out in a parameter space box defined by using the diagonal elements of the inverse Fisher matrix in each dimension around the starting point. These elements are described by Eq. (16) in [48] and have been computed from the inverse of the semi-coherent parameter space metric computed at the candidate point, re-scaled by the measured SNR at the same point; for more details we refer the reader to [48].
In both the semi-coherent and fully-coherent NOMAD searches, we ran multiple instances of the algorithm iterating over the mesh coarsening exponent using both deterministic [46] and stochastic [49] methods for the choice of search directions. Based on Monte-Carlo studies, the false-dismissal probability of the follow-up procedure is found to be less than 10%.
As said earlier, in the presence of a real signal we expect the significance of a candidate to increase as the template grid becomes finer because there will be a template with a smaller mismatch with respect to the real signal. At that template the signature of the signal should be more evident and all the consistency tests should continue to hold. If the candidate signal detected on the finer grid does not pass a consistency test, this indicates that it is not behaving as we would expect from the signals that we are targeting. The candidates at ∼96.6 Hz, ∼144.7 Hz, ∼932.4 Hz, ∼1030.2 Hz and ∼1142 Hz fail a multi-detector versus single-detector F-statistic consistency test (see Sec. IV C) after performing the semicoherent NOMAD search; therefore, they cannot be considered defensible CW signals. Moreover, line artifacts appear in the average power spectrum of S5 H data at ∼932.4 Hz, ∼1030.2 Hz and ∼1142 Hz.
The remaining candidates, namely at ∼434.1 Hz, ∼677.5 Hz and ∼984.4 Hz, survive the F-statistic consistency test on the finer grid and are followed-up with the fully-coherent F-statistic NOMAD search. However, for each of them, the maximum value of the detection statistic over the parameter space searched is much lower than would be expected based on the original candidate parameters, and in fact is consistent with the expectation for Gaussian noise. Hence, also these three candidates do not survive a more sensitive inspection and cannot be considered viable detection candidates. Thus, we see that all the candidates listed in Table II are inconsistent with the properties of a true CW signal.

V. UPPER LIMIT ESTIMATION
The analysis of the Einstein@Home searches presented here has not identified any convincing CW signal. Hence, we proceed to set upper limits on the maximum intrinsic gravitational wave strain h 0 that is consistent with our observations for a population of CW signals described by Eq. (8), from random positions in the sky, in the gravitational wave frequency range [50.5, 1 190] Hz, and with spindown values in the range of ∼ [−20, 1.1] × 10 −10 Hz s −1 . The nuisance parameters cos ι, φ 0 , and ψ are assumed to be uniformly distributed. As commonly done in all-sky, all-frequency searches, the upper limits are given in different frequency sub-bands and here we have chosen these to be 0.5 Hz wide. Each upper limit is based on the most significant event from the S5R5 search in its 0.5 Hz band.

A. Monte-Carlo upper-limit estimates
Our procedure for setting upper limits uses Monte-Carlo signal injection studies using the same search and post-processing pipeline (except for the S5R3 and fullycoherent follow-ups) that we have described above. In every 0.5 Hz band, our goal is to find the value of h 0 (denoted h 90% 0 ) such that 90% of the signal injections at this amplitude would be recovered by our search and are more significant than the most significant candidate from our actual search in that band. We can thus exclude, with 90% confidence, the existence of sources (from our specific population) that have an amplitude h 0 > h 90% 0 .
In order to estimate h 90% 0 , for each injection at a randomly chosen parameter space point, a hierarchical search is performed over a small parameter space region, which consists of: • a 0.8 mHz frequency band centered at the S5R5 frequency grid-point closest to the randomly chosen source frequency; • 4 spindown values around the S5R5 frequency derivative grid-point closest to the randomly chosen spindown; • a sky-patch consisting of 10 S5R5 coarse sky-gridpoints closest to the randomly chosen sky location.
At the end of each hierarchical search, the most significant candidate is selected and post-processed as described in the previous section. The vetoes for excluding known instrumental lines are not required because we have already excluded them from the upperlimit analysis. We first check if this candidate is significant enough to be part of the 100 Einstein@Home loudest candidates originally selected in its corresponding 0.5 Hz band. Then we perform the F-statistic consistency check and finally we compare the computed average multi-interferometer 2F-value ( 2F Cand HL ) with the maximum 2F HL value we have in the corresponding 0.5 Hz band. If 2F Cand HL is greater than the maximum 2F HL , then the simulated source is considered to be recovered and more significant than the most significant candidate of the search. The confidence level is defined as C = n rec /n tot , where n rec is the number of recovered candidates, and n tot = 100 is the total number of injections performed.
After some preliminary tuning to determine a range of h 0 values close to the 90% confidence level, we use an iterative procedure to determine the confidence as a function of the injected population h 0 until we hit a confidence value close to 90%, within the expected 1 σ fluctuations. Since we use 100 injections, from a binomial statistic we estimate the 1 σ fluctuation to be 3% and hence we associate the h 90% 0 value to any measured confidence in the range 87% − 93%. The 3% uncertainty in confidence translates to an uncertainty in h 90% 0 smaller than 5%, as can be seen from Fig. 8, which shows a typical confidence versus injected h 0 behavior. Each point in Fig. 8 was derived with 1 000 injections and hence is affected by fluctuations smaller than 1%. The lower (red) curve in Fig. 9 shows the resulting upper limits as a function of the frequency. The upper (blue) curve shows the upper limit values from the previous Einstein@Home search in early S5 data [24]. The current upper limit values are about a factor 3 more constraining than the previous Einstein@Home ones. In particular, the most constraining upper limit falls in the 0.5 Hz-wide band at 152.5 Hz, where we can exclude the presence of signals with h 0 greater than 7.6 × 10 −25 . The three stars shown in Fig. 9 correspond to the simulated pulsars 2, 3 and 5, i.e. the hardware injections recovered in the S5R5 search (discussed in Sec. VI).
The numerical data for the plot in Fig. 9 can be obtained separately 6 . A conservative estimate of the overall uncertainty on the h 90% FIG. 9. Upper limits for the S5R5 Einstein@Home search (red dots) as well as the previous Einstein@Home search, called S5R1, which used early S5 data (blue dots) [24]. The three stars correspond to hardware-injected simulated pulsars which were recovered in the S5R5 search. The curves represent the source strain amplitude h0 at which 90% of simulated signals would be detected. The vertical bars represent 156 half-Hz frequency bands contaminated by instrumental disturbances for which no upper limits are provided. The upper limits for the 0.5 Hz-wide bands starting at 69.5 Hz, 139.5 Hz, 335.5 Hz and 648 Hz are fairly high due to significant partial contamination in these bands by lines listed in Tables VI and VII. Note that the broadness of the red curve is due to the 5% steps used to vary the injected population h0 values in the Monte-Carlo signal software-injections until a confidence value close to 90% is reached. In addition, less than 1/4 of the spectral range shown was excluded in many narrow bands because of known instrumental artifacts, as described in Sec. IV B. The cyan curve shows the predicted h 90% 0 upper limits according to Eq. (43).
having added to the 1 σ statistical upper limit estimation procedure uncertainty the 10% amplitude calibration uncertainties for the data used in this Einstein@Home run [50].
As we have excluded from the search those frequency bands hosting spectral artifacts (Tables VI and VII) and the cleaned noise bands (Table I), we therefore also exclude these frequency intervals from the upper limit statements. Vertical bars in Fig. 9 represent 156 half-Hz frequency bands for which no upper limits are provided because the entire half-Hz band has been excluded.
As shown in Fig. 9, the upper limits on h 0 provided in the 0.5 Hz-wide bands starting at 69.5 Hz, 139.5 Hz, 335.5 Hz and 648 Hz are fairly high, roughly equal to 3 × 10 −23 , 8.8×10 −24 , 1.7×10 −23 and 2×10 −23 , respectively. This is due to significant partial contamination in these bands by lines listed in Tables VI and VII; the upper limit is given for the remaining, clean part of the band, but loud candidates from the disturbed part make up the loudest 100 candidates selected in the processing, so a simulated signal must be especially loud to surpass those. Note that, for the same reason, if we had set upper limits in the 156 half-Hz bands shown in Fig. 9 we would have obtained similarly high upper limits on h 0 .

B. Analytic sensitivity estimates
The h 90% 0 upper limits can be independently predicted using the method in [51], adapted to the Hough-on-F-statistic search method (see [52] for details). The upper limit procedure described above is modelled by a simple threshold on the number count, where the thresholds are the largest number counts observed in each 0.5 Hz upper limit band. The probability that, in the neighbourhood of an injected signal, the number count n c will exceed a threshold n c,th is denoted by P [n c > n c,th |ρ(h 0 , cos ι, φ 0 , ψ, m)]; this probability can be calculated analytically from the known distribution of n c . The recovered SNR, ρ, is a function of the nuisance parameters and of the mismatch m between the injected signal and the nearest template. In addition, note that in the presence of a signal, 2F follows a non-central χ 2 distribution with 4 degrees of freedom; ρ 2 is the non-centrality parameter of this distribution. In the presence of a signal, averaging P over the parameters of ρ (except h 0 ) gives P (n c > n c,th |ρ(h 0 )) , which equals the confidence of recovering a population of injections with amplitude h 0 . In each 0.5 Hz band, we determine the value of h 0 such that P (n c > n c,th |ρ(h 0 )) = 90%; this value is then the predicted value of h 90% 0 , given by as a function of the detector noise S h and the total data volume T data = N seg T seg . The factor H varies from ∼141 to ∼ 150 over the range of search frequencies, and is plotted in [52]. It is given by H = 2.5ρ 90% N seg , wherê ρ 90% is the mean injected SNR per segment of a population of signals as described above, and is itself a function of the false alarm and false dismissal probabilities, and N seg [51,52]. The variation of H as a function of frequency arises from the variation ofρ 90% as a function of the false alarm probability in each upper limit band, which are calculated from the largest number counts n c plotted in Fig. 3.
The predicted values are shown by the cyan curve in Fig. 9. The root-mean-square error between the Monte-Carlo estimated and predicted h 90% 0 values (∼ 7% over all frequencies) is comparable with the uncertainties due to calibration (10%) and the finite stepping of the Monte Carlo procedure (5%). This demonstrates that the sensitivity of the Hough-on-F-statistic search method, as a function of the search parameters, is well understood. Figure 10 shows the maximum reach of our search. The top panel shows the maximum distance at which we could have detected a source emitting a CW signal with strain amplitude h 90% 0 . The source is assumed to be spinning down at the maximum rate considered in the search, ∼ −2 nHz s −1 , and emitting at the spindown limit, i.e. with all of the lost rotational energy going into gravitational waves. The intrinsic gravitational wave strain from a source at a distance d, with frequency f , frequency FIG. 10. Panel (a) and (b) represent the distance range (in kpc) and the maximum ellipticity, respectively, as a function of the frequency. Both the panels are valid for neutron stars spinning down solely due to gravitational radiation and assuming a spindown value of ∼ −2 nHz s −1 . In these plots, the 156 half-Hz frequency bands for which no upper limits are provided have not been considered. derivativeḟ , and emitting at the spindown limit is

C. Astrophysical reach
where the canonical value of 10 38 kg m 2 is assumed for I zz in Eq. (45) Around the frequency of greatest sensitivity, 152.5 Hz, we are sensitive to objects as far as 3.8 kpc and with an ellipticity ε ∼ 10 −4 . Normal neutron stars are expected to have ε less than a few times 10 −6 based on theoretical predictions [31]. A plausible value of ε ∼ 3.5 × 10 −6 could be detectable by a search like this if the object were emitting at 625 Hz and at a distance no further than 500 pc.

VI. STUDY OF HARDWARE-INJECTED SIGNALS
As part of the testing and validation of search pipelines and analysis codes, simulated signals are added into the interferometer length control system to produce mirror motions similar to what would be generated if a gravitational wave signal were present. Table III shows the parameters of the set of simulated CW signals injected into the LIGO detectors; we shall often refer to these injections also as "fake pulsars". These injections were active from the GPS epoch 829 412 600 s until 875 301 345 s. Of these ten hardware-injected CW signals, eight had frequencies covered by the S5R5 search frequency band: the fake pulsars 4 and 7 have frequencies outside this band and thus have not been taken into account during this analysis. The fake pulsars 6 and 8 have spindown values outside the S5R5 search frequency derivative range.
As a minor complication, the hardware injections were not active all the time. In the S5R5 data set, their duty cycle was ∼ 63% and ∼ 60% in L and H, respectively. The hardware injections were active in 76 of the 121 S5R5 data segments, and they were completely absent in the remaining 45 segments. Their expected value of the Hough number count is thus: where the superscript "HI" refers to Hardware Injection, j runs over the number of data segments where the hardware injections were active and l runs over the remaining segments, and the w j are the Hough weights given by Eq. (18); η j and p are the probabilities that the estimated value of 2F (for a given data segment) crosses the threshold 2F th = 5.2 in the presence and absence of a signal, respectively (expressions for η j and p can be found in [12,19]). Figure 11 shows, for each pulsar hardware injection, the histograms of the expected 2F-values, computed for all the segments where a particular hardware injection was active. From this figure one can infer that the fake pulsars expected to be the loudest in terms of 2F-values are pulsars 2, 3, 5 and 8. The pulsar hardware injections go through the normal search and post-processing pipelines in the usual way as described in the previous sections. As discussed in the following, the fake pulsars recovered in the S5R5 search are pulsars 2, 3 and 5. On the other hand, fake pulsars 6 and 8 are missed because their spindown values are outside our search range. As expected, fake pulsars 0, 1 and 9 are not recovered because their amplitudes are too weak and they do not pass the 2F HL > 6.5 cut, as shown in Table IV. In this table, the observed number counts of 69 to 72 are consistent with noise fluctuations in those half-Hz bands (see Fig. 3) superseding the weak injected signals. Note that the expected 2F-values shown in Fig. 11 were computed assuming a search using exactly the correct signal parameters provided in Table III, while the observed 2F HL -values in Table IV were obtained from our search which uses a grid of templates, so significant mismatch is to be expected. Table IV compares the expected and observed values of the number counts associated with the S5 hardware injections and the surviving S5R5 candidate events closest to them. As in Sec. IV E, the measure of distance used here is a Euclidean distance, expressed in bins, in the four dimensions (f,ḟ , α, δ). Table V lists the parameters of pulsars 2, 3 and 5 and the parameters of the corresponding recovered candidates. We successfully find candidates near the correct signal parameters. The mismatch in spindown might seem large, but in fact the injections were found at the nearest spindown template. The number count values show consistency within the 3σ range.  , δ), the polarization angle ψ, the initial phase φ0, the inclination parameter cos ι and the dimensionless strain amplitude h0. These parameters correspond to the only set of hardware injections, injected into the S5 LIGO data, that fall within the GPS times of the S5R5 data. respectively) associated with the hardware injections and the recovered S5R5 candidates closest to these injections. The 2FHL -values for each of the candidates are also listed. The fake signals labeled pulsar 4 and 7 were not taken into account in this analysis since they have frequencies outside the S5R5 search frequency range. The "expected" values marked by asterisks are not actually expected to be obtained because the spindown rates for those signals lie outside the range of this search.

VII. CONCLUSIONS
No evidence for continuous gravitational waves has been observed in the search presented here. Upper limits on the intrinsic gravitational wave strain have been derived using standard population-based methods, and are shown in Fig. 9. These results are about a factor of 3 more constraining than those from the previous Ein-stein@Home search in early S5 data [24]. This improvement comes from using more data (a year versus two months), from using a multi-detector coherent statistic (versus a single-detector statistic), from a lower effective threshold in the coherent detection stage, and from a more sensitive incoherent method to combine the information from the coherently analyzed segments. The largest effect comes from lowering the effective threshold. Indeed much of the improvement in sensitivity is attributable to improved data analysis methods (as opposed to improved detector sensitivity). If we had used the much higher threshold of 25 on 2F, as in [24], our sensitivity would have been a factor of ∼ 2.5 worse than our final upper limits, thereby undoing almost all of the factor of 3 improvement mentioned above. We have not included second time-derivatives of the frequency in our search. This could be astrophysically significant in some regions of parameter space, as discussed in Sec. III D. It is important to keep this caveat in mind while interpreting our results. This is the most sensitive wide-frequency-range, allsky search for CW signals performed to date. The upper limit values are comparable to those obtained recently using the PowerFlux method [22,26] on the entire S5 data set (S5R3+S5R5). Ref. [26] searched for CW signals over the whole sky, in a smaller frequency band (up to 800 Hz versus 1 190 Hz here), but a broader spindown range up to −6 nHz/s. Strain upper limits were set at the 95% confidence in 0.25-Hz wide sub-bands. In particular, near 152 Hz, the PowerFlux strict, all-sky upper limit on worst-case linearly (best-case circularly) polarized strain amplitude h 0 is ∼ 1 × 10 −24 (3.5 × 10 −25 ). As a comparison, at the same frequency, this search constrains the strain to h 0 7.6×10 −25 (as shown in Fig. 9), 9.2×10 −25 and 3.2×10 −25 for the case of average, linear and circular polarization, respectively, with a 90% confidence level in a 0.5-Hz wide band.
It has long been expected that searching a large parameter space for CW signals will require hierarchical semi-coherent searches. This analysis is a milestone towards that goal, and we expect that future analyses will build on the tools developed here.
There are a number of areas where further improvements are possible. In the latest round of analysis (an Einstein@Home processing run that began in March 2012), some of the post-processing techniques developed for this analysis have been "moved upstream" to the hosts. One example is the generalized F-statistic consistency test [53]. This continues the pattern of moving analyses formerly carried out in the post-processing stage (for example, the incoherent combination step) onto host machines. Another step forward is in the semi-coherent algorithm that combines the coherent analyses from the different segments. The Hough algorithm described here turned out to be rather cumbersome, and does not combine the coarse and fine grids in an optimal way. The latest round of Einstein@Home processing makes use of a simpler optimal semi-coherent method, which allows longer coherent time baselines to be used. This method, based on a detailed analysis of correlations in parameter space [54], is described in [55]. Looking farther forward, we expect to use higher-order spindown parameters both to search for a broader class of signals as well as to be able to employ longer coherent time baselines in the analysis.
The Advanced LIGO and Advanced Virgo detectors are currently under construction, and should begin operations around 2016. In comparison with the current generation, these instruments will provide an order-ofmagnitude improvement in strain sensitivity, increasing the volume of space observed by a factor of a thousand.
These and other improvements in data analysis methods and instrumentation make us optimistic that we will eventually be able to make direct detections of CW signals. Such detection will provide new insights into the internal structure, formation history and population statistics of neutron stars.
is not necessary to step through parameter space pointby-point to calculate the final number count. For the vast majority of cases, our implementation of the Houghtransform agrees with the brute force approach for calculating the number count, but the two can occasionally differ. If we were not using weights, these differences would have a minor effect on the number count. Occasionally, however, these floating point errors coincide with the cases when we assign excessively large weights to particular segments as discussed above. In these cases, the discrepancies in the number count can be large and in some cases may even exceed the number of segments, which is in principle a strict upper bound on the number count. Note that our upper limits remain valid because they consistently use the same search code, and any candidates are followed-up by independent codes thereby avoiding spurious false alarms. Using the modified weights based on Eq. (19) fixes this problem as well.  [50, 1 200] Hz, with a frequency resolution δf ∼ 6.7 × 10 −6 Hz and a mismatch value m = 0.3, leading to the spindown resolution given by Eq. (16). In the S5R3 run, each WU analyzed a constant frequency bandwidth B 16 mHz, the full spindown interval, ranging roughly from −1.6 nHz s −1 to −3.1 × 10 −11 Hz s −1 in steps of 3.8 × 10 −10 Hz s −1 , and a certain region of the sky, as already discussed in Sec. III C. The sky-grids and output file formats are identical to those used in S5R5. A major difference with S5R5 are the weights: no weighting scheme was used in S5R3 and all the weights w i appearing in Eq. (17) were set to unity. In total, ∼ 7 × 10 10 S5R3 candidates were sent back to the Einstein@Home server. The S5R3 post-processing was performed before that for S5R5, and it was used as a "test-bed" for the latter. It consists of the same series of steps described previously in Sec. IV. The significance values of the 100 loudest candidates in 0.5 Hz-wide frequency bands are plotted in Fig. 12 as a function of the frequency. These represent a total of 230 000 candidates. This set is reduced by ∼ 27% via the removal of instrumental noise artifacts (listed in Tables VI and VII) and of the candidates from search frequency bands close to the fake noise (according to Table I). For all the surviving 167 779 candidates, the F-statistic consistency veto described in Sec. IV C removes an additional 3.6%. Figure 13 shows the values of 2F H , 2F L and 2F HL versus 2F HL  Candidates whose 2F HL -value is greater than 6.5, i.e. 1 465 out of 161 785, are followed-up by performing a hierarchical search using the S5R5 data set. This run, whose details are described in Section IV, consists of ∼ 46% more data than was used for S5R3 and is thus more sensitive than S5R5. The followed-up candidates are plotted in Fig. 14 as a function of the frequency: 87 candidates are clustered at ∼ 108.857 Hz and represent the contribution of the hardware injection pulsar 3 [56]; 73 candidates have frequencies peaked around ∼ 1081 Hz, 59 candidates at ∼ 1198.9 Hz and 28 at ∼ 1141 Hz. Except for the candidates that represent the contribution from the hardware injected signal, none of the S5R5 candidates is consistent with a signal at a level consistent with the observed S5R3 excess. We then conclude that no gravitational wave signal is observed in the S5R3 data.

Appendix C: Instrumental noise artifacts
This appendix contains lists of the main known spectral lines of instrumental origin in the LIGO detectors during the S5 run. They were individually identified with a particular source, or were members of identified combs of lines found in many channels (e.g. 60 Hz combs), or the source was unknown, but they were found in frequency coincidence between the gravitational wave channel and auxiliary channels. In the latter case, to ensure that actual gravitational wave signals were not rejected, the SNR in the auxiliary channel had to be at least 5 times larger than in the gravitational wave channel and the density of lines in the auxiliary channel had to be low enough that accidental coincidence with a line in the gravitational wave channel was highly unlikely. The spectral lines and harmonics detected in H and L are listed in Tables VI and VII, respectively. As mentioned earlier, on the basis of these tables, about 22% and 25% of candidates have been excluded from the analysis in the S5R5 and S5R3 post-processing, respectively.
For each candidate with frequency f c and derivativė f c , the candidate was rejected if a band ∆f c f c × 10 −4 + |ḟ c | × 10 7 s on either side of the signal had any overlap with an instrumental line band. This was done in order to take into account the maximum possible Doppler shift due to the Earth's orbital velocity, which is roughly 10 −4 in units of the speed of light, and the maximum frequency shift due to the spindown of the source over the ∼ ±10 7 s time span relative to the reference time during the Einstein@Home run.