Hierarchical follow-up of sub-threshold candidates of an all-sky Einstein@Home search for continuous gravitational waves on LIGO sixth science run data

We report results of an all-sky search for periodic gravitational waves with frequency between 50 and 510 Hz from isolated compact objects, i.e. neutron stars. A new hierarchical multi-stage approach is taken, supported by the computing power of the Einstein@Home project, allowing to probe more deeply than ever before. 16 million sub-threshold candidates from the initial search [LVC,arXiv:1606.09619] are followed up in three stages. None of those candidates is consistent with an isolated gravitational wave emitter, and 90% confidence level upper limits are placed on the amplitudes of continuous waves from the target population. Between 170.5 and 171 Hz we set the most constraining 90% confidence upper limit on the strain amplitude h0 at 4.3x10-25 , while at the high end of our frequency range we achieve an upper limit of 7.6x10-25. These are the most constraining all-sky upper limits to date and constrain the ellipticity of rotating compact objects emitting at 300 Hz at a distance D to less than 6x10-7 [d/100pc].


I. INTRODUCTION
The beauty of continuous signals is that, even if a candidate is not significant enough to be recognised as a real signal after a first semi-coherent search, it is still possible to improve its significance to the level necessary to claim a detection after a series of follow-up searches. Hierarchical approaches were first proposed in the late 90s and developed over a number of searches on LIGO data: [2] and [3] detail a semi-coherent search plus a three-stage follow-up of order 100 candidates; [4] and [5] detail a semi-coherent search plus a series of vetoes and a final coherent follow-up of over 1000 candidates. The search detailed here follows up 16 million candidates and is the first large-scale hierarchical search ever done.
We use a hierarchical approach consisting of 4 stages applied to the processed results ("Stage 0") of an initial search [1]. At each stage a semi-coherent search is performed and the top ranking cells in parameter space (also referred to as "candidates") are marked and are searched in the next stage. At each stage the significance of a cell harbouring a real signal would increase with respect to the significance it had in the previous stage. The significance of a cell that did not contain a signal on the a email: maria.alessandra.papa@aei.mpg.de other hand is not expected to increase consistently over the different stages. In the first three stages the thresholds that define the top ranking cells are low enough that many false alarms are expected over the large parameter space that was searched. And indeed at the end of the first stage we have 16 million candidates. At the end of the second stage we have 5 million . At the end of the third stage we have 1 million . At the end of the fourth stage we are left with only 10 candidates.
The paper is organised very simply: Section II introduces the quantities that characterise each stage of the follow up. Section III illustrates how the different stages were set up for the S6 LIGO Einstein@Home candidates follow-ups. Section IV present the gravitational wave amplitude and ellipticity upper limit results. In the last section, Section V, we summarise the main findings and discuss prospects for this type of search.

II. QUANTITIES DEFINING EACH STAGE
From one stage to the next in this hierarchical scheme, the number of surviving candidates is reduced, the uncertainty over the signal parameters for each candidate is also reduced and the significance of a real signal increased. This latter effect is due both to search being intrinsically more sensitive and to the trials' factor de-arXiv:1608.08928v1 [astro-ph.IM] 31 Aug 2016 creasing for every search from one stage to the next.
Each stage performs a stack-slide type of search using the GCT method and implementation of [6,7]. Important variables are: the coherent time baseline of the segments, the number of segments used (N seg ), the total time spanned by the data, the grids in parameter space and the detection statistic used to rank the parameter space cells. All stages use the same data set. The first three follow-up searches are performed on the Ein-stein@Home volunteer computing platform [8], the last on the Atlas computing cluster [9].
FIG. 1. These are the mismatch histograms of the four followup searches, so the y-axis represents normalised counts. For a given search and search set-up, the mismatch distribution depends of the template grid. The injection-and-recovery Monte Carlo studies to determine these distributions were performed without noise.
The parameters for the various stages are summarised in Table I. The grids in frequency and spindown are each described by a single parameter, the grid spacing, which is constant over the search range. The same frequency grid spacings are used for the coherent searches over the segments and for the incoherent summing. The spindown spacing for the incoherent summing step is finer than that used in for the coherent searches by a factor γ. The notation used here is consistent with that used in previous observational papers [1,4,12] and in the GCT methods papers [6,7].
The sky grids for stages 1 to 4 are approximately uniform on the celestial sphere projected on the ecliptic plane. The tiling is an hexagonal covering of the unit circle with hexagons' edge length d: with τ E 0.021 s being half of the light travel time across the Earth and m sky the so-called mismatch parameter. As was done in previous searches [1,2] the sky-grids are constant over 10Hz bands and the spacings are the ones associated through Eq. 1 to the highest frequency in the range. The sky grid of stage 0 is the union of two grids: one is uniform on the celestial sphere after projection onto the equatorial plane and the tiling (in the equatorial plane) is approximately square with edge d(0.3) from Eq.1; the other grid is limited to the equatorial region (0 ≤ α ≤ 2π and −0.5 ≤ δ ≤ 0.5), with constant actual α and δ spacings equal to d(0.3) -see Fig.1 of [1]. The reason for the equatorial "patching" with a denser sky grid is to improve the sensitivity of the search.
After each stage a threshold is set on the detection statistic to determine what candidates will be searched by the next stage. We set this detection threshold to be the highest such that the weakest signal that survived the first stage of the pipeline would, with high confidence, not be lost.
The set-up for each stage is determined at fixed computation cost. The computational cost is mostly set by practical considerations such as the time-frame on which we'd like to have a result, the number of stages that we envision in the hierarchy and the availability of Ein-stein@Home.
Since an analytical model that predicts the sensitivity of a search with the current implementation of the GCT method does not exist, we consider different search setups and for every set-up we perform fake-signal injection and recovery Monte Carlos. From these we determine the detection efficiency and the signal parameter uncertainty for signals at the detection threshold. We pick the search set-up based on these. Typically the search setup with the lowest parameter uncertainty volume (R given below) also has the highest detection efficiency and we pick that. As a further cross-check we also determine the mismatch distributions for the detection statistic. We define the mismatch µ as where F signal is the value of the detection statistic that we measure when we search the data with a template that is perfectly matched to the signal and F candidate is the value of the detection statistic that we obtain when running a search on a set of templates none of which, in general, will perfectly coincide with the signal waveform. The mismatch is hence a measure of how fine the grid that we are using is. As expected, Fig. 1 shows that the grids of subsequent stages get finer and finer. At each stage we determine the signal parameter uncertainty for signals at least at the detection threshold, in each search dimension: the distance in parameter space around a candidate that with high confidence (at least 90%), includes the signal parameter values. The uncertainty region around each candidate associated with stage i is searched in stage i + 1. The uncertainty volume at stage i is smaller than the containment volume of stage i − 1.

III. THE S6 SEARCH FOLLOW-UP
A series of all-sky Einstein@Home searches looked for signals with frequencies from 50 Hz through 510 Hz and frequency derivatives from 3.1 × 10 −10 Hz/s through −2.6 × 10 −9 Hz/s. Results from these were combined and analysed as described in [1]: no significant candidate was found and upper limits were set on the gravitational wave signal amplitude in the target signal parameter space. The data set that we begin with, is that described in Section III.1 and III.2 of [1]: a ranked-list of 3.8 × 10 10 candidates each with an associated detection statistic value 2F. We now take the 16 million most promising regions in parameter space from that search and inspect them more closely. This is done in four stages which we describe in the next subsections.
We remind the reader that some of the input data to this search was treated by substituting the original frequency-domain data with fake Gaussian noise at the same level as that of the neighbouring frequencies. This is done in frequency regions affected by well-known artefacts, as described in [1]. Results stemming entirely from this fake data are not considered in any further stage. Moreover, after the initial Einstein@Home search, the results in 50 mHz bands were visually inspected and those 50 mHz bands that present obvious noise disturbances are also removed from the analysis. A complete list of the excluded bands is given in the Appendixes of [1]. We will come back to this point as we present the results of this search.

A. Stage 0
This is the most complex stage of the hierarchy and determines the sensitivity of the search: if a signal does not pass this initial stage it will be lost forever. So we try here to keep the threshold that candidates have to exceed to be considered further as low as possible, compatibly with the feasibility of the next stage with the available computing resources. Such threshold was set at 2F = 6.109.
The identification of correlated candidates saves compute cycles in the next steps of the search. As was done in [4] the clustering procedure aims at bundling together candidates that could be ascribed to the same signal. In fact, a loud signal as well as a loud disturbance would produce high values of the detection statistic at a number of different template grid points and it would be a waste to follow up each of these independently. As described in [4,5] we begin with the loudest candidate, i.e. the candidate with the highest value of 2F . This is the seed for the first cluster. We associate with it close-by candidates in parameter space. Together, the seed and the nearby candidates, constitute the first cluster. We remove the candidates from the first cluster from the candidate list. The loudest candidate on the resulting list is the seed of the second cluster. We proceed in the same way as for the first cluster and reiterate the procedure until no more seeds with 2F values equal to or larger than 6.109 remain.
Monte Carlo studies are conducted to determine the cluster box size, i.e. the neighbourhood of the seed that determines the cluster occupants. We inject signals in Gaussian noise data at the level of our detectors' noise, search a small parameter space region around the signal parameters and use the resulting candidates as representative of what we would find in an actual search, after comparing them with the candidates obtained in the case of 'noise-only'. For signals at the detection threshold the 90% confidence cluster box is: If we consider as cluster occupants only those with 2F values greater or equal to 5.9, we observe that signals tend to produce slight over-densities in the clusters with respect to noise. This feature is exploited with an occupancy veto that discards all clusters with less than 2 occupants. We find that the false dismissal for signals at threshold is hardly affected (∼ 0.02% of signal clusters ) whereas the noise rejection is quite significant: we exclude 45% of signal clusters .
This same set of injection-data is utilised to characterise the false dismissals and the parameter uncertainty regions for all the stages of the hierarchy.
To summarise: the total number of candidates returned by the Einstein@Home searches is 3.8 × 10 10 . Of these we consider the ones with 2F above 6.109 , excluding frequency bands with obvious noise disturbances. There are 21.6 million such candidates. After clustering and occupancy veto we have reduced this number to 16.23 ×10 6 . The distribution of the detection statistic values 2F for these candidates is shown in Fig. 2 as well as their distribution in frequency. The maximum value is 8.6 and occurs at ∼ 52 Hz. All remaining values are smaller than 7.1 .

B. Stage 1
In this stage we search a volume of parameter space (Eqs. 3) around each candidate (around each cluster seed) equal to the cluster box defined in 3. We fix the total run time to be 4 months on Einstein@Home and this yields an optimal search set-up having the same coherent time baseline as stage 0, 60 hours, with the same number of segments N seg = 90 and the grid spacings shown in Table I. We use the same ranking statistic as in the original search [1], the O SGL [10], with the same tunings. The 90% uncertainty regions for this search set-up for signals just above the detection threshold are The search is divided among 16.23 ×10 6 work-units (WUs) each lasting about 2 hours and performed by one of the Einstein@Homevolunteer computers. From each follow-up search we record the most significant candidate. The distribution of these is shown in Fig.3. A threshold at F = 6.109 has a ∼ 9% false dismissal for signals at threshold and a 70% noise rejection. Using this threshold to determine what candidates to consider in the next stage yields 5.3 ×10 6 candidates.

C. Stage 2
In this stage we search a volume of parameter space around each candidate defined by Eq. 4. As shown in Table I As done in stage 1 we record the most significant candidate from each search. The distribution is shown in Fig.5. In the next stage we follow-up the top 1.1 million candidates, corresponding to a threshold on 2F at 7.38 . This threshold has a ∼ 0.6% false dismissal for signals at threshold and a 79% noise rejection. In this stage we search a volume of parameter space around each candidate defined by Eq.5. As shown in Table I, the coherent timebaseline is as long as that used in the previous stage but the grid spacings are finer. The search is divided among 1.1 million WUs each lasting about 2 hours.
The > 99% uncertainty regions for this search set-up for signals close to the detection threshold are As done in previous stages we record the most significant candidate from each search. The distribution is shown in Fig.7. In the next stage we follow-up the top 10 candidates, corresponding to a threshold on 2F at 8.82 . This threshold has a ∼ 4×10 −4 false dismissal for signals at threshold and a 99.9991% noise rejection.

E. Stage 4
In this stage we search a volume of parameter space around each candidate defined by Eq.7. The set-up of choice has a coherent time-baseline of 280 hrs, twice as long as that used in stage-3, and the grid spacings shown in Table I. The search has a relatively modest cost and is performed on the Atlas cluster: each follow-up lasts about 14 hours. The ranking statistic is O SGL with a retuned c * = 96.1. We consider the loudest candidate from each of the 10 follow-ups. In our Monte Carlo studies None of those injections has a 2F below 16.2. Conservatively, we pick a threshold at 15.0 . The Gaussian false alarm at 2F = 15.0 is very low (≈ 2 × 10 −20 ) and hence we do not expect any candidate from random Gaussian noise fluctuations. Since we only follow-up 10 candidates we report our findings explicitly for each of them in Tab.II. Candidates 3, 4 and 6 satisfy this condition but unfortunately they are ascribable to fake signals hardwareinjected in the detector to test the detection pipelines. The search recovers all fake signals in the data with parameters within its search range and not absurdly loud 1 . We note that candidates 3 and 4 come from the same fake signal. For a complete list of the fake signals present in the data -see Table 6 of [11]. In Table III we show the signal parameters and report the distance with respect to the candidate parameter values. These distances are all 1 A fake signal was injected at about 108 Hz at such a high amplitude that it saturates the Einstein@Home toplists across the entire sky. Upon visual inspection it is immediately obvious that the f −ḟ morphology is that of a signal, albeit an unrealistically loud one. We categorised the associated band as disturbed because the data is corrupted by this loud injection and it is impossible to detect any real signal in its frequency neighbourhood.
within the Stage-4 uncertainties of 7. We do not follow up these candidates any further because we know that they are associated with the hardware injections. The remaining candidates are below the threshold of 15.0 which is the minimum value of 2F that we demand candidates to pass before we inspect them further.
FIG. 9. Fraction of signals that are recovered with a detection statistic value larger than or equal to the threshold value (vertical line) after the stage-4 follow-up. The dashed line is a linear extrapolation based on the last two data points to guide the eye to the false dismissal value for signals at threshold. This line is a conservative estimate in the sense that it over estimates the false dismissal.

IV. RESULTS
The search did not reveal any continuous gravitational wave signal in the parameter volume that was searched. We hence set frequentist upper limits on the maximum gravitational wave amplitude consistent with this null result in 0.5 Hz bands : h 90% 0 (f). h 90% 0 (f) is the GW amplitude such that 90% of a population of signals with parameter values in our search range would have been detected by our search, i.e. would have survived the last 2F threshold at 15.0 at stage-4. Since an actual full scale injection-and-recovery Monte Carlo for the entire set of follow-ups in every 0.5 Hz band is prohibitive, in the same spirit as [1,12], we perform such study in a limited set of trial bands. We pick 100. For each of these we determine the sensitivity depth of the search corresponding to the detection criterion stated above. As representative of the sensitivity depth D 90% of this hierarchical search, we take the average of these depths, 46.9 [1/ √ Hz]. Given the noise level of the data as a function of frequency, S h (f ) , we then determine the 90% upper limits as  III. Columns 2-5 show the parameters of the fake injected signal closest to the candidate whose ID identifies it in Table  II. Columns 6-9 display the distance between the candidates' and the signals' parameters (candidate parameter minus signal parameter). The reference time eference time (GPS s): 960541454.5 Figure 10 shows these upper limits as a function of frequency. They are also presented in tabular form in the Appendix with the associated uncertainties which amount to 20%, including calibration uncertainties. The most constraining upper limit is in the band between 170.5 and 171 Hz and it is 4.3 ×10 −25 . At the upper end of the frequency range, around 510 Hz, the upper limit rises to 7.6 ×10 −25 .
The upper limits can be recast as exclusion regions in the signal frequency-ellipticity plane parameterised by the distance, for an isolated source emitting continuous gravitational waves due to its shape presenting an ellipticity where I are the principal moments of inertia and the coordinate system is taken so that the z-axis is aligned with the spin axis of the star. Fig.11 shows these upper limits. Above 200 Hz we can exclude sources with ellipticities larger than 10 −6 within 100 pc of Earth and above 400 Hz ellipticities above 4 × 10 −7 .

V. CONCLUSIONS
With a hierarchy of five semi-coherent searches at increasing coherent time baselines and resolutions in parameter space we searched over 16 million regions in few hundred Hz around the most sensitive frequencies of the LIGO detectors during the S6 science run. All stages but the very last ran on the Einstein@Home distributed computing project lasting a few to several weeks. This is the first large-scale hierarchical search for gravitational wave signals ever performed.
Having carried out this search proves that one can successfully perform deep follow-ups of marginal candidates and elevate their significance to the level necessary to be able to claim a detection. This paper proves that searches with thresholds at the level of the Einstein@Home search described in [16] are possible; [16] demonstrates that they are the most sensitive and our observational results confirm this.
The sensitivity of broad surveys for continuous gravitational wave signals is computationally limited. For this reason we employ Einstein@Home to deploy our searches. However, following up tens of millions of candidates is not just a matter of having the computational power. This paper illustrates how to perform and optimise the different stages, factoring in all the practical aspects of a real analysis.
None of the investigated candidates survived the five stages apart from those arising from the two fake signals injected in the detector for control purposes. These fake signals were recovered with the correct signal parameters. Candidate 6 comes from a hardware injection weak enough that no other search on this data set was ever able to detect it. This search recovers it well above the detection threshold.
The gravitational wave amplitude upper limits that we set improve on existing ones [1] by about 30%. This corresponds to an increase in accessible space volume of 2.
We excluded 10% of the original data from this analysis because the stage-0 results had different statistical properties than the bulk of the results and the automated frequency Hough search results [14] Powerflux search results [13] results from Stage-0 [1] results from this search  [1]. These lie on the curve above the lowest one and are marked by dark blue diamonds. The results from a previous broad all-sky survey [13] are the top curve (lighter circles and crosses) above 100 Hz. In the lower frequency we compare with a search on Virgo data contemporary to the LIGO S6 data [14].
methods employed here which are necessary in order to deal with a large number of candidates, would not have yielded meaningful statistical results. We might go back to these excluded parameter space regions and attempt to extract information. This is a time-consuming process and the odds of finding a signal versus the odds of missing one by not analysing more sensitive data might well indicate that we shouldn't pursue this.
The optimal set-up for the various stages and the upper limits were determined at the expense of signal injectionand-recovery Monte Carlo studies. This is due to the fact that the implementation of a stack-slide search that we are using does not allow an analytical prediction of the sensitivity of a search with a given set-up (coherent segments and grid spacings). This major drawback will soon be overcome by a new implementation of stack-slide searches based on [17][18][19][20]. Such a search is being characterised and tuned at the time of writing and we hope to employ it in the context of our contributions to the LIGO Scientific Collaboration for searches on data from the next LIGO run (O2).
In principle we would like to carry out the entire hierarchy of stages on Einstein@Home. For this to happen two aspects of the search presented here need to be automated: the visual inspection and the follow-up stages. The first is underway [21]. The second will be significantly eased by the new stack-slide search that we alluded to, above.

VI. ACKNOWLEDGMENTS
Our sincere gratitude firstly goes to the Einstein@Home volunteers who have made these searches possible: thank you, thank you, thank you. Maria Alessandra Papa, Sinéad Walsh, Bruce Allen and Xavier Siemens gratefully acknowledge the support from NSF PHY Grant 1104902. All the tuning and preparatory computational work for this search was carried out on the ATLAS super-computing cluster at the Max-Planck-Institut für Gravitationsphysik/ Leibniz Universität Hannover. We also acknowledge the Contin-  11. Ellipticity of a source at a distance d emitting continuous gravitational waves that would have been detected by this search. The dashed line shows the spin-down ellipticity for the highest magnitude spindown parameter value searched: 2.6 ×10 −9 Hz/s. The spin-down ellipticity is the ellipticity necessary for all the lost rotational kinetic energy to be emitted in gravitational waves. If we assume that the observed spin-down is all actual spin-down of the object, then no ellipticities could be possible above the dashed curve. In reality the observed and actual spindown could differ due to radial motion of the source. In this case the actual spin-down of the object may even be larger than the apparent one. In this case our search would be sensitive to objects with ellipticities above the dashed line.
uous Wave Group of the LIGO Scientific Collaboration for useful discussions and Graham Woan for reviewing the manuscript on behalf of the Collaboration. This document has been assigned LIGO Laboratory document number LIGO-P1600213.