A search for periodic gravitational waves from three recycled pulsars using the AURIGA detector

A. Mion; M. J. Benacquista; M. Kramer; P. P. C. Freire; A. Possenti

doi:10.1051/0004-6361/200911946

Home

All issues

Volume 504 / No 2 (September III 2009)

A&A, 504 2 (2009) 673-679

Full HTML

Free Access

Issue		A&A Volume 504, Number 2, September III 2009


Page(s)		673 - 679
Section		Astronomical instrumentation
DOI		https://doi.org/10.1051/0004-6361/200911946
Published online		15 July 2009

An implementation of a modified version of the unified approach method

A. Mion^1, - M. J. Benacquista² - M. Kramer³ - P. P. C. Freire⁴ - A. Possenti⁵

1 - Dipartimento di Fisica, Università di Trento and INFN, Gruppo Collegato di Trento, Sezione di Padova, 38050 Povo, Trento, Italy
2 - Department of Physics and Astronomy, University of Texas at Brownsville, Brownsville, TX, USA
3 - University of Manchester, Jodrell Bank Observatory, Macclesfield, Cheshire, SK11 9DL, UK
4 - West Virginia University, Department of Physics, PO Box 6315, Morgantown, WV, 26506, USA
5 - INAF-Osseratorio di Cagliari, loc. Poggio dei Pini, strada 54, 09012 Capoterra, Italy

Received 25 February 2009 / Accepted 4 July 2009

Abstract
Aims. We report on a search for continuous gravitational wave emission from three recycled radio pulsars, performed by using the data of the resonant detector AURIGA. Given the spin rate of the selected targets - the isolated pulsar PSR J1939+2134 and the binary pulsars PSR J0024-7204J and PSR J0218+4232 - the expected frequency of the emitted gravitational waves falls in the high sensitivity range of the detector.
Methods. The main topic is the method, meaning that the statistical analysis is performed by implementing a slightly modified version of the Feldman and Cousins Unified Approach.
Results. By using ephemerides provided by suitable radio observations of the targets, we were able to demodulate the Doppler shifts within a coherence time of 1 day, and then incoherently sum 10 daily spectra collected from December 8th to December 17th, 2006. We have found upper limits on the gravitational wave amplitudes in the order of a few units of 10^-23 at 90% Confidence Level (C.L.), which translate to limits in the ellipticity of the targeted pulsars of $\varepsilon < 10^{-4}$ at 90% C.L.
Conclusions. The same framework can then be applied to data coming from most sensitive experiments as VIRGO or LIGO; moreover, an application to recently discovered transients in X-ray pulsars is discussed.

Key words: gravitational waves - pulsars: general - methods: data analysis - methods: statistical - techniques: radial velocities

1 Introduction

Continuous, high-frequency gravitational waves (GWs) sources are good candidates for detection with resonant mass detectors. Even though their signals are expected to have smaller amplitude with respect to other GW emission phenomena in the Universe, such as supernovae or gamma-ray bursts (GRBs) (Cutler & Thorne 2001), the data can be folded over long time intervals in order to increase the signal to noise ratio.

In particular, a promising class of sources for the resonant mass detectors are the so called millisecond radio pulsars (MSPs). According to common wisdom, they are relatively old (age typically longer than at least 10⁸ yr) neutron stars (NSs) spun up to high rotational rate (typically larger than 100 Hz) during a stage of mass and angular momentum transfer from a companion star in a binary system (Alpar et al. 1982). During this process - often referred to as ``recycling'' - the NS is believed to appear as an accreting X-ray source belonging to a low mass X-ray binary system (LMXB) (Battacharya & van den Heuvel 1991).

The occurrence of GW emission from MSPs and from NSs in LMXBs has both theoretical and observational support. As to the theory, there are several mechanisms for a rapidly rotating NS to emit GWs. The two most relevant (Jones 2002; Owen 1999) are (i) the presence of a time-varying NS quadrupole moment (in a description considering the NS as a rotating rigid body); and (ii) the occurrence of the so-called r-modes. In the first model, the NS is approximated as a non-spherical rigid object, rotating at an angular velocity $\nu_0\approx 1~{\rm KHz}$ ; the ellipticity parameter $\varepsilon$ is defined as $\varepsilon = \left(I_1-I_2\right)/I_3$ where the I_k are the principal inertia moments of the star, and presumably ranges from 10^-7to 10^-9 (Jones 2002; Ruderman 2006). In this case, one of the possible GW emission channels is characterized by the angular frequency $\omega_{\rm gw}=2\omega_0$ , where $\omega_0=2\pi\nu_0.$ The relation between the gravitational strain amplitude h₀ and the pulsar's parameters (angular frequency $\omega$ , distance d, moment of inertia I and ellipticity $\varepsilon$ ) is

$\begin{displaymath}h_0 = \frac{4G}{c^4 d} \varepsilon I \omega^2. \end{displaymath}$

(1)

The second emission mechanism involves the so-called r-modes, which have been widely studied in the literature (Anderson et al. 2001). Assuming that the NS is spherically symmetric and that second order terms can be neglected in the Poisson equation describing the stellar fluid (Anderson & Kokkotas 2001), r-modes exhibit linear relations between the rotational and gravitational frequencies, i.e. $\omega_{\rm gw}=(4/3)\omega_0$ , $(2/3)\omega_0$ . These GW emission frequencies are exact only if the linearized equations of fluid dynamics hold, and this can in principle be the case, for binary systems where accretion occurs at a stable rate.

On the side of the observations, it has been noticed that the spins of both the known accretion powered NSs included in LMXBs (Bildsten 1998) and the spins of the known rotation powered MSPs (the most rapidly spinning of which is rotating at 716 Hz, Hessels et al. 2006) tend to be clustered at frequencies significantly below the break up frequency imposed by most of the proposed Equations of State for the nuclear matter. Binary evolution effects may play a role for that (Possenti et al. 2000), but another viable possibility is that GW emission sets in at high spin rates, preventing the achievement of a rotational period close to the mass shedding spin limit (Bildsten 1998; Levin & Ushomirsky 2001).

In this paper we report on a search for continuous GW emission from three recycled MSPs, performed by using the data of the resonant detector AURIGA. The selection of the sample of targeted sources accounts for the bandwidth of the detector and is described in Sect. 2.

As anticipated above, a data folding procedure is applied in order to extract a GW signal from the detector noise. In Sect. 3 it is first reported on the algorithm used to take into account the AURIGA motion in the Solar System Barycenter frame (SSB), and the relative motion of the neutron star with respect to the barycenter of the binary system. In the same section the implemented method for the data analysis is also presented.

Section 4 contains the discussion of the statistical frame (the unified approach) used to interpret the results. We show here, as a new feature, as a slightly modified version of the Feldman and Cousins method can be used to search for continuous signals. As we were not able to reject the null hypothesis, we can only set upper limits on the associated GW amplitudes. The statistical methods used to determine the upper limits involve the construction of a modified version of the Feldman and Cousins Confidence Belt, and this is also discussed in some detail in Sect. 4.

In Sect. 5 we will finally draw our conclusions and discuss future prospects in this field opened by the upcoming Advanced LIGO/VIRGO detectors. In particular, we suggest the application of the algorithms presented in this paper to the sample of the accreting X-ray millisecond pulsars (AMXPs), belonging to transient LMXBs.

2 The selection of the sample

AURIGA is a resonant mass gravitational wave detector with an usable bandwidth going from $850~{\rm Hz}$ to $960~{\rm Hz}$ (Baggio et al. 2005). We have inspected the ATNF pulsar catalogue (Manchester et al. 2005) searching for sources whose gravitational radiation emission could occur at a frequency $\omega_{\rm gw}$ in the range above. For the GW emission mechanisms described in Sect. 1, we found five MSPs with a suitable spin frequency.

The folding procedure described below requires the knowledge of the positional, rotational and, when appropriate, binary parameters of the targeted NSs, valid over the time interval of collection of the AURIGA data (i.e. from December 8th, 2006 until December 17th, 2006). These parameters were retrieved from timing observations in the radio band performed at the 64-m Parkes Telescope (Australia) and at the 76-m Lovell Telescope at Jodrell Bank (UK), over a $\sim$ 1 yr dataspan including the 10 days of the AURIGA runs.

We note that the simple extrapolation of the parameters from the ephemeris reported in the aforementioned pulsar catalog, which is built on the basis of published observations usually taken many years before the AURIGA runs, may not always be a viable choice. This is particularly true for sources orbiting in a tight orbit with a non degenerate companion. In this case, matter irregularly released from the companion and/or tidal effects can significantly affect the orbital parameters of the binary, hampering the extrapolation of a reliable timing solution over long time intervals. Despite MSPs usually being much less affected by intrinsic timing noise than other kinds of pulsars, some of them also displayed rotational irregularities, like small glitches (Cognard & Backer 2004). In summary, monitoring of the sources (both the binary and the isolated MSPs) in a time interval bracketing the runs of the GW detector is the safest procedure for optimizing the capability of the GW search analysis.

In view of that, we had to exclude two of the five objects in the original list: PSR J1701-3006F and PSR J0024-7204W, for which no viable timing solution for the dates including the AURIGA runs was available. In the end we have been left with PSRs J0024-7204J, J0218+4232 and J1939+2134. We checked that their spin frequencies fall in frequency sub-bands where the AURIGA noise is well behaved (Gaussian and stationary) (Vinante 2006), which is a requirement for the following statistical analysis of the results. The spin frequency $\nu _0$ of the 3 selected targets and the frequency $\nu_{\rm gw}$ of the gravitational waves that we searched for are reported in Table 1.

Table 1: $\nu _0$ is the spin frequency of the targeted MSPs at the time of the AURIGA run, whereas $\nu _{\rm gw} {\rm [Hz]}$ is the frequency of the GW signal we searched for.

In particular, PSR J0024-7204J is a binary pulsar discovered on 1991 (Manchester et al. 1991) in the globular cluster 47 Tucanae. It orbits a very low mass (minimum mass 0.021 ${M_\odot}$ ) companion in a ${\sim}0.12$ day almost circular (eccentricity ${\sim}4\times 10^{-5}$ ) orbit. The ephemerides for the epoch of interest have been retrieved from observations performed at the Parkes radio telescope at a center frequency of 1390 GHz. The total 256 MHz bandwidth has been split in 512 0.5-MHz wide channels per polarization, in order to minimize the effects of the interstellar dispersion. After having been summed in polarization and digitized every 80 $\mu$ s, the resulting 512 data streams (accurately tagged in time) have been de-dispersed and folded off-line in order to produce pulse profiles. Topocentric pulse time of arrivals (ToAs) were determined by convolving these profiles with a template pulse profile of high signal to noise ratio and then analyzed using TEMPO. It converts the topocentric ToAs to solar-system barycentric ToAs at infinite frequency (using the DE405 solar-system ephemeris ftp://ssd.jpl.nasa.gov/pub/eph/export/) and then determines the pulsar (positional, rotational and binary) parameters using a multi-parameter fit.

PSR J0218+4232 is a binary millisecond pulsar discovered on 1995 Navarro et al. (1995) in the Galactic field. Also in this case the orbit is nearly circular (eccentricity ${\sim}7\times 10^{-6}),$ but the orbital period is significantly longer ( ${\sim}2.03$ days) and the companion (likely a white dwarf) more massive (minimum mass of 0.167 ${M_\odot}$ ) than for the case of J0024-7204J. The ephemeris were obtained from regular timing observations carried on at a central frequency of 1400 MHz at the Jodrell Bank Observatory. A cryogenic receiver mounted on the 76-m Lovell Telescope provided 32 MHz bandwidth (split in 32 1-MHz wide channels) over two hands of circular polarization. Pulsar profiles were formed every 1-min sub-integrations and then were added in polarization pairs and combined to produce a single total-intensity profile. This was then convolved with a template derived from a single high signal-to-noise ratio profile at the same frequency to give a topocentric ToA. The pulsar parameters were finally determined from the ToAs with the same procedure described above.

Contrary to the other two selected targets, PSR J1939+2134 (also known as PSR B1937+21) is an isolated MSP, the first detected ever (Backer et al. 1982) and ranked as the most rapidly rotating pulsar until 2006 (Hessels et al. 2006). The positional and rotational parameters of the source at the time of the AURIGA runs resulted from timing observations conducted at the Jodrell Bank Observatory with the same instrumental set up as for PSR J0218+4232.

At the time when the paper was prepared, the ephemerides for these PSRs where not yet published in the ATNF catalogue, and this is why we don't show them here.

3 Data analysis

The instantaneous frequency and phase of a GW impinging upon the GW detector are of course affected by both the detector and source motions. A basic step in the data analysis is then to report frequency and phase to a suitable reference frame, which has been chosen to be the Solar System Barycenter (SSB).

Firstly, we computed the detector motion by means of a code which makes use of the freely available NOVAS routine package http://aa.usno.navy.mil/software/novas/novas_c/novasc_info.php and the ephemeris file DE405 by JPL. Then, the source motion along its orbit - when necessary - was calculated from the radio ephemeris, using procedures largely applied in the radio astronomy community (Manchester & Taylor 1974).

The problem is basically to exactly implement the matched filter where the prior knowledge of pulsar parameters is used to extract the signal power. The GW phase doesn't increase linearly with the time, but can be conveniently expressed as $\phi_0 + \omega_0 t + \hat{\phi}(t)$ , where $\phi_0$ is the unknown initial GW phase, $\omega_0$ is the intrinsic frequency at some initial time t₀, t is the time elapsed from time t₀, and $\hat{\phi}(t)$ represents the phase shift due to both intrinsic frequency derivatives from the ephemerides files and the Doppler effects. For each MSP, the GW signals to extract from AURIGA data read

$\begin{displaymath}s(t) = M(t) \cos{(\phi_0 + \phi(t))} , \end{displaymath}$

(2)

where $\phi(t)=\omega_0 t + \hat{\phi}(t)$ is the estimated phase evolution (i.e. the known part of the phase) and M(t) is an amplitude modulation due to the AURIGA antenna pattern. Moreover, in Eq. (2) there will generally be a suppressing factor due to the generally non optimal orientation of the pulsar rotation axis. This representation in which s(t) is factorized in two terms is interesting when the timescales of variation of the two terms are much different one from each other. M(t) is the slow-varying term and takes account of the changes in the orientation between the source and the detector.

The angle between this axis and the direction of sight is unknown and this uncertainty results in a scale factor to be applied to Eq. (2). The signal template in Eq. (2) reproduces the physical GW signals when quasi-stationary approximation holds, i.e. one has also to require slowness of the function M(t) and smoothness of the functions $\phi(t)$ and $\omega(t)=\omega_0+ {\rm d} \hat{\phi}(t)/{\rm d}t$ . For our targets it results $\dot{\phi}(t) - \omega_0< 2 \times 10^{-5}~\rm s^{-1}$ . Moreover, from the bar antenna pattern (Misner et al. 1973), we can write

$\begin{displaymath}M(t)=\sum_{k=0}^{4} M_k \cos(k \omega_\oplus t+c_k) \end{displaymath}$

(3)

where $\omega_\oplus$ is the Earth rotation frequency, whereas amplitudes M_k, and phase shifts c_k are constants which depend on the detector latitude and longitude and on the source right ascension and declination. Clearly, as the antenna pattern varies with the period of 1 sidereal day we have $\dot{M}(t)/M < \omega_\oplus = 1.16 \times 10^{-5} \ \rm s^{-1}$ . Thus, the amplitude modulation due to the antenna pattern factors can be taken into account by the calculation of averaged values of the antenna pattern during the observation time.

Table 2: Antenna Pattern factors.

The next step of the analysis is to represent s(t) as

$\begin{displaymath}s(t) = A(t) \cos \phi (t) - B(t) \sin \phi (t) \end{displaymath}$

(4)

where we have introduced the two time-dependent functions

$\begin{displaymath}A(t) = M(t) \cos \phi_0 \end{displaymath}$

(5)

and

$\begin{displaymath}B(t) = M(t) \sin \phi_0 \end{displaymath}$

(6)

where A(t) and B(t) are the two components of the signal in quadrature with respect to each other. The optimal filtering problem, i.e. the matching of the template s(t) to the detector data is then reduced to the application to the data of a locking filter at the variable frequency $\omega(t)$ . What we did is to consider M(t) as a slow modulation multiplying the frequency and amplitude modulation, and this results in a small signal to noise ratio (SNR) loss, as we now show. In fact, the maximum SNR would be obtaining applying the exact matched filter. In case of Gaussian and stationary noise it is

$\begin{displaymath}{\it SNR}^2 = \int \frac{\vert s(\nu)\vert^2}{S_h (\nu)} {\rm d}\nu \end{displaymath}$

(7)

where $\nu$ is the frequency and $S_h(\nu)$ is the noise power spectrum in equivalent gw strain. In our case the signal always stays in a very narrow band around a frequency $\nu _0$ where S_h is constant, so

$\begin{displaymath}{\it SNR}^2 = \frac{\int \vert s(\nu)\vert^2 {\rm d}\nu}{S_h(\nu_0)}\cdot \end{displaymath}$

(8)

The integrals over frequencies can be substituted with integrals over the time. Assuming stationary, Gaussian, and white noise within the signal bandwidth ( ${<}1~{\rm Hz}$ ) we calculated the SNR loss due to our approximation and result are reported in Table 2.

$\begin{figure} \par\includegraphics[width=14cm,clip]{aa11946-09f1.eps} \end{figure}$

Figure 1:

Complete data-analysis pipeline. The h-reconstructed signal is passed trough a band-pass digital filter before the Doppler correction. This allows to decimate the samples to deal with a reasonable amount of data. Then, after the Fourier transformation, all the signal, if present, is in the bin 0 of the spectrum. The other bins are used to infer the statistics of the noise, thus allowing to calculate the confidence intervals for the GW amplitude.

Open with DEXTER

It is worth noticing that we are not able to estimate the initial phase $\phi_0$ by means of radio timing observations and so both components in phase and in quadrature are present.

The Fourier transform of s(t) is the convolution of the Fourier transforms of A(t) and $\cos\phi(t)$ , minus the convolution of the Fourier transforms of B(t) and $\sin\phi (t)$ , and within the constant phase approximation it reads

$\displaystyle \widetilde{s}(\omega) \simeq \frac{1}{2} [\widetilde{A}(\omega - ... ...2\rm i} [\widetilde{B}(\omega - \omega_0) - \widetilde{B} (\omega + \omega_0)].$

(9)

In order to extract the signal from the detector's output, we provided the signal s(t) as input to a code that first creates a copy of the time series. Then it puts (see the pipeline in Fig. 1) the signal multiplied by $\cos\phi(t)$ in a first channel, whereas a second channels stores the signal multiplied by $\sin\phi (t)$ . The first channel s₁ will contain the time series

s₁(t)	=	$\displaystyle A (t)\cos^2\phi (t) -B(t)\sin\phi (t) \cos\phi (t)$
	=	$\displaystyle \frac{1}{2}{A(t)[1+\cos 2 \phi (t)] - \frac{1}{2} B(t)\sin 2\phi (t)}$	(10)

and its Fourier transform is

$\displaystyle s_1(\omega) = \frac{1}{2}A(\omega) + \frac{1}{4}A(\omega - 2 \ome... ...ga + 2 \omega_0) - \frac{1}{4}B(\omega - 2 \omega_0) - B (\omega + 2 \omega_0).$

(11)

Equation (11) shows that it is now possible to extract the signal: remembering that A and B are functions that slowly vary with time, the only important components of their Fourier transforms $A(\omega)$ and $B(\omega)$ are the one characterized by the condition $\omega \ll \omega_0$ . For these frequencies, $\vert\omega \pm 2 \omega_0\vert \gg \omega$ . So, if we apply a suitable numerical low-pass filter, with a time constant

$\begin{displaymath}\tau \gg \frac{1}{2 \omega_0}~ \end{displaymath}$

(12)

the product of Eq. (11) and the low pass filter, performed in the frequency domain, is

$\begin{displaymath}s_1^f (\omega) = \frac{A(\omega)}{2 (1 + i \omega \tau)}\cdot \end{displaymath}$

(13)

If the filter is chosen to have a rapid decay, i.e. the filter is a square box with a certain amplitude around the 0 frequency, leads to

$\begin{displaymath}s_1^f (\omega) \simeq \frac{A(0)}{2} \cdot \end{displaymath}$

(14)

Similar considerations would lead to

$\begin{displaymath}s_2^f (\omega) \simeq \frac{B(0)}{2} \end{displaymath}$

(15)

for the second channel s₂. It's clear that this approach completely solves the problem of Doppler modulations. The coorections of each sample for its proper Doppler factor is unnecessary. We compute the instantaneous frequency at the start time of each frame of data, and then, inside each frame, the phase is taken to evolve linearly with that frequency. This does not lead to errors, provided the duration of each frame is short enough (in our case, about 80.5 s). However, the duration of each coherent sub-search is 1 day; the frame duration is not the time of each coherent search but only the dimension of the buffers taken by the acquisition. For the following steps required for the analysis, we implemented the code using the MATLAB environment. Finally, the amplitude of the signal H was simply estimated as

$\begin{displaymath}H = 2 \sqrt{s_1^f(0)^2 + s_2^f(0)^2} . \end{displaymath}$

(16)

The outline of the method we have now discussed is represented in Fig. 1.

$\begin{figure} \par\includegraphics[width=11.5cm,]{11946fg2.eps} \end{figure}$	Figure 2: The AURIGA strain sensitivity (gray) compared with the expected (black curve) sensitivity at $4.5\rm~ K$ .
Open with DEXTER

4 Upper limits

We chose to run the analysis on 10 segments - 1 day long each - of contiguous data (MJD between 54077-54086), during the scientific run numbered 852. This period was chosen because the behavior of AURIGA was satisfactory, namely the noise was low and stable. This means that during the run of our choice the white noise level is better compared with the ones of other scientific runs. To extend the duration of this search to more segments would not improve so much the results, given the fact that in this case the result in terms of upper limit on GW emission would scale only with the fourth root of the number of segments. Also, we checked that in this period the rotational behavior of our sources was well reproduced by the radio timing data. A coherent search over the whole 10 days period was impossible, because the stability of the poles of the transfer function of AURIGA is not guaranteed for such a long time. So we opted for performing 10 coherent sub-searches (each lasting 1 day) and then incoherently summing (i.e. averaging) the results, thus ending with a unique spectrum for each target. The continuous component of the final spectrum (i.e. the bin 0 of a Fast Fourier Transform) holds the signal, if present. The noise spectral density is shown in Fig. 2. The plot represents the quantity

$\begin{displaymath}S^{1/2}_{hh}(\omega) = \sqrt{\frac{\vert V(\omega)\vert^2}{\vert T(\omega)\vert^2}} , \end{displaymath}$

(17)

where $V(\omega)$ is the Fourier transform of the output data v(t) and $T(\omega)$ is the transfer function.

All the bins in the final spectrum except the bin 0 (which is special because it could hold the signal) can be used in order to test the null hypothesis ${\cal H}_0$ that the signal is not present in the data. Because of the performed operation of spectrum averaging, the distribution of the frequency counts we find is a non-centered $\chi^2$ , with some non-centrality parameter $\theta$ holding the information about the signal energy. So, $\theta$ is the parameter we want to estimate. The first part of the hypothesis test consists of arbitrarily setting a false alarm probability, in our case, we choose a false alarm probability of 10^-2. When the signal is present our bin 0 is made by the sum of the energies of the noise and of the quantity h₀² T² where T is the coherence-time $T = 86~400~\rm s$ and h₀ is the average signal amplitude. Since the standard deviation $\sigma$ assumes a different value for each MSP (because each MSP belongs to a different band in the spectrum which has its own specific variance) we end up with a different value of the threshold for each MSP. Let x be the generic energy in some bin. Let's call $x_{\rm FA}$ this calculated threshold on x. Finally, we have to define a procedure in order to set a confidence interval, either upper-limit or two-sided, on the measured GW amplitude h₀. The way we choose, to measure h₀ and to know the statistical meaning of our conclusions, is to build the so-called ``confidence belt'', in the plane $(x,\theta)$ , where x is the generic result of our experiment and $\theta$ is the parameter we want to estimate. There are several ways to construct a confidence belt. Here, we decided to follow the recipe given by Feldman and Cousins (Feldman & Cousins 1998). We require our confidence belt to have the property to guarantee a selected coverage over all the parameters region. This selected coverage is, for us, C = 0.9. In reality, our method is a little different from Feldman and Cousins' one (Feldman & Cousins 1998), because we also choose to set a small false alarm probability. This choice, in fact, causes the coverage to be more than the goal-coverage in the upper-limit region, namely for $x < x_{\rm FA}$ . This over-coverage is the price we need to pay in order to have a small false-alarm probability. The confidence belt that we find can be seen in Fig. 3. The construction of our confidence belt proceeds as follows. For each fixed value $\bar{x}$ , of x we define $\theta_{\rm best}$ to be the value of $\theta$ that maximizes the likelihood $f(x,\theta)$ , requiring the physical constraint that $\theta_{\rm best}\geq 0$ . In particular, if the measured x is less than its average value $\bar{x}$ , we impose $\theta_{\rm best} = 0$ , because if the result of the measurement is less than the mean value, the best estimator of the signal is 0. Now, for each possible value $\bar{\theta}$ , we calculate the likelihood ratio given by

$\begin{displaymath}R(\bar{x},\bar{\theta})= \frac{f(\bar{x},\bar{\theta})}{f(\bar{x},\theta_{\rm best})}\cdot \end{displaymath}$

(18)

This ratio of likelihoods is the function that we use to choose the confidence intervals. In fact, for each choice of $\bar{\theta}$ , the confidence interval (x₁,x₂) is uniquely defined by the requirements (19) and (20):

R(x₁)=R(x₂)

(19)

and

$\begin{displaymath}\int_{x_1}^{x_2} f(x,\bar{\theta}) = C. \end{displaymath}$

(20)

$\begin{figure} \par\includegraphics[width=10cm,clip]{aa11946-09f3.eps} \end{figure}$	Figure 3: The construction of the confidence intervals. Once we have decided confidence belt parameters (goal coverage and false alarm probability), the shape of the belt is the same for all the targets. However, due to different noise standard deviation and bin 0 value, the target sources have different abscissa.
Open with DEXTER

Table 3: Values of the statistical quantities for 3 targeted MSPs.

If the condition (19) cannot be satisfied within the condition (20), we consider as good the confidence interval also if R(x₂) < R(x₁). The choice of this ordering principle for the choice of the confidence intervals, will result in a more regular behavior of the confidence belt in the regions of the parameter space where x is very low. The pairs of values (x₁,x₂) are taken starting from the value $x_{\max}$ which maximizes $R(x,\bar{\theta})$ for a given $\bar{\theta}$ . Taking all the different values of $\theta$ , we cover all the parameter space and so we can trace the confidence belt. We start by selecting a grid of values in the parameters space $(x,\theta)$ and calculating the value of the likelihood ratio R over all the points of the grid. It's important to notice that, unfortunately, the problem of finding $\theta_{\rm best}$ can be solved only numerically, because in our case the probability density function f is very complicated: it is in fact expressed in terms of the regularized hypergeometric $_0\bar{F}_1$ functions, and the problem of finding an always viable relationship $\theta_{\rm best}=\theta_{\rm best}(f,x)$ is not analytically solvable. Then the program takes a value $\bar{\theta}$ and, for the section of the parameter grid at $\theta=\bar{\theta}$ , searches for $x_{\max}$ . Next, we move from $x_{\max}$ to larger values of x, and the program shows all the pairs of values (x₁,x₂) that best satisfy Eq. (19). For each pair, the program integrates the probability density function to find the coverage. Finally, between all the coverage values, the program extracts the pair for which the coverage is the most similar to C and not less than it, thus computing the confidence interval at $\bar{\theta}$ . The computed confidence belt is shown in picture 3. The confidence belt must be interpreted in this way: given a result $\bar{x}$ of the experiment, we trace a vertical line, which intercepts the edges of the confidence belt, thus giving the extreme values $\theta_1$ and $\theta_2$ of the confidence interval.

The final values of the statistical quantities are shown in Table 3. These results produce the upper limits summarized in Table 4 if we require the different coverages 0.68, 0.9 and 0.95.

Table 4: Measured upper limits for the 3 targeted MSPs for 3 different goal coverages.

5 Discussion

We have performed a search for GW emission from 3 MSPs using the resonant mass detector AURIGA.

It's important to notice here that this is the first attempt to look at GW emission from known binary pulsars using a resonant detector. Moreover, this is done by implementing a modified version of the Feldman and Cousins confidence belt method, which is quite a new feature for searches for continuous quasi-periodic signals. The upper limits calculated in Sect. 4 translate to limits on the neutron star ellipticity $\varepsilon$ . We adopted the distances calculated from the dispersion measure of the targets and a model for the distribution of the electrons in the interstellar medium (Taylor & Cordes 1993) in the cases of J0218+4232 (5.8 kpc) and PSR J1939+2134 (3.6 kpc), whereas for PSR J0024-7204J we used the distance of the related cluster 47 Tucanae (4.5 kpc, from feb 2003 revision of the Harris catalog, Harris 1996). According to the considerations in (Abbott et al. 2005) we find the following upper limits on $\varepsilon$ : $3.7 \times 10^{-5}$ for PSR J0024-7204J and $4.1 \times 10^{-5}$ for PSR J0218+4232. The inferred upper limit on $\varepsilon$ scales with the adopted distance d of the source like d^-1.

Unfortunately, given the current theoretical predictions for $\varepsilon$ (see e.g. Jones 2002; Ruderman 2006), these values are not very restrictive yet, reflecting the still relatively low sensitivity of the class of the detectors like AURIGA with respect to the interferometric GW detectors.

In particular, the most sensitive search to date performed by using interferometric data exploited the LIGO runs S3 and S4 (Abbott et al. 2007). For the 3 MSPs that we looked for, the resulting upper limits on h were: $h^{\rm u.l.}_{95 \%}= 7.41 \times 10^{-25}$ for J0024-7204J, $h^{\rm u.l.}_{95 \%}= 1.11 \times 10^{-24}$ for J0218+4232 and $h^{\rm u.l.}_{95 \%}= 1.65 \times 10^{-24}$ for J1939+2134. Comparing these values with the ones in Table 4, we see that the capabilities of the interferometers allows for reaching upper limits about 20 up to 50 times better than what a detector like AURIGA can do. The most recent improvement in this field is an all sky search using LIGO S5 data (Abbot et al. 2008).

However, the procedure introduced in this paper is suitable to be applied to any other more sensitive and broadband GW detectors, such as Advanced LIGO and Advanced VIRGO (Viceré 2007): firstly, the method is not affected by cumulative phase errors because we recalculate the phase every time a new frame begins, using the associated GPS time and the pulsar ephemeris. Secondly, the noise statistics is well known because we have a background given by many bins which don't carry the signal, and so the accuracy in the noise variance estimate is very satisfatory.

We finally note that a very promising application of the procedure described in this work involves the class of the so called Accreting X-ray millisecond pulsars (AMXPs). They are NSs undergoing transients phases of intense X-ray activities (known as outbursts) due to accretion of matter from the companion star in a LMXB system. During the outbursts, these sources display coherent pulsations in the X-ray band at frequency of order hundreds of Hz. These pulsations are believed to reflect the spin frequency of the neutron star. This opened the possibility of performing an accurate timing of the NSs, tracking their rotational phase for the duration (typically at least few weeks long) of the outburst (Falanga et al. 2005; Burderi et al. 2007). Since the X-ray emissions is due to matter falling onto the star from the accretion disk, it is possible that, during the outbursts, the ellipticity of the star changes. Moreover, also emitting mechanisms such the r-modes are likely to become instable and so available as GW emission channels. The amplitude of these effects is still uncertain (see e.g. (Watts et al. 2008), for limits and caveats); however a folding procedure like that described in this paper and based on ephemeris provided by X-ray data, could be a very effective method for revealing this putative GW emission.

References

Abbott, B., Abbott, R., Adhikari, R., et al. (the LIGO Scientific Collaboration) 2005, Phys. Rev. Lett., 94, 181103 (In the text)
Abbott, B., et al. (the LIGO Scientific Collaboration) 2007,[ArXiv:gr-qc/0702039v5] (In the text)
Abbott, B., et al. (the LIGO Scientific Collaboration) 2008,[arXiv:0810.0283v1] (In the text)
Alpar, M. A., Cheng, A. F., Ruderman, M. A., & Shaham, J. 1982, Nature, 300, 728 [NASA ADS] [CrossRef] (In the text)
Andersson, N., & Kokkotas, K. D. 2001, Int. J. Mod. Phys. D, 10, 381 [NASA ADS] [CrossRef] (In the text)
Andersson, N., Kokkotas, K. D., & Ferrari, V. 2001, Int. J. Mod. Phys. D, 10, 381 [NASA ADS] [CrossRef] (In the text)
Backer, D. C., Kulkarni, S. R., Heiles, C., Davis, M. M., & Goss, W. M. 1982, Nature, 300, 615 [NASA ADS] [CrossRef] (In the text)
Baggio, L., Bignotto, M., Bonaldi, M., et al. 2005, Phys. Rev. Lett., 94, 241101 [NASA ADS] [CrossRef] (In the text)
Bhattacharya, D., & van den Heuvel, E. P. J. 1991, Phys. Rep., 203, 1 [NASA ADS] [CrossRef] (In the text)
Bildsten, L. 1998, ApJ, 501, L89 [NASA ADS] [CrossRef] (In the text)
Burderi, L., Di Salvo, T., Lavagetto, G., et al. 2007, ApJ, 657, 961 [NASA ADS] [CrossRef]
Cognard, I., & Backer, D. 2004, ApJ, 612, 125 [NASA ADS] [CrossRef] (In the text)
Cutler, C., & Thorne, K. S. 2001, Proceedings of GR16 (In the text)
Falanga, M., di Salvo, T., Burderi, L., et al. 2005, A&A, 436, 313 [NASA ADS] [CrossRef] [EDP Sciences]
Feldman, G. J., & Cousins, R. D. 1998, Phys. Rev. D, 57, 3873 [NASA ADS] [CrossRef] (In the text)
Harris, W. E. 1996, AJ, 112, 1487 [NASA ADS] [CrossRef] (In the text)
Hessels, J. W. T., Ransom, S. M., Stairs, I. H., et al. 2006, Science, 311, 1901 [NASA ADS] [CrossRef] (In the text)
Jones, D. I. 2002, Classical and Quantum Gravity, 1255 (In the text)
Levin, Y., & Ushomirsky, G. 2001, MNRAS, 324, 917 [NASA ADS] [CrossRef]
Manchester, R. N., & Taylor, J. 1974, in Pulsars (Wiley: San Francisco) (In the text)
Manchester, R. N., Lyne, A. G., Robinson, C., et al. 1991, Nature, 352, 219 [NASA ADS] [CrossRef] (In the text)
Manchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, AJ, 129, 1993 [NASA ADS] [CrossRef]; http://www.atnf.csiro.au/research/pulsar/psrcat (In the text)
Misner, C. W., Thorne, K. S., & Weeler, J. A. 1973, Gravitation (W.H. Freeman, New York), 1023 (In the text)
Navarro, J., de Bruyn, G., Frail, D., Kulkarni, S. R., & Lyne, A. G. 1995, ApJ, 455, L55 [NASA ADS] [CrossRef] (In the text)
Owen, B. J. 1999, 3rd Edoardo Amaldi Conference on Gravitational Waves (In the text)
Possenti, A., et al. 2000, ApJS, 125, 463 [NASA ADS] [CrossRef] (In the text)
Ruderman, M. 2006, proceedings of the 363rd Heraeus Seminar on Neutron Stars and Pulsars
Taylor, J. H., & Cordes, J. M. 1993, ApJ, 411, 674 [NASA ADS] [CrossRef] (In the text)
Viceré, A., for the VIRGO Collaboration 2007, talk at conference PAFT (In the text)
Vinante, A. (for the AURIGA collaboration) 2006, Class. Quantum Grav., 23, S103-S110 (In the text)
Watts, A. L., Krishnan, B., Bildsten, L., & Schutz, B. F. 2008, MNRAS, 389, 893 [NASA ADS] [CrossRef] (In the text)
Weinberg, S. 1972, Gravitation and Cosmology: principles and applications of the general theory of relativity

Footnotes

...: For the AURIGA Collaboration.
... detector: See http://www.auriga.lnl.infn.it/ for a detailed description of the instrument).
... catalogue: http://www.atnf.csiro.au/research/pulsar/psrcat
... TEMPO: http://www.atnf.csiro.au/research/pulsar/timing/tempo

All Tables

Table 1: $\nu _0$ is the spin frequency of the targeted MSPs at the time of the AURIGA run, whereas $\nu _{\rm gw} {\rm [Hz]}$ is the frequency of the GW signal we searched for.

Table 2: Antenna Pattern factors.

Table 3: Values of the statistical quantities for 3 targeted MSPs.

Table 4: Measured upper limits for the 3 targeted MSPs for 3 different goal coverages.

All Figures

$\begin{figure} \par\includegraphics[width=14cm,clip]{aa11946-09f1.eps} \end{figure}$	Figure 1: Complete data-analysis pipeline. The h-reconstructed signal is passed trough a band-pass digital filter before the Doppler correction. This allows to decimate the samples to deal with a reasonable amount of data. Then, after the Fourier transformation, all the signal, if present, is in the bin 0 of the spectrum. The other bins are used to infer the statistics of the noise, thus allowing to calculate the confidence intervals for the GW amplitude.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=11.5cm,]{11946fg2.eps} \end{figure}$	Figure 2: The AURIGA strain sensitivity (gray) compared with the expected (black curve) sensitivity at $4.5\rm~ K$ .
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=10cm,clip]{aa11946-09f3.eps} \end{figure}$	Figure 3: The construction of the confidence intervals. Once we have decided confidence belt parameters (goal coverage and false alarm probability), the shape of the belt is the same for all the targets. However, due to different noise standard deviation and bin 0 value, the target sources have different abscissa.
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.