BEATING THE SPIN-DOWN LIMIT ON GRAVITATIONAL WAVE EMISSION FROM THE VELA PULSAR

We describe here a search for continuous gravitational radiation from the Vela pulsar (PSR B0833−45, PSR J0835−4510) in data from the Virgo detector VSR2 run, which began on 2009 July 7 and ended on 2010 January 8. Continuous gravitational waves (CWs) can be emitted by a rotating neutron star through a variety of possible mechanisms, including non-axisymmetry of its mass distribution, giving rise to a time-varying quadrupole moment. Such emission would imply loss of rotational energy and decrease in spin frequency. Hence, a pulsar's observed frequency spin-down can be used to place an indirect upper limit on its gravitational wave (GW) emission, named spin-down limit. While a recent search for CW radiation using LIGO data has been carried out for more than 100 known pulsars (Abbott et al. 2010), the resulting upper limits have beaten the spin-down limit for only the Crab pulsar (Abbott et al. 2008, 2010). A search over LIGO data for CW signals from the non-pulsing neutron star in the supernova remnant Cassiopeia A has established an upper limit on the signal amplitude over a wide range of frequencies which is below the indirect limit derived from energy conservation (Abadie et al. 2010). In this paper we present upper limits on CW emission from the Vela pulsar that lie below its spin-down limit, making Vela only the second pulsar for which this experimental milestone has been achieved. The only previous targeted search for CW emission from the Vela pulsar was in CLIO data over the period 2007 February 12–28, which produced an upper limit of ∼5.3 × 10⁻²⁰, several orders of magnitude above the spin-down limit (Akutsu et al. 2008).

Vela is observed to pulsate (f_rot ≃ 11.19 Hz) in radio, optical, X-ray, and γ-ray radiation and is associated with the Vela supernova remnant. The association of the pulsar to the supernova remnant was made in 1968 (Large et al. 1968) and was the first direct observational proof that supernovae can produce neutron stars. The Vela spin-down rate is $\dot{f}_{{\rm rot}}\simeq -1.56\!\times \!10^{-11}$ Hz s⁻¹, corresponding to a kinetic energy loss of $\dot{E}_{{\rm sd}}\simeq 6.9\!\times \!10^{29}$ W, where the canonical value for a neutron star's moment of inertia, I = 10³⁸ kg m², has been assumed. This loss of energy is due to various mechanisms, including magnetic dipole radiation, acceleration of charged particles in the pulsar magnetosphere and possibly the emission of GWs. In this analysis we assume a triaxial neutron star rotating around a principal axis of inertia, so that the GW signal frequency is f = 2f_rot (see Section 2). With an estimated distance from the Earth of ∼290 pc (Dodson et al. 2003), Vela is one of the nearest known pulsars. Assuming that all the observed spin-down is due to the emission of GWs, we obtain the spin-down limit h^sd₀ = 3.29  ×  10⁻²⁴ for GW tensor amplitude at the Earth. With an estimated age of ∼11,000 yr (Caraveo & Bignami 1989), Vela is relatively young and could, in principle, have a significant residual non-axisymmetry from its formation. The spin-down limit on the signal amplitude can be converted into an upper limit on the star's equatorial ellipticity (see Equation (15)). For Vela we have ^sd = 1.8 × 10⁻³. This value is far larger than the maximum allowed by standard equations of state for neutron star matter (Horowitz & Kadau 2009), but is comparable to the maximum value foreseen by some exotic equations of state (Owen 2005; Lin 2007; Haskell et al. 2007). Because of very effective seismic isolation (Acernese et al. 2010), Vela's GW emission frequency (f ≃ 22.38 Hz) is within the sensitive band of the Virgo detector; this frequency range is inaccessible to all other GW detectors to date.

Vela is a particularly glitchy pulsar, with an average glitch rate of ∼1/3 yr⁻¹, making it important to know whether or not a glitch occurred during the VSR2 run. Vela is regularly monitored by both the Hobart radio telescope in Tasmania and the Hartebeesthoek radio telescope in South Africa. According to their observations, no glitch occurred during the time span of VSR2. Prior to VSR2 it last glitched on 2007 August 1, and it has since glitched on 2010 July 31 (Buchner 2010). Observations from the Hobart and Hartebeesthoek telescopes have also been used to produce updated ephemerides for Vela, which are important given Vela's relatively large timing noise. If timing noise is a consequence of fluctuations in the star's rotation frequency, not taking it into account would result in an increasing mismatch over time between the signal and template phases, thus producing a sensitivity loss in a coherent search. In this search updated ephemerides have been computed using the pulsar software TEMPO2 starting from the set of times of arrival of the electromagnetic pulses observed by the Hobart and Hartebeesthoek telescopes covering the whole duration of the VSR2 run. Including in the fitting process up to the second derivative of frequency is enough in order to have flat post-fit residuals. The post-fit position and frequency parameters are shown in Tables 1 and 2, respectively. The corresponding post-fit residuals rms amounts to a negligible 100 μs.

Table 1. Position and Estimated Distance of the Vela Pulsar

α	δ	d (pc)
$08^{{\rm h}} 35^{{\rm m}} 20\mbox{$.\!\!^{\mathrm s}$}75438(3)$	−45°10'329507(7)	287 (− 17, +19)

Notes. For the position, parentheses give the 1σ error on the final digit as produced by the TEMPO2 fit; for the distance, the uncertainty estimated in Dodson et al. (2003) is quoted. Positional parameters refer to epoch (MJD) 54620.

Download table as:  ASCIITypeset image

Table 2. Spin Frequency, Spin-down Rate, and Estimated Age of the Vela Pulsar

frot (Hz)	$\dot{f}_{{\rm rot}}$ (Hz s−1)	$\ddot{f}_{{\rm rot}}$ (Hz s−2)	Age (yr)
11.19057302331(9)	−1.5583876(4) × 10−11	4.9069(9) × 10−22	11000

Notes. Parentheses give the 1σ error on the final digit of spin frequency and spin-down rate estimations as produced by the TEMPO2 fit. Rotational parameters refer to epoch (MJD) 54620. The quoted precision is enough to determine the rotational phase to within about 0.012 cycles.

Download table as:  ASCIITypeset image

Recent Chandra X-ray observations provide accurate determination of the orientation of the Vela spin axis. In Ng & Romani (2008), estimates of the pulsar wind nebula's "position angle," ψ_P, and inclination ι_P are given:

$$\begin{eqnarray} \nonumber \psi _P & = & 130\mbox{$.\!\!^\circ $}63 \pm 0\mbox{$.\!\!^\circ $}05, \\ \iota _P & = & 63\mbox{$.\!\!^\circ $}6 \pm 0\mbox{$.\!\!^\circ $}6. \end{eqnarray} \tag{ 1 }$$

The "position angle" is related to the GW polarization angle ψ (see Section 2) by either ψ = 180° + ψ_P or ψ = ψ_P depending on the unknown spin direction. Our analyses are insensitive to rotations of ψ by integer multiples of 90°, so the spin direction is not needed. The inclination angle calculated from the pulsar wind nebula ι_P is taken to be the same as that of the pulsar ι. The physics of pulsar wind nebulae is complex, and a model leading to the above fits has several uncertainties. Thus, we perform separate searches for the GW signal from Vela, both assuming that the angles ψ and ι are known within the above uncertainties, and assuming that they are unknown.

The remainder of this paper is organized as follows. In Section 2 we summarize the characteristics of the GW signals for which we search. In Section 3 we describe the data set used for the analysis. In Section 4 we briefly describe the three analysis methods used. In Section 5 we present the results of the analysis. In Section 6 we provide conclusions. Some more details on the analysis methods are given in the Appendices.

The continuous GW signal emitted by a triaxial neutron star rotating around a principal axis of inertia as seen from Earth is described by the following tensor metric perturbation:

$$\begin{equation} {\bf h}(t)=h_+(t)\, {\bf e_+} + h_{\times }(t)\, {\bf e_{\times }}, \end{equation} \tag{ 2 }$$

where

$$\begin{eqnarray} h_+(t) &=& h_0\left(\frac{1+\cos ^2{\iota }}{2}\right)\cos {\Phi (t)} \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} h_{\times }(t) &=& h_0\cos {\iota }\sin {\Phi (t)}, \end{eqnarray} \tag{ 4 }$$

and e₊ and e_× are the two basis polarization tensors. They are defined, see, e.g., Misner et al. (1973), in terms of unit orthogonal vectors e_x and e_y where e_x is along the x-axis of the wave frame, defined as the cross product $\hat{s}\times \hat{n}$ between the source spin direction $\hat{s}$ and the source direction $\hat{n}$ in the solar system barycenter (SSB).

The angle ι is the inclination of the star's rotation axis with respect to the line of sight and Φ(t) is the signal phase function, where t is the detector time, while the amplitude h₀ is given by

$$\begin{equation} h_0=\frac{4\pi ^2G}{c^4}\frac{I_{zz}\epsilon f^2}{d}, \end{equation} \tag{ 5 }$$

where I_zz is the star moment of inertia with respect to the rotation axis, the equatorial ellipticity is defined, in terms of principal moments of inertia, as $\epsilon =\frac{I_{xx}-I_{yy}}{I_{zz}}$, d is the star distance, and f is the signal frequency. As the time-varying components of the mass quadrupole moment tensor are periodic with period half the star rotation period, it follows that f = 2f_rot.

The GW strain at the detector can be described as

$$\begin{equation} h(t)=h_+(t)F_+(t;\psi)+h_{\times }(t)F_{\times }(t;\psi), \end{equation} \tag{ 6 }$$

where the two beam-pattern functions, which are periodic functions of time with period of one sidereal day, are given by

$$\begin{eqnarray} F_+(t;\psi) &=& a(t)\cos {2\psi }+b(t)\sin {2\psi } \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} F_{\times }(t;\psi) &=& b(t)\cos {2\psi }-a(t)\sin {2\psi }. \end{eqnarray} \tag{ 8 }$$

The two functions a(t), b(t) depend on the source position in the sky and on the detector position and orientation on the Earth. Their time dependency is sinusoidal and cosinusoidal with arguments Ω_⊕ t and 2 Ω_⊕ t, where Ω_⊕ is the Earth angular rotation frequency; ψ is the wave polarization angle defined as the angle from $\hat{z} \times \hat{n}$ to the x-axis of the wave frame, measured counterclockwise with respect to $\hat{n}$, where $\hat{z}$ is the direction of the north celestial pole (see, e.g., the plot in Prix & Krishnan 2009). The effect of detector response on a monochromatic signal with angular frequency ω₀ is to introduce an amplitude and phase modulation which determine a split of the signal power into five frequencies, ω₀, ω₀ ± Ω_⊕, ω₀ ± 2Ω_⊕. The distribution of power among the five bands depends on the source and detector angular parameters. In Figure 1 the power spectrum at the Virgo detector of a hypothetical monochromatic signal coming from the location of the Vela pulsar is shown for two assumed polarizations (pure "+" linear polarization and circular left-handed polarization).

Zoom In Zoom Out Reset image size

To a very good approximation, the SSB can be used as an inertial reference frame in which to define the signal phase. In this frame, with barycentric time T, the signal phase is

$$\begin{equation} \Phi (T)=\Phi _0+2\pi f_0 \left(T-T_0\right), \end{equation} \tag{ 9 }$$

where the signal intrinsic frequency f₀ is a function of time due to the spin-down:

$$\begin{equation} f_0(T)=f^{(0)}+\sum _{n=1}^{2}\frac{f^{(n)}}{n!}\left(T-T_0\right)^n, \end{equation} \tag{ 10 }$$

where $f^{(n)}=\frac{d^nf_0}{dT^n}|_{T=T_0}$. The time at the detector, t, differs from T due to the relative motion between the source and the detector and to some relativistic effects. Considering only isolated neutron stars, we have the well-known relation (Lyne & Graham-Smith 1998; Hobbs et al. 2006; Edwards et al. 2006)

$$\begin{equation} T=t+\Delta _{R}+\Delta _E+\Delta _S, \end{equation} \tag{ 11 }$$

where

$$\begin{equation} \Delta _{R}=\frac{\vec{r}\cdot \hat{n}}{c} \end{equation} \tag{ 12 }$$

is the classical Roemer delay, which gives the main contribution ($\vec{r}$ is the vector identifying the detector position in the SSB, while $\hat{n}$ is the unit vector toward the source). The term Δ_E is the Einstein delay which is the sum of two contributions, one due to the gravitational redshift produced by the Sun and the other due to the time dilation produced by Earth's motion. Δ_S is the Shapiro delay due to the curvature of spacetime near the Sun. Expressing the signal phase in the detector frame, by using Equation (11), we can write the signal frequency at the detector as

$$\begin{equation} f(t)=\frac{1}{2\pi }\frac{d\Phi (t)}{dt}\simeq f_0(t) \left(1+\frac{\vec{v} \cdot \hat{n}}{c}\right)+{\rm rel.\ corr.}, \end{equation} \tag{ 13 }$$

where $\vec{v}$ is the detector velocity vector and terms of order $|f^{(1)}\frac{\vec{r}\cdot \hat{n}}{c}|$ or smaller have been omitted from the equation (though they are included in the analyses).

A useful quantity to which to compare the upper limit on signal strength set in a given analysis is the so-called spin-down limit. It is computed (Abbott et al. 2007) assuming that all the observed spin-down is due to the emission of GWs:

$$\begin{equation} h^{{\rm sd}}_0= 8.06\!\times \!10^{-19}\,I_{38}\,d^{-1}_{\rm kpc}\,\sqrt{\frac{|(\dot{f}_{{\rm rot}}/{\rm Hz}\,{\rm s}^{-1})|}{(f_{{\rm rot}}/{\rm Hz})}}, \end{equation} \tag{ 14 }$$

where I₃₈ is the star's moment of inertia in units of 10³⁸ kg m² and d_kpc is the star's distance from the Sun in kiloparsecs. It is an absolute upper limit to the amplitude of the GW signal that could be emitted by the star, where electromagnetic radiation is neglected. The spin-down limit on the signal amplitude corresponds to an upper limit on the star's ellipticity given by

$$\begin{equation} \epsilon ^{{\rm sd}}=0.237 \left(\frac{h^{{\rm sd}}_{0}}{10^{-24}}\right)I^{-1}_{38}(f_{{\rm rot}}/{\rm Hz})^{-2}\,d_{\rm kpc}. \end{equation} \tag{ 15 }$$

The Vela pulsar has a measured braking index n ≃ 1.4 (Lyne et al. 1996) and this, together with the estimation of its age, can be used to compute a stricter indirect limit on the signal amplitude (Palomba 2000), which only holds under the assumption that the spin-down is due to the combination of emission of GW and magnetic dipole radiation, about four times lower than the spin-down limit.

Achieving sensitivity better than the spin-down limit is an important milestone toward probing neutron star structure via GWs.

We have analyzed calibrated strain data from the Virgo VSR2 run. This run (started in coincidence with the start of the LIGO S6 data run) began on 2009 July 7 21:00:00 UTC (GPS 931035615) and ended on 2010 January 8 22:00:01 UTC (GPS 947023216). The duty cycle was 80.4%, resulting in a total of ∼149 days of science mode data, divided among 379 segments. Science mode is a flag used to indicate when the interferometer is locked and freely running at its working point, with all the controls active and no human intervention. In Figure 2, the fraction of total time covered by science data segments with duration longer than a given value is plotted. The longest segment lasts ∼88 hr.

Zoom In Zoom Out Reset image size

The detector showed a good sensitivity around the expected Vela signal frequency during the entire run. The sensitivity was typically within a factor of two of the target Virgo design sensitivity (Accadia et al. 2010). Figure 3 shows the estimation of the power spectrum of the data, computed through an average of ∼1000 s periodograms after removal of some large outliers (see Section 4.3.1), on a 0.8 Hz frequency band around the expected frequency of the GW signal from the Vela pulsar for the entire VSR2 run. An instrumental disturbance right at the Vela signal frequency degraded the sensitivity by ∼20% with respect to the background. The source of this disturbance was seismic noise produced by the engine of the chiller pumps that circulate coolant fluid for the laser of the mirror thermal compensation system and it has been removed during the next Virgo VSR3 run.

Zoom In Zoom Out Reset image size

The data used in the analysis have been produced using the most up-to-date calibration parameters and reconstruction procedure. The associated systematic error amounts to 5.5% in amplitude and ∼50 mrad in phase (Accadia et al. 2010) over the frequency range between ∼10 Hz and ∼1 kHz, with lower uncertainties at the Vela frequency. The reconstructed data have a sampling rate of 20000 Hz. However, two more reconstructed data streams, sampled respectively at 16384 Hz and 4096 Hz, were also produced to be consistent with LIGO/GEO sampling rates.

Three different and largely independent analysis methods have been applied to this search: (1) a complex heterodyne method using Bayesian formalism and a Markov Chain Monte Carlo (MCMC; Abbott et al. 2010), (2) a time-domain matched filter method using the ${\mathcal F}$ statistic (Jaranowski et al. 1998) and a new extension known as the ${\mathcal G}$ statistic (Jaranowski & Królak 2010), and (3) a matched filter method applied to the signal's Fourier components at the five frequencies to which the signal is spread by the sidereal modulation (Astone et al. 2010).

There are several reasons to use different methods in the search for CW signals, provided they have comparable performance. First, it makes it easier to cross-check each method by comparing the analysis outputs, even at intermediate steps. Second, different methods can be more suitable, or efficient, for given characteristics of the data to be analyzed, or for given characteristics of the signal emitted by a source, e.g., a method can be more robust against noise non-stationarity with respect to another. Third, in case of detection with a given analysis it will be of paramount importance to confirm the detection with one or more independent analyses.

In the analyses described in this paper, we observe consistent results from the three methods, which provide valuable cross-checks.

All the analyses clean the data in some way to remove large transient outliers. This is necessary, as large short-duration transients will skew noise estimates and adversely affect results. The amount of data removed during cleaning is negligible compared to the total data span and would produce a decrease of the signal-to-noise ratio (S/N) of a signal present in the data of less than 1%.

Among the three methods, two different approaches have been used toward setting upper limits. In the heterodyne method, the posterior probability for the signal parameters is calculated, from which degree-of-belief (or credibility) regions can be set to give limits on particular parameters (e.g., an upper limit on h₀ can be set by finding the value that bounds a given percentage of the probability). In the two other analyses, a frequentist approach is used and upper limits are set through Monte Carlo methods where many simulated signals with different amplitude and randomly varying parameters and frequency near the expected one from the Vela are added to the data. These two approaches should produce quantitatively similar results, see, for example, Abbott et al. (2004), but they are answering different questions and therefore cannot be meaningfully combined.

The three analysis methods are described in the following sections of this paper.

4.1. Complex Heterodyne

This method, developed in Dupuis & Woan (2005), provides a way to reduce the search data set to a manageable size, and use it to perform Bayesian parameter estimation of the unknown signal parameters.

4.1.1. Data Reduction

The known signal phase evolution (Equation (9)) is used to heterodyne the data, changing the time series detector data x(t) = h(t) + n(t), where h(t) is the signal given by Equation (6) and n(t) is the noise, to

$$\begin{equation} x^{\prime }(t) = x(t)e^{-i[\Phi (t) - \Phi _0]}, \end{equation} \tag{ 16 }$$

giving a complex data set in which the signal is given by $h^{\prime }(t) = h(t)e^{-i[\Phi (t) - \Phi _0]}$. The now complex signal is

$$\begin{eqnarray} h^{\prime }(t) & = & h_0\left(\frac{1}{4}F_+(1+\cos {}^2\iota)\cos {\Phi _0} + \frac{1}{2}F_{\times }\cos {\iota }\sin {\Phi _0} \right)\nonumber \\ & &+ i h_0 \left(\frac{1}{4}F_+(1+\cos {}^2\iota)\sin {\Phi _0} - \frac{1}{2}F_{\times }\cos {\iota }\cos {\Phi _0} \right),\nonumber\\ \end{eqnarray} \tag{ 17 }$$

where F₊ and F_× are given by Equations (7) and (8). This heterodyne therefore removes the fast-varying part of the signal (the time dependent part of Equation (9)) leaving a complex data stream with the signal shifted to zero frequency (setting aside small offsets due to the diurnal amplitude modulation of the signal from the detector beam pattern). In practice this heterodyne is performed in a two-stage process. First, a coarse heterodyne is performed using the phase evolution calculated assuming a stationary frame. These data are then low-pass filtered (in this case using a ninth-order Butterworth filter with a 0.25 Hz knee frequency) and heavily downsampled from the original rate of 16384 Hz to 1 Hz. A second stage of heterodyne takes into account the signal's modulation due to Earth's motion and relativistic effects (see Equation (13)). The data are then further downsampled from 1 Hz to 1/60 Hz by taking the mean of 60 samples, which has the effect of an additional low-pass filter.

4.1.2. Data Cleaning

The fully heterodyned data are cleaned to remove the largest outliers, by discarding points with absolute values greater than five times the standard deviation of the data. This cleaning is performed twice to combat the effect of extreme outliers (many order of magnitude larger than normal) skewing the standard deviation estimate. This removes ∼0.05% of the data.

For the parameter estimation, as in Abbott et al. (2007, 2010), the likelihood calculation assumes the data are stationary for contiguous 30 minute segments, although shorter segments of 5 minutes or more are also included to account for shorter stretches of data at the end of longer contiguous segments. This contiguity requirement removes a further ∼0.2% of the heterodyned data, which is within segments shorter than 5 minutes.

4.1.3. Parameter Estimation and Upper Limits

This new, and far smaller, 1/60 Hz sampled data set is then used to estimate the four unknown signal parameters h₀, Φ₀, cos ι, and ψ. These are estimated using a Bayesian formalism, with a Students-t-like distribution for the likelihood (formed by marginalizing a Gaussian likelihood over an unknown noise standard deviation) given the heterodyned data and a signal model from Equation (17), and specific priors (see below) on these parameters. This posterior probability volume is explored using an MCMC (Abbott et al. 2010), which gives posterior probability distribution functions (PDFs) on each parameter marginalized over the three others.

In this analysis two different sets of independent priors are used for the parameters. In one case uniform priors on all four parameters are set—for the angular parameters this means that they are uniform across their allowable ranges, but for h₀ the lower bound is zero, and the upper bound is set at a level well above any values that could be consistent with the data. For reasons set out in Section 1, the other case sets the priors on ψ and cos ι to be Gaussians given by Equation (1), whilst keeping the h₀ and Φ₀ priors as uniform.

The marginalized h₀ posterior, p(h₀|d, I), can be used to set an upper limit in the amplitude by finding the value of h^ul₀ that bounds (from zero) the cumulative probability to a given degree of belief, B,

$$\begin{equation} B =\int _0^{h_0^{\rm ul}} p(h_0|d,I) dh_0. \end{equation} \tag{ 18 }$$

Here we set 95% degree-of-belief upper limits. Due to the fact that the MCMC is finite in length there will be small statistical uncertainties between different MCMC runs, which for cleaned data we find to be ≲ 1 × 10⁻²⁶. The difference in results between using cleaned and non-cleaned data, as above, is within the statistical uncertainty from the MCMC.

4.2. ${\mathcal F}$ and ${\mathcal G}$ Statistics Method

4.3. Matched Filter on the Signal Fourier Components

The third search method uses the Fourier amplitudes computed at five frequencies where the signal would appear due to sidereal amplitude modulation and applies a matched filter to this five-point complex data vector. Further details can be found in Appendix A and in Astone et al. (2010).

4.3.1. Data Reduction

The starting point for this method is a short Fourier transform database (SFDB) built from calibrated strain data sampled at 4096 Hz (Astone et al. 2005). The FFTs have a duration of 1024 s and are interlaced by 50% and windowed with a flat top—cosine edges window. From the SFDB a small band (0.2 Hz in this analysis) around the frequency of interest is extracted from each FFT. The SFDB contains, among other information, the position and the velocity of the detector in the SSB at the center time of each FFT. Each frequency domain chunk is zero-padded and inversely Fourier-transformed to obtain a complex time series with the same sampling time of the original time series, but with a spectrum different from zero only in the selected band (i.e., it is an analytical signal, see, e.g., Astone et al. 2002). Then, for each sample, the detector position in the SSB is computed, by interpolating with a third-degree polynomial. The Doppler and Einstein effects can be seen as a varying time delay Δ(t). A new non-uniformly sampled time variable t' with samples t'_i = t_i + Δ(t_i) is computed. The spin-down is corrected by multiplying each data chunk by $e^{-i \Delta \phi _{{\rm sd}}(t^{\prime })}$ where $\Delta \phi _{{\rm sd}}(t^{\prime })=2\pi (\dot{f}\frac{t^{\prime 2}}{2}+\ddot{f}\frac{t^{\prime 3}}{6})$. Then the data are resampled at equal intervals in t'. The final complex time series has a sampling frequency of 1 Hz. At this point, a true GW signal would be sinusoidal with a sidereally modulated amplitude and phase, as described in Section 2, containing power at the nominal source frequency and in lower and upper sidebands of ±Ω_⊕, ±2Ω_⊕. The Fourier coefficients at these five frequencies are taken to form a complex data 5-vector $\boldsymbol{X}$.

The detection method described here relies on a description of the GW signal given in Appendix A. When the polarization angle ψ and the inclination angle of the star rotation axis ι are unknown, we use a procedure that we denote 4 degrees of freedom (dof) detection, in which the two signal 5-vectors $\boldsymbol{A^+},\boldsymbol{A^{\times }}$, corresponding to the + and × polarizations and defined in Appendix A, are numerically computed and projected onto the data 5-vector $\boldsymbol{X}$:

$$\begin{eqnarray} &&\hat{H}_+=\frac{\boldsymbol{X}\cdot \boldsymbol{A^+}}{|\boldsymbol{A^+}|^2} \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} &&\hat{H}_{\times }=\frac{\boldsymbol{X}\cdot \boldsymbol{A^{\times }}}{|\boldsymbol{A^{\times }}|^2}. \end{eqnarray} \tag{ 20 }$$

The output of the two matched filters are the estimators of the amplitudes $H_0e^{i\Phi _0}H_+,\, H_0e^{i\Phi _0}H_{\times }$. The final detection statistic is defined by

$$\begin{equation} \textsl {S}=|\boldsymbol{A^+}|^4 |\hat{H}_+|^2+|\boldsymbol{A^{\times }}|^4 |\hat{H}_{\times }|^2. \end{equation} \tag{ 21 }$$

More details can be found in Astone et al. (2010).

If estimations of ψ and ι provided by X-ray observations (Section 1) are used, we can apply a simpler procedure that we call a 2 dof detection. In this case the signal is completely known, apart from an overall complex amplitude $H=H_0 e^{i\Phi _0}$. Then, the template consists of just one 5-vector $\boldsymbol{A}=H_+\boldsymbol{A^+}+H_{\times } \boldsymbol{A^{\times }}$, where H₊, H_× are given by Equation (A3), and only one matched filter must be applied to the data 5-vector $\boldsymbol{X}$:

$$\begin{equation} \hat{H}=\frac{\boldsymbol{X}\cdot \boldsymbol{A}}{|\boldsymbol{A}|^2}, \end{equation} \tag{ 22 }$$

which provides an estimation of the signal complex amplitude. The detection statistic is then given by $\textsl {S}=|\hat{H}|^2$.

4.3.2. Data Cleaning

In addition, various cleaning steps were applied to the data. The data can be modeled as a Gaussian process, with slowly varying variance, plus some unmodeled pulses affecting the tails of data distribution. The cleaning procedure consists of two parts. First, before the construction of the SFDB, high-frequency time-domain events are identified after applying to the data a first-order Butterworth high-pass bilateral filter, with a cutoff frequency of 100 Hz. These events are then subtracted from the original time series. In this way we do not reduce the observation time because we are simply removing from the data the high-frequency noisy component. The effect of this kind of cleaning has been studied in data from Virgo Commissioning and Weekly Science runs and typically reduces the overall noise level by up to 10%–15%, depending on the quality of the data (Acernese et al. 2009). After Doppler and spin-down correction, further outliers that appear in the small band to be analyzed are also removed from the data set by using a threshold of ±5  ×  10⁻²¹ on the data strain amplitude, reducing the amount of data by ∼1.3%. Slow non-stationarity of the noise is taken into account by applying a Wiener filter to the data, in which we estimate the variance of the Gaussian process over periods of ∼1000 s, and weight the data with its inverse in order to de-emphasize the more disturbed periods.

4.3.3. Parameter Estimation and Upper Limits

Following the frequentist prescription, the value of $\textsl {S}$ obtained from the search is compared with a threshold $\textsl {S}^*$ corresponding to a given false alarm probability (1% in this analysis). If $\textsl {S} > \textsl {S}^*$, then one has a potential signal detection deserving deeper study. In the case of signal detection, the signal parameters can be estimated from $\hat{H}_+,\, \hat{H}_{\times }$, using the relations shown in Appendix B. If the measured $\textsl {S}$ value lies below the threshold, we can set an upper limit on the amplitude of a possible signal present.

The determination of upper limits is carried out via Monte Carlo simulations similar to the limit determination described in Section 4.2.3. In the case of 4 dof, the unknown parameters, ψ and cos ι were taken to be uniformly distributed. The analysis method allows us to establish an upper limit for the wave amplitude H₀ defined in Appendix A. This was translated into an upper limit on h₀, under the assumption that the source is a triaxial neutron star, using Equation (A5) after maximizing the factor under the square root with respect to the inclination angle. In this way the upper limit we obtain is conservative. In the 2 dof case, we compute the upper limit by using for ψ and ι the values given in Equation (1).

The statistical error associated with the Monte Carlo simulations is estimated as half of the difference between the two signal amplitudes that bound the 95% confidence level. The grid in the amplitude of the injected signals has been chosen fine enough that the resulting statistical error is about one order of magnitude smaller than the systematic error coming from calibration and actuation uncertainty.

In the analyses, all available science mode data recorded by Virgo were used. No evidence for a CW signal was seen using any of the three analysis methods described in Section 4. We have therefore used the data to set upper limits on the GW amplitude.

For the complex heterodyne method (Section 4.1) the marginalized posteriors for the four parameters, using the two different priors, are shown in Figures 4 and 5. The presence of a detectable signal would show up as a posterior distribution in h₀ that is peaked away from h₀ = 0. The observed distributions are consistent with no signal being present. The 95% credible limits on h₀ are shown and have values 2.4 × 10⁻²⁴ and 2.1 × 10⁻²⁴, respectively (note that the strongly peaked distributions for cos ι and ψ in Figure 4 are simply the restricted priors placed on those parameters).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the ${\mathcal F}$ and ${\mathcal G}$ statistics (Section 4.2), the values obtained were consistent with false alarm probabilities of 22% and 35%, respectively. Since these probabilities are far above our 1% false alarm threshold, we conclude that the data are consistent with the absence of a signal. Using the Monte Carlo method described in Section 4.2.3, we set 95% confidence upper limits on h₀ of 2.4 × 10⁻²⁴ and 2.2 × 10⁻²⁴, respectively.

For the matched filter (MF) on Fourier components (Section 4.3), the values computed for the 4 dof and 2 dof statistics were consistent with false alarm probabilities of 46% and 40%, respectively. Again we conclude that the data are consistent with the absence of a signal. We obtain 95% confidence upper limits on h₀ of 2.2 × 10⁻²⁴ and 1.9 × 10⁻²⁴, respectively.

The results for all three analyses are summarized in Table 3, which also includes the systematic uncertainty in the upper limit from calibration and actuation uncertainties. For each analysis, results are given both for the case in which ψ and cos ι are assumed to be known (i.e., with restricted priors) and unknown (i.e., with unrestricted priors).

Table 3. Estimated 95% Upper Limit on h₀ for PSR J0835−4510 from the Three Different Analysis Methods (the Horizontal Line Separates Bayesian from Frequentist Results)

Analysis Method	95% Upper Limit for h0
Heterodyne, restricted priors	(2.1 ± 0.1) × 10−24
Heterodyne, unrestricted priors	(2.4 ± 0.1) × 10−24
${\mathcal G}$ statistic	(2.2 ± 0.1) × 10−24
${\mathcal F}$ statistic	(2.4 ± 0.1) × 10−24
MF on signal Fourier components, 2 dof	(1.9 ± 0.1) × 10−24
MF on signal Fourier components, 4 dof	(2.2 ± 0.1) × 10−24

Notes. The systematic error on amplitude from calibration and actuation amounts to ∼5.5%, as discussed in Section 3. This corresponds to an uncertainty on the upper limits of about ±0.1 × 10⁻²⁴. For all upper limits the statistical error, associated with the Monte Carlo simulations used to establish the limit itself, is about one order of magnitude smaller.

Download table as:  ASCIITypeset image

We emphasize once again that the two results for the complex heterodyne method are Bayesian 95% credible limits on h₀, while the ${\mathcal G}$, ${\mathcal F}$, 2 dof, and 4 dof results are frequentist 95% confidence upper limits. While we would expect the two types of upper limit to be similar in value, they are not directly comparable, because they address different questions. The Bayesian question asks: "Given our priors and our data, for what value of h₀ are we 95% certain that any true signal lies below that value?" The frequentist question asks: "Above what value of h₀ would a signal produce a larger value of our statistic 95% of the time?" The subtle difference between these questions means that they may give different answers for the same data, and we should not read too much into the fact that in this search the two approaches gave very similar numbers.

5.1. Validation with Hardware Injections

All three pipelines used in the analysis have been tested with both software and hardware injections of CW signals in the VSR2 data. In particular, we discuss here hardware injections. For the entire duration of the run, 13 CW signals (named Pulsar0–12) have been injected in the Virgo detector by sending the appropriate excitations to the coils used to control one mirror's position. These signals were characterized by various amplitudes, spanned a frequency range from ∼20 Hz to ∼1400 Hz, and covered a range of values for the spin-down $\dot{f}$ from ∼ − 4 × 10⁻¹⁸ Hz s⁻¹ to ∼ − 2.5 × 10⁻⁸ Hz s⁻¹. The corresponding source position (α, δ), inclination ι of the source spin axis, and polarization angle ψ were chosen randomly. All the injected signals have been generated using the same software as the signals injected in LIGO S5 and previous runs. Injected signals, Pulsar0–9 have also the same parameters as the LIGO injections, while Pulsar10–12 have very low frequency and have been injected in Virgo only. The three pipelines were exercised on several of these simulated signals. The pipelines have been able to detect the signals and to estimate their parameters with good accuracy when the S/N is sufficient. In particular, in Tables 5–7 we report the results obtained for Pulsar3, characterized by a very small spin-down and high S/N, Pulsar5 with low frequency, very small spin-down, and relatively low S/N, and Pulsar8 with high spin-down and S/N. The frequency parameters for these three injections are given in Table 4. There is good agreement between the true and recovered signal parameters. With the method based on matched filtering on the signal Fourier components, the estimation of the signal absolute phase is not straightforward.

Table 4. Frequency and Positional Parameters for the Hardware Injections ($\ddot{f}=0$ for all the Injections)

Name	f (Hz)	$\dot{f}$ (Hz s−1)	α (deg)	δ (deg)	S/N
Pulsar3	108.8571594	−1.46 × 10−17	178.372574	−33.436602	192
Pulsar5	52.80832436	−4.03 × 10−18	302.626641	−83.8391399	40
Pulsar8	194.3083185	−8.65 × 10−9	351.389582	−33.4185168	197

Notes. The reference time epoch for the source frequency is MJD = 52944 for all the injections. The optimal signal-to-noise ratio (S/N) is also given.

Download table as:  ASCIITypeset image

Table 5. Estimated Parameters for Hardware Injection Pulsar3 from the Three Different Analysis Methods

Method	$\frac{h_{0,{\rm found}}}{h_{0,{\rm inj}}}$	ι [ιinj = 1.651]	ψ [ψinj = 0.444]	Φ0 [Φ0, inj = 5.53]
Heterodyne	0.97	1.67	0.43	5.55
${\mathcal F}$ statistic	0.96	1.65	0.44	5.54
MF on signal Fourier comp., 4 dof	0.96	1.66	0.44	*

Download table as:  ASCIITypeset image

Table 6. Estimated Parameters for Hardware Injection Pulsar5 from the Three Different Analysis Methods

Method	$\frac{h_{0,{\rm found}}}{h_{0,{\rm inj}}}$	ι(η) [ιinj = 1.089]	ψ [ψinj = −0.364]	Φ0 [Φ0, inj = 2.23]
Heterodyne	0.90	0.99	−0.27	2.05
${\mathcal F}$ statistic	0.89	0.98	−0.27	2.10
MF on signal Fourier comp., 4 dof	0.97	0.96	−0.26	*

Download table as:  ASCIITypeset image

Table 7. Estimated Parameters for Hardware Injection Pulsar8 from the Three Different Analysis Methods

Method	$\frac{h_{0,{\rm found}}}{h_{0,{\rm inj}}}$	ι(η) [ιinj = 1.497]	ψ [ψinj = 0.170]	Φ0 [Φ0, inj = 5.89]
Heterodyne	0.97	1.49	0.18	5.90
${\mathcal F}$ statistic	0.95	1.49	0.17	6.07
MF on signal Fourier comp., 4 dof	0.98	1.50	0.17	*

Download table as:  ASCIITypeset image

In this paper we present the results of the analysis of Virgo VSR2 run data for the search of continuous GW signals from the Vela pulsar. The data have been analyzed using three largely independent methods and assuming that the GW emission follows the radio timing. For an assumed known orientation of the star's spin axis and value of the polarization angle, two methods have determined frequentist upper limits at 95% confidence level of, respectively, 1.9 × 10⁻²⁴ and 2.2 × 10⁻²⁴. The third method has determined a Bayesian 95% degree-of-belief upper limit of 2.1 × 10⁻²⁴. The lowest of these is about 41% below the indirect spin-down limit. It corresponds to a limit on the star ellipticity of 1.1 × 10⁻³, which is well above the maximum equatorial ellipticity that a neutron star with a "standard" equation of state can sustain, but comparable to the maximum value permitted by some exotic equations of state (Owen 2005; Lin 2007; Haskell et al. 2007). Given that the power emitted in GW is $\dot{E}_{{\rm GW}}=-\frac{32\pi ^6G}{5c^5}I_{zz}^2\epsilon ^2 f^6$, our results constrain the fraction of spin-down energy due to the emission of GW to be below 35%. For an unknown orientation of the star's spin axis and polarization angle, the two frequentist upper limits are, respectively, 2.2 × 10⁻²⁴ and 2.4 × 10⁻²⁴ while the Bayesian upper limit is 2.4 × 10⁻²⁴. The lowest of these is about 33% below the spin-down limit. In this case the limit on the star ellipticity is 1.2 × 10⁻³, while the corresponding limit on the fraction of spin-down energy emitted through GW is 45%. These numbers assume the canonical value for the star moment of inertia, I = 10³⁸ kg m². However, the theoretically predicted values of I vary in the range ∼1–3 × 10³⁸ kg m² (Abbott et al. 2010), so our upper limit on the ellipticity can be considered as conservative. Such ellipticities could also be sustained by internal toroidal magnetic fields of order 10¹⁶ G, depending on the field configuration, equation of state, and superconductivity of the star (Akgun & Wassermann 2007; Haskell et al. 2008; Colaiuda et al. 2008; Ciolfi et al. 2010). Then, our results have constrained the internal toroidal magnetic field of the Vela to be less than of the order of that value (it must be stressed, however, that the stability of a star with an internal field much larger than the external one is still an open issue). Vela is the second young pulsar for which the spin-down limit has now been beaten.

A more stringent constraint on the emission of GW from the Vela pulsar may be established by analyzing data of the next Virgo+ run (VSR4) which is tentatively scheduled for summer 2011 and should last a few months. This run, assuming the planned sensitivity is reached, could be able to probe values of the Vela pulsar ellipticity below a few units in 10⁻⁴, corresponding to a fraction of spin-down energy emitted through the emission of GW below a few percent. We note that this run will also provide interesting results for several other low-frequency pulsars. In particular, it could allow detection of GW from the Crab pulsar and J1952+3252 if their ellipticities are larger than ∼10⁻⁵, a value nearly compatible with the maximum deformation allowed by standard neutron star equations of state.

Second-generation detectors are expected to have a still better sensitivity at low frequency. Advanced Virgo (Acernese et al. 2009) and Advanced LIGO (Harry & the LIGO Scientific Collaboration 2010), which should come into operation around 2014–2015, in one year could detect a GW signal from the Vela pulsar if its ellipticity is larger than a few times 10⁻⁵, the corresponding fraction of spin-down energy emitted through GW being below a few times 10⁻⁴ in this case.

The possibility of building a third-generation GW detector, with a sensitivity a factor of 10 or more better than advanced detectors in a wide frequency range, is also being studied. The Einstein Telescope (Punturo et al. 2010), which is currently at the stage of design study, is expected to release its first science data around 2025–2027. It should be able to detect GWs from the Vela pulsar, using one year of data, for ellipticity larger than 4 × 10⁻⁷ to 10⁻⁶, depending on the detector configuration that will be chosen.

We dedicate this paper to the memory of our friend and colleague Stefano Braccini, who made very important contributions to the development of the Virgo detector and, more recently, contributed to the search effort for CW signals with his usual enthusiasm and skillfulness.

The authors gratefully acknowledge the support of the United States National Science Foundation for the construction and operation of the LIGO Laboratory, the Science and Technology Facilities Council of the United Kingdom, the Max-Planck-Society, the State of Niedersachsen/Germany for support of the construction and operation of the GEO600 detector, and the Italian Istituto Nazionale di Fisica Nucleare and the French Centre National de la Recherche Scientifique for the construction and operation of the Virgo detector. The authors also gratefully acknowledge the support of the research by these agencies and by the Australian Research Council, the Council of Scientific and Industrial Research of India, the Istituto Nazionale di Fisica Nucleare of Italy, the Spanish Ministerio de Educación y Ciencia, the Conselleria d'Economia Hisenda i Innovació of the Govern de les Illes Balears, the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, the Polish Ministry of Science and Higher Education, the FOCUS Programme of Foundation for Polish Science, the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, The National Aeronautics and Space Administration, the Carnegie Trust, the Leverhulme Trust, the David and Lucile Packard Foundation, the Research Corporation, and the Alfred P. Sloan Foundation.

The continuous GW signal emitted by a generic rotating rigid star can be described by a polarization ellipse. The polarization ellipse is characterized by the ratio $\eta =\frac{a}{b}$ of its semi-minor to its semi-major axis and by the angle ψ defining the direction of the major axis. The angle ψ is the same introduced in Section 2. The ratio η varies in the range [ − 1, 1], where η = 0 for a linearly polarized wave and η = ±1 for a circularly polarized wave. The (complex) signal can be expressed as

$$\begin{equation} {\bf h}(t)=H_0\left(H_+ {\bf e_+} + H_{\times } {\bf e_{\times }}\right) e^{i\Phi (t)}, \end{equation} \tag{ A1 }$$

where e₊ and e_× are the two polarization tensors and the plus and cross amplitudes are given by

$$\begin{eqnarray} &&H_+=\frac{\cos {2\psi }-i \eta \sin {2\psi }}{\sqrt{1+\eta ^2}} \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} &&H_{\times }=\frac{\sin {2\psi }+i \eta \cos {2\psi }}{\sqrt{1+\eta ^2}}. \end{eqnarray} \tag{ A3 }$$

If we consider, as in Section 2, a triaxial neutron star rotating around a principal axis of inertia, the following relations among H₀, η and h₀, ι hold:

$$\begin{eqnarray} \eta & = & {-}\frac{2\cos {\iota }}{1+\cos {}^2{\iota }} \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} H_0 & = & \frac{h_0}{2}\sqrt{1+6\cos {}^2{\iota }+\cos {}^4{\iota }}. \end{eqnarray} \tag{ A5 }$$

In terms of + and × components we have

$$\begin{eqnarray} H_{+,\psi =0} & = & \frac{h_+}{H_0}=h_0 \frac{(1+\cos {}^2{\iota })}{2H_0} \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} \Im (H_{\times,\psi =0}) & = & {-}\frac{h_{\times }}{H_0}=-\frac{h_0 \cos {\iota }}{H_0}. \end{eqnarray} \tag{ A7 }$$

In this formalism the complex gravitational strain at the detector is given by

$$\begin{equation} h(t)=H_0 \left(A_+(t)H_++A_{\times }(t)H_{\times }\right) e^{i\Phi (t)}, \end{equation} \tag{ A8 }$$

where

$$\begin{eqnarray} &&A_+=F_+(\psi =0) \end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray} &&A_{\times }=F_{\times }(\psi =0). \end{eqnarray} \tag{ A10 }$$

After Doppler and spin-down corrections, as described in Section 4.3, we have

$$\begin{equation} h(t)=H_0\left(A_+(t)H_++A_{\times }(t)H_{\times }\right) e^{i \left(\omega _0t+\Phi _0\right)}. \end{equation} \tag{ A11 }$$

We now introduce the signal 5-vectors for the + and × components, $\boldsymbol{A^+}, \boldsymbol{A^{\times }}$, given by the Fourier components, at the five frequencies produced by the amplitude and phase modulation, of the detector response functions A₊, A_×. It is straightforward to see that the signal in the antenna is completely defined by the 5-component complex vector

$$\begin{equation} \boldsymbol{A}=H_0e^{i\Phi _0}(H_+ \boldsymbol{A^+}+H_{\times } \boldsymbol{A^{\times }}). \end{equation} \tag{ A12 }$$

More details can be found in Astone et al. (2010).

Once the two estimators $\hat{H}_+, \hat{H}_{\times }$ have been computed from the data, if a detection is claimed, the signal parameters H₀, η, ψ can be estimated using the following relations. The estimator of the signal amplitude is given by

$$\begin{equation} \hat{H_0}=\sqrt{|\hat{H}_+|^2+|\hat{H}_{\times }|^2}. \end{equation} \tag{ B1 }$$

Introducing the quantities

$$\begin{eqnarray} &&\hat{H}_+ \cdot \hat{H}_{\times }=A+i B, \end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray} &&|\hat{H}_+|^2-|\hat{H}_{\times }|^2=C, \end{eqnarray} \tag{ B3 }$$

where the scalar product is between two complex numbers and includes a complex conjugation of one, the estimation of the ratio between the axes of the polarization ellipse is

$$\begin{equation} \hat{\eta }=\frac{-1+\sqrt{1-4B^2}}{2B}, \end{equation} \tag{ B4 }$$

while the estimation of the polarization angle can be obtained from

$$\begin{eqnarray} &&\cos {(4\hat{\psi })}=\frac{C}{\sqrt{4A^2+B^2}} \end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray} &&\sin {(4\hat{\psi })}=\frac{2A}{\sqrt{4A^2+B^2}}. \end{eqnarray} \tag{ B6 }$$

We assume that the noise in the data is Gaussian and uncorrelated. In order to take into account non-stationarity of the data, we assume that each noise sample n(l) in the data time series is drawn from a Gaussian distribution with a variance σ²(l). We assume that the Gaussian distributions in question have zero means. Thus, the autocorrelation function K(l, l') for the noise is given by

$$\begin{equation} K(l,l^{\prime }) = \sigma ^2(l)\,\delta _{ll^{\prime }}, \end{equation} \tag{ C1 }$$

where l, l' are integers and where $\delta _{ll^{\prime }}$ is Kronecker's delta function. Let us first assume that the signal h(l) is completely known and that the noise is additive. Thus, when the signal is present the data take the following form:

$$\begin{equation} x(l) = n(l) + h(l). \end{equation} \tag{ C2 }$$

For Gaussian noise the optimal filter q(l) is the solution of the following (integral) equation (see Jaranowski & Królak 2009, p. 72),

$$\begin{equation} h(l) = \sum ^N_{l^{\prime }=1} K(l,l^{\prime }) q(l^{\prime }), \end{equation} \tag{ C3 }$$

where N is the number of data points. Consequently, we have the following equation for the filter q(l):

$$\begin{equation} q(l) = \frac{h(l)}{\sigma ^2(l)}, \end{equation} \tag{ C4 }$$

and the following expression for the log likelihood ratio ln Λ:

$$\begin{equation} \ln \Lambda [x] = {\left\langle h\, x\right\rangle } - \frac{1}{2} {\langle h^2\rangle }, \end{equation} \tag{ C5 }$$

where the operator 〈 · 〉 is defined as

$$\begin{equation} {\left\langle g\,f\right\rangle } = \sum ^N_{l=1}\frac{g(l) f(l)}{\sigma ^2(l)}. \end{equation} \tag{ C6 }$$

Thus, we see that for non-stationary Gaussian noise with the autocorrelation function (C1) the optimal processing is identical to matched filtering for a known signal in stationary Gaussian noise, except that we divide both the data and the filter by time-varying standard deviation of the noise. This may be thought as a special case of whitening the data and then correlating it using a whitened filter. The method is essentially the same as the Wiener filter introduced in Section 4.3. The generalization to the case of signal with unknown parameters is immediate.

In the analysis we use complex heterodyne data x_het,

$$\begin{equation} x_{{\rm het}}(l) = x(l)e^{-i\Phi _{{\rm het}}(l)}, \end{equation} \tag{ C7 }$$

where Φ_het is the heterodyne phase (Φ_het can be an arbitrary real function). Thus, we rewrite the ${\mathcal F}$ and ${\mathcal G}$ statistics and amplitude parameter estimators using complex quantities. We introduce complex amplitudes A_a and A_b,

$$\begin{eqnarray} A_a &=& A_1 + i A_3, \end{eqnarray} \tag{ C8 }$$

$$\begin{eqnarray} A_b &=& A_2 + i A_4, \end{eqnarray} \tag{ C9 }$$

where the amplitudes A_k, k = 1, 2, 3, 4 are defined by Equation (23) of Jaranowski & Królak (2010) and we also introduce the complex filters:

$$\begin{eqnarray} h_a(l) &=& a(l) e^{-i [\Phi (l) - \Phi _{{\rm het}}(l)]}, \nonumber\\ h_b(l) &=& b(l) e^{-i [\Phi (l) - \Phi _{{\rm het}}(l)]}, \end{eqnarray} \tag{ C10 }$$

where a and b are amplitude modulation functions (see Equations (7) and (8)) defined by Equations (12) and (13) in Jaranowski et al. (1998), and Φ(l) is the phase defined by Equation (9).

The ${\mathcal F}$ statistic takes the following form:

$$\begin{eqnarray} &&\fontsize{13}{15}\selectfont\hbox{${\mathcal F}=\!\frac{{\langle b^2\rangle }|{\langle x_{{\rm het}}\, h_a\rangle }|^2 + {\langle a^2\rangle } |{\langle x_{{\rm het}}\, h_b\rangle }|^2 - 2 {\langle a\, b\rangle } \Re ({\langle x_{{\rm het}}\, h_a\rangle } {\langle x_{{\rm het}}\, h_b\rangle }^*)}{{\langle a^2\rangle }{\langle b^2\rangle } - {\langle a\, b\rangle }^2},$}\nonumber\\ \end{eqnarray} \tag{ C11 }$$

and the complex amplitude parameter estimators are given by

$$\begin{eqnarray} \widehat{A}_a &=& {\displaystyle 2\frac{{\langle b^2\rangle } {\langle x_{{\rm het}}\, h_a\rangle }^* - {\langle a\, b\rangle }{\langle x_{{\rm het}}\, h_b\rangle }^*}{{\langle a^2\rangle }{\langle b^2\rangle } - {\langle a\, b\rangle }^2}},\nonumber\\ \widehat{A}_b &=& {\displaystyle 2\frac{{\langle a^2\rangle } {\langle x_{{\rm het}}\, h_b\rangle }^* - {\langle a\, b\rangle }{\langle x_{{\rm het}}\, h_a\rangle }^*}{{\langle a^2\rangle }{\langle b^2\rangle } - {\langle a\, b\rangle }^2}}. \end{eqnarray} \tag{ C12 }$$

In the case of the ${\mathcal G}$ statistic, it is useful to introduce a complex amplitude A,

$$\begin{equation} A = A_c + i A_s, \end{equation} \tag{ C13 }$$

where real amplitudes A_c and A_s are defined by Equation (20) of Jaranowski & Królak (2010) and a complex filter h_g,

$$\begin{equation} h_g = (h_c + i h_s)e^{i\Phi _{{\rm het}}}, \end{equation} \tag{ C14 }$$

where real filters h_c and h_s are defined by Equation (7) of Jaranowski & Królak (2010). In complex notation, the ${\mathcal G}$ statistic assumes the following simple form (cf Equation (18) of Jaranowski & Królak 2010):

$$\begin{equation} {\mathcal G}= \frac{|{\langle x_{{\rm het}}\, h_g\rangle }|^2}{2 D}, \end{equation} \tag{ C15 }$$

where

$$\begin{equation} D = {\langle |h_g|^2\rangle }. \end{equation} \tag{ C16 }$$

The estimator of the complex amplitude A is given by

$$\begin{equation} \widehat{A} = \frac{{\langle x_{{\rm het}}\, h_g\rangle }}{D}. \end{equation} \tag{ C17 }$$

The Grubbs test (Grubbs 1969) is used to detect outliers in a univariate data set. Grubbs' test detects one outlier at a time. This outlier is removed from the data set and the test is iterated until no outliers are detected.

Grubbs' test is a test of the null hypothesis:

H0: There are no outliers in the data set x_i.

against the alternate hypotheses:

H1: There is at least one outlier in the data set x_i.

Grubbs' test assumes that the data can be reasonably approximated by a normal distribution.

The Grubbs test statistic is the largest absolute deviation from the sample mean in units of the sample standard deviation and it is defined as

$$\begin{equation} G = \frac{\mbox{max}|x_i-\mu |}{\sigma }, \end{equation} \tag{ D1 }$$

where μ and σ denote the sample mean and standard deviation, respectively.

The hypothesis of no outliers is rejected if

$$\begin{equation} G > \frac{n - 1}{\sqrt{n}}\sqrt{\frac{t^2_{\alpha /(2 n), n - 2}}{n - 2 + t^2_{\alpha /(2 n), n - 2}}}, \end{equation} \tag{ D2 }$$

with t_{α/(2n), n − 2} denoting the critical value of the t-distribution with n − 2 dof and a significance level of α/(2n).

We have applied the Grubbs test to the coarse heterodyne data before analyzing them with ${\mathcal F}$ and ${\mathcal G}$ statistics. We have applied the test to segments of 2¹⁶ data points and we have assumed false alarm probability of 0.1%. This resulted in identification of 13,844 outliers from the original data set containing 12,403,138 points. We replaced these outliers with zeros. The time series before and after the removal of the outliers are presented in Figure 6. The number of outliers constitutes 0.1% of the total data points in input data resulting in a negligible loss of S/N of any continuous signal present in the data. With different methods to identify the outliers used by other searches the number of outliers was similar.

Zoom In Zoom Out Reset image size