Detection and localization of single-source gravitational waves with pulsar timing arrays

Abstract

Pulsar timing arrays (PTAs) can be used to search for very low frequency (10⁻⁹–10⁻⁷ Hz) gravitational waves (GWs). In this paper, we present a general method for the detection and localization of single-source GWs using PTAs. We demonstrate the effectiveness of this new method for three types of signals: monochromatic waves as expected from individual supermassive binary black holes in circular orbits, GWs from eccentric binaries and GW bursts. We also test its implementation in realistic data sets that include effects such as uneven sampling and heterogeneous data spans and measurement precision. It is shown that our method, which works in the frequency domain, performs as well as published time-domain methods. In particular, we find it equivalent to the |$\mathcal {F}_{e}$|-statistic for monochromatic waves. We also discuss the construction of null streams – data streams that have null response to GWs, and the prospect of using null streams as a consistency check in the case of detected GW signals. Finally, we present sensitivities to individual supermassive binary black holes in eccentric orbits. We find that a monochromatic search that is designed for circular binaries can efficiently detect eccentric binaries with both high and low eccentricities, while a harmonic summing technique provides greater sensitivities only for binaries with moderate eccentricities.

gravitational waves, methods: data analysis, pulsars: general

1 INTRODUCTION

The direct detection of gravitational waves (GWs) will have profound implications in physics and astronomy. This may become possible within this decade. In the audio band (10–1000 Hz), second-generation km-scale laser interferometers, such as Advanced LIGO (Harry 2010), Advanced Virgo (Acernese et al. 2015) and KAGRA (Somiya 2012), are about to start scientific observations as early as the second half of 2015 and are expected to detect dozens of compact binary coalescence events per year when they achieve their design sensitivities in around 2020 (Abadie et al. 2010; Aasi et al. 2013). In the nanohertz frequency range (1–100 nHz), high-precision timing observations of millisecond pulsars (pulsar timing arrays – PTAs) provide a unique means of detecting GWs. With the concept proposed decades ago (Sazhin 1978; Detweiler 1979; Hellings & Downs 1983; Foster & Backer 1990), there are now three major PTA projects around the globe, namely, the Parkes Pulsar Timing Array (PPTA; Hobbs 2013; Manchester et al. 2013), the European Pulsar Timing Array (Kramer & Champion 2013), and NANOGrav (McLaughlin 2013). While these PTAs have individually collected high-quality data spanning ≳ 5 yr for ∼20 pulsars and produced some astrophysically interesting results (e.g. Shannon et al. 2013), they have also been combined to form the International Pulsar Timing Array (IPTA; Hobbs et al. 2010; Manchester 2013) aiming at significantly enhanced sensitivities.

The primary sources in the PTA frequency band are inspiralling supermassive binary black holes (SMBBHs). It is widely considered that a stochastic GW background due to the combined emission from a large number of individual SMBBHs over cosmological volume (e.g. Sesana 2013b; Ravi et al. 2014, for recent work) provides the most promising target; indeed, some studies suggested that a detection of this type could occur as early as 2016 (Siemens et al. 2013). Analyses of actual PTA data have previously focused on a search for such a GW background, leading to more and more stringent constraints on the fractional energy density of the background (Jenet et al. 2006; Yardley et al. 2011; van Haasteren et al. 2011; Demorest et al. 2013; Shannon et al. 2013). Over the past few years interest has grown substantially regarding the prospects of detecting single-source GWs using PTAs, for example, for individual SMBBHs (Sesana, Vecchio & Volonteri 2009; Lee et al. 2011; Mingarelli et al. 2012; Ravi et al. 2015), for GW memory effects associated with SMBBH mergers (Seto 2009; Cordes & Jenet 2012; Madison, Cordes & Chatterjee 2014), for GW bursts (Pitkin 2012) and for unanticipated sources (Cutler et al. 2014). In the meantime, many data analysis methods have been proposed in the context of PTAs for single-source GW detection, for example, for monochromatic GWs emitted by SMBBHs in circular orbits (Yardley et al. 2010; Babak & Sesana 2012; Ellis, Siemens & Creighton 2012; Ellis 2013; Wang, Mohanty & Jenet 2014; Taylor, Ellis & Gair 2014; Zhu et al. 2014), for GW memory effects (van Haasteren & Levin 2010; Wang et al. 2015), and GW bursts (Finn & Lommen 2010; Deng 2014). A few papers have presented searches in real PTA data for GWs from circular SMBBHs and yielded steadily improved upper limits on the GW strain amplitude (Yardley et al. 2010; Arzoumanian et al. 2014; Zhu et al. 2014). While circular binaries emit GWs at the second harmonic of the orbital frequency, eccentric binaries radiate at multiple harmonics. Jenet et al. (2004) first derived the expression of pulsar timing signals produced by eccentric binaries and developed a framework in which pulsar timing observations can be used to constrain properties of SMBBHs.

In our previous work (Zhu et al. 2014), we performed an all-sky search for GWs from SMBBHs in circular orbits using the latest PPTA data set. Here, we adapt the network analysis method used in the context of ground-based interferometers (e.g. Pai, Dhurandhar & Bose 2001; Wen & Schutz 2005, 2012; Klimenko et al. 2008; Wen 2008; Sutton et al. 2010) as a general method for detection and localization of single-source GWs using PTAs. In particular, we consider the following types of sources: (1) SMBBHs in circular orbits; (2) eccentric SMBBHs; and (3) GW bursts. We demonstrate the effectiveness of this method using synthetic data sets that contain both idealized and realistic observations.

The organization of this paper is as follows. In Section 2, we briefly review the features of pulsar timing signals caused by single-source GWs and discuss the signal models used in this work. In Section 3, we present the mathematical framework of our method and propose practical detection statistics. The relation to the method used in our previous work is also discussed. Using idealized simulations, we show examples of detection, localization and waveform estimation in Section 4. We demonstrate the implementation of the method in realistic data sets in Section 5. In Section 6, we present sensitivities to eccentric SMBBHs. In particular, we compare two detection strategies towards the detection of eccentric binaries – a monochromatic search and a harmonic summing technique. Finally, we summarize and outline future work in Section 7.

2 SINGLE-SOURCE GWs AND PULSAR TIMING

In pulsar timing, the times of arrival (ToAs) of radio pulses from millisecond pulsars are measured and compared with predictions based on a timing model that describes the pulsar's rotation (e.g. spin period and spin-down rate), the relative geometry between the pulsar and the observer, and other effects such as dispersion of radio waves due to interstellar medium (Lorimer & Kramer 2005; Edwards, Hobbs & Manchester 2006). The differences between the measurements and predictions are called timing residuals. These residuals may be partly attributed to effects of GWs as they are not normally included in a timing model. Observations of a single pulsar cannot make a definite detection of a GW, but can lead to constraints on the strength of potential GWs (see e.g. Yi et al. 2014, for a recent work). By timing an array of pulsars, GWs can be unambiguously detected by searching for correlated signatures among different pulsars. Typical PTA observations have a sampling interval of weeks and span over ∼10 yr, implying a sensitive GW frequency range of ∼1–100 nHz.

In this work, we are interested in single-source GWs, i.e. those coming from some particular directions in the sky. Timing residuals of such an origin can be generally written as

\begin{equation} r(t,\hat{\Omega }) = F_{+}(\hat{\Omega })\Delta A_{+}(t) + F_{\times }(\hat{\Omega }) \Delta A_{\times }(t), \end{equation}

(1)

where |$\hat{\Omega }$| is a unit vector pointing from the GW source to the observer. The two functions |$F_{+}(\hat{\Omega })$| and |$F_{\times }(\hat{\Omega })$| are the geometric factors as given by (e.g. Lee et al. 2011)

\begin{eqnarray} F_{+}(\hat{\Omega }) &=& \frac{1}{4(1-\cos \theta )}\left\lbrace (1+\sin ^2 \delta )\cos ^2 \delta _p \cos [2(\alpha -\alpha _p)] \right.\nonumber\\ &&\left.- \sin 2\delta \sin 2\delta _p\cos (\alpha -\alpha _p) + \cos ^2 \delta (2-3\cos ^2 \delta _p)\right\rbrace \nonumber\\ && \\ \nonumber \end{eqnarray}

(2)

\begin{eqnarray} F_{\times }(\hat{\Omega }) &=& \frac{1}{2(1-\cos \theta )}\lbrace \cos \delta \sin 2\delta _p \sin (\alpha -\alpha _p) \nonumber\\ & & - \sin \delta \cos ^2 \delta _p \sin [2(\alpha -\alpha _p)]\rbrace , \end{eqnarray}

(3)

where cos θ = cos δcos δ_pcos (α − α_p) + sin δsin δ_p with θ being the opening angle between the GW source and pulsar with respect to the observer, α (α_p) and δ (δ_p) are the right ascension and declination of the GW source (pulsar), respectively. Note that we separate the polarization angle from |$F_{+}(\hat{\Omega })$| and |$F_{\times }(\hat{\Omega })$| and put it into ΔA_{+, ×}(t) for convenience of later analyses. The two functions |$F_{+}(\hat{\Omega })$| and |$F_{\times }(\hat{\Omega })$| are analagous to the antenna pattern functions used in the context of laser interferometric GW detection (Thorne 1987).

Because the wavelengths of GWs in the PTA band are much smaller than the pulsar-Earth distance, GW induced timing residuals are the combination of two terms – the Earth term A_{+, ×}(t) and the pulsar term A_{+, ×}(t_p) (see e.g. Jenet et al. 2004):

\begin{eqnarray} \Delta A_{+,\times }(t) &=& A_{+,\times }(t)-A_{+,\times }(t_{\rm p})\\ && \nonumber \end{eqnarray}

(4)

\begin{equation} t_{\rm p} = t-d_{\rm p}(1-\cos \theta )/c, \end{equation}

(5)

where d_p is the pulsar distance and we have adopted the plane wave approximation. A₊(t) and A_×(t) are source-dependent functions. Throughout this paper, we only consider the correlated Earth-term signals. In fact for the case of continuous waves as expected from SMBBHs, pulsar terms will act as an extra source of uncorrelated noise for different pulsars. For GW bursts, whose duration is smaller than the data span, it is very unlikely that pulsar terms and Earth terms are simultaneously present in timing residuals for one particular pulsar unless the source sky direction is very close to that pulsar. Below we briefly discuss the signal properties for the three types of sources considered in this paper.

At the leading Newtonian order SMBBHs in circular orbits are expected to generate pulsar timing signals with the following forms (see e.g. Zhu et al. 2014):

\begin{eqnarray} A_{+}(t) &=& \frac{h_0}{2\pi f} \left[(1+\cos ^{2}\iota ) \cos 2\psi \sin (2\pi f t+\phi _0)\right. \nonumber\\ &&\left.+\,2\cos \iota \sin 2\psi \cos (2\pi f t+\phi _0)\right] \end{eqnarray}

(6)

\begin{eqnarray} A_{\times }(t) &=& \frac{h_0}{2\pi f} \left[(1+\cos ^{2}\iota ) \sin 2\psi \sin (2\pi f t+\phi _0)\right.\nonumber\\ &&\left.-\;2\cos \iota \cos 2\psi \cos (2\pi f t+\phi _0)\right], \end{eqnarray}

(7)

where the GW strain amplitude is |$h_0=2 (G M_{\rm c})^{5/3} (\pi f)^{2/3} c^{-4} d_{\rm L}^{-1}$| with d_L being the luminosity distance of the source and M_c being the chirp mass defined as |$M_{\rm c}^{5/3} = m_1 m_2(m_1+m_2)^{-1/3}$| where m₁ and m₂ are the binary component masses, f is the GW frequency, ι is the inclination angle of the binary orbit, ψ is the GW polarization angle, ϕ₀ is a phase constant. Note that we have neglected the frequency evolution over the observation span (typically ∼10 yr). This represents the majority of circular binaries that are observable for PTAs as shown in Sesana & Vecchio (2010).

Although it is well known that radiation of GWs circularizes the binary orbits (Peters & Mathews 1963), the assumption of circular orbits is not always appropriate. Recent models for the SMBBH population including the effects of binary environments on orbital evolution suggest that eccentricity is important for GW frequencies ≲ 10⁻⁸ Hz (Sesana 2013a; Ravi et al. 2014). In this work, we use the expressions of GW-induced timing residuals for eccentric SMBBHs as given in Jenet et al. (2004). It is interesting to note that eccentric binaries emit GWs at multiple harmonics of the binary orbital frequency. At low eccentricities the emission is dominated by the second harmonic, while for high eccentricities the orbital frequency itself will dominate.

Generally a GW burst is defined as a transient signal with a duration smaller than the observation span. This may be the only information we know about the source. For this reason, we use a simple but general sine-Gaussian model for timing residuals induced by GW bursts:

\begin{eqnarray} A_{+}(t) &=& A\exp \left(-\frac{(t-t_0)^2}{2\tau ^2}\right) \nonumber\\ &&\left\lbrace (1+\cos ^{2}\iota )\cos 2\psi \cos [2\pi f_{0} (t-t_0)+\phi _0] \right.\nonumber\\ &&\left.-\,2\cos \iota \sin 2\psi \sin [2\pi f_{0} (t-t_0)+\phi _0]\right\rbrace \end{eqnarray}

(8)

\begin{eqnarray} A_{\times }(t) &=& A\exp \left(-\frac{(t-t_0)^2}{2\tau ^2}\right) \nonumber\\ &&\left\lbrace (1+\cos ^{2}\iota )\sin 2\psi \cos [2\pi f_{0} (t-t_0)+\phi _0] \right.\nonumber\\ &&\left.+\,2\cos \iota \cos 2\psi \sin [2\pi f_{0} (t-t_0)+\phi _0]\right\rbrace , \end{eqnarray}

(9)

where A is the signal amplitude (in seconds), τ is the Gaussian width, ι is the source inclination angle, f₀ is the central frequency, t₀ and ϕ₀ is the time and phase at the mid-point of the burst, respectively. The sine-Gaussian model used here can represent qualitatively the signals for parabolic encounters of two massive black holes as studied in Finn & Lommen (2010) and more recently in Deng (2014). The detection method that we will describe in the next section should apply equally well to other burst sources such as cosmic (super)string cusps and kinks (e.g. Vilenkin 1981; Damour & Vilenkin 2000; Siemens et al. 2006) and triplets of supermassive black holes (Amaro-Seoane et al. 2010).

3 THE DATA ANALYSIS METHOD

In this section, we describe how the singular value decomposition (SVD) can be used in a general method for the detection and localization of single-source GWs using PTAs. The method is adapted from the coherent network analysis method used in the context of ground-based GW interferometers. There are some important features that are unique to PTA observations, e.g. (1) PTA data are irregularly sampled in contrast to nearly continuous sampling for ground-based experiments; and (2) a least-squares fitting process is performed to obtain estimates of timing parameters such as the pulsar's spin period and its first time derivative, pulsar position and proper motion, etc. As we show later in this section, the SVD method proposed here relies on transforming the timing residuals of each pulsar to the frequency domain. For the idealized simulations (assuming even sampling and white Gaussian noise with equal error bars) used in Section 4, a discrete Fourier transform was used. In Section 5 for more realistic data sets, we adopt a maximum-likelihood-based method to estimate Fourier components of the timing residuals making use of the noise covariance matrix (e.g. section 3 in Zhu et al. 2014). This is equivalent to the least-squares spectral analysis method (see e.g. Coles et al. 2011). Our method works with post-fit timing residuals, i.e. after fitting ToAs for timing models of individual pulsars. The effects of the fitting process on our results will be discussed in Section 5.

For a given source direction, the timing residuals from an array of N_p pulsars can be generally written in the frequency domain as

\begin{equation} {\bf d}_{\rm k}={\bf F}_{\rm k}{\bf A}_{\rm k}+{\bf n}_{\rm k}, \end{equation}

(10)

where the index k denotes the kth frequency bin, |${\bf d}_{\rm k}$| are timing residual data, |${\bf F}_{\rm k}$| is the response matrix, |${\bf A}_{\rm k}$| are GW waveforms and |${\bf n}_{\rm k}$| is the timing noise. The data are whitened¹ so that

\begin{equation} {{\bf d}}_{\rm k}= {\left[\begin{array}{c}d_{1k}/\sigma _{1k}\\ d_{2k}/\sigma _{2k}\\ \vdots \\ d_{N_pk}/\sigma _{N_pk}\end{array}\right]}, {{\bf A}}_{\rm k} = {\left[\begin{array}{c}A_{+k}\\ A_{\times k}\end{array}\right]}, {{\bf n}}_{\rm k}={\left[\begin{array}{c}n_{1k}/\sigma _{1k}\\ n_{2k}/\sigma _{2k}\\ \vdots \\ n_{N_pk}/\sigma _{N_pk}\end{array}\right]}, \end{equation}

(11)

where |$\sigma ^2_{ik}$| is the noise variance of the ith pulsar at the kth frequency bin. Here |${\bf d}_{\rm k}$|⁠, |${\bf A}_{\rm k}$| and |${\bf n}_{\rm k}$| are all complex vectors, while the real whitened response matrix |${\bf F}_{\rm k}$| is defined as

\begin{equation} {{\bf F}}_{\rm k} = {\left[\begin{array}{cc}F^+_1/\sigma _{1k} &\quad F^\times _1 /\sigma _{1k} \\ F^+_2/\sigma _{2k} &\quad F^\times _2/\sigma _{2k} \\ \vdots &\quad \vdots \\ F^+_{N_p}/\sigma _{N_p k} &\quad F^\times _{N_p} /\sigma _{N_p k} \end{array}\right]}. \end{equation}

(12)

For simplicity, we will hereafter suppress the index k while keeping in mind that equations (10)–(12) apply to each frequency bin in the analysis. Data for all frequencies can be stacked (e.g. in order of increasing frequency) to preserve the format of equation (10), in which case |${\bf F}$| is a block-diagonal matrix (with a dimension of N_pN_k × 2N_k where N_k is the number of frequency bins) of |${\bf F}_{\rm k}$|⁠.

The SVD of the response matrix |${\bf F}$| can be written as

\begin{equation} {{\bf F}}={{\bf U}}{{\bf S}}{{\bf V}}^{\ast }, {{\bf S}} = {\left[\begin{array}{cc}s_{1} &\quad 0 \\ 0 &\quad s_{2} \\ 0 &\quad 0 \\ \vdots &\quad \vdots \\ 0 &\quad 0 \end{array}\right]}, \end{equation}

(13)

where |${\bf U}$| and |${\bf V}$| are unitary matrices with dimensions of N_p × N_p and 2 × 2, respectively, and the symbol * denotes the conjugate transpose. Singular values in |${\bf S}$| are ranked such that s₁ ≥ s₂ ≥ 0. We can then construct new data streams as follows:

\begin{equation} {\tilde{\bf d}}={\bf U}^{\ast }{\bf d}, \,\, {\tilde{\bf A}}={\bf V}^{\ast } {\bf A}, \,\, {\tilde{\bf n}}={\bf U}^{\ast }{\bf n}. \end{equation}

(14)

One can find that |${\tilde{\bf d}}= {\bf S}{\tilde{\bf A}}+{\tilde{\bf n}}$|⁠, and explicitly

\begin{equation} {\tilde{\bf d}}={\left[\begin{array}{c}s_{1}\left({{\bf V}}^{\ast } {{\bf A}}\right)_{1}+{\tilde{\bf n}}_{1} \\ s_{2}\left({{\bf V}}^{\ast } {{\bf A}}\right)_{2}+{\tilde{{\bf n}}}_{2} \\ {\tilde{{\bf n}}}_{3} \\ \vdots \\ {\tilde{{\bf n}}}_{N_{\rm p}}\end{array}\right]}. \end{equation}

(15)

In the absence of GW signals, the real and imaginary parts for each element of |${\tilde{\bf d}}$| independently follow the standard Gaussian distribution. Here, |${\tilde{\bf d}}_1$| and |${\tilde{\bf d}}_2$| are referred to the two signal streams² since they contain all information about GW signals if present, while the remaining terms are ‘null streams’ (denoted by |${\tilde{\bf d}}_{\rm {null}}$|⁠) since they have a null response to GWs. It has been shown that, in the context of GW detection with ground-based laser interferometers, null streams can be used as a consistency check on whether a candidate GW event is produced by detector noise or by a real GW (Wen & Schutz 2005) and ‘semi-null streams’ (i.e. |${\tilde{\bf d}}_2$| if s₁ ≫ s₂) can be included to improve the angular resolution (Wen 2008; Wen, Fan & Chen 2008). In this paper, we only consider Earth terms in our signal model so the null streams (as constructed in the current way) are indeed ‘null’ but in reality they would have some response to the pulsar-term signals. In a future study pulsar terms will be incorporated in the detection framework with the addition of N_p free parameters for pulsar distances in the response matrix |${\bf F}$|⁠.

It is straightforward to show that the maximum likelihood estimator for the GW waveform is

\begin{eqnarray} \hat{{{\bf A}}}=\bar{{{\bf F}}}{{\bf d}},\, {\rm {where}} \,\, \bar{{{\bf F}}}={{\bf V}}\bar{{{\bf S}}}{{\bf U}}^{\ast }, \,\bar{{{\bf S}}} = {\left[\begin{array}{lllll}1/s_{1} &\quad 0 &\quad 0 &\quad \cdots &\quad 0 \\ 0 &\quad 1/s_{2} &\quad 0 &\quad \cdots &\quad 0 \end{array}\right]}.\nonumber\\ \end{eqnarray}

(16)

The matrix |$\bar{{\bf F}}$| is the Moore-Penrose pseudo-inverse of the response matrix |${\bf F}$|⁠. The covariance matrix for the estimated waveform is

\begin{equation} {\rm {var}}(\hat{{\bf A}})= \left({\bf V}{\bf S}^{{\rm T}}{\bf S}{\bf V}^{\ast }\right)^{-1}. \end{equation}

(17)

It is interesting to note that the statistical uncertainties of the estimated waveforms are a linear combination of |$s_{1}^{-2}$| and |$s_{2}^{-2}$|⁠. Estimation of physical parameters can be obtained based on the estimated waveforms |$\hat{{\bf A}}$| and we leave this to a future study.

Now we propose our detection statistics for the three types of signals discussed in the previous section. For monochromatic GWs, the detection statistic can be written as

\begin{equation} \mathcal {P}_{\rm {mon}}= |{\tilde{\bf d}}_1|^{2}+|{\tilde{\bf d}}_2|^{2}. \end{equation}

(18)

It is important to note that: (1) in the absence of GW signals, |$\mathcal {P}_{\rm {mon}}$| follows a χ² distribution with 4 degrees of freedom; and (2) the detection statistic applies to a given GW frequency and source sky location. In practice when such information is unknown, a search is usually performed to find the maximum statistic.

For signals produced by eccentric binaries, we use a harmonic summing technique for which the detection statistic is given by

\begin{equation} \mathcal {P}_{\rm {ecc}}= \sum ^{N_{\rm h}}_{j=1} \left(|{\tilde{\bf d}}_{1,jk_{0}}|^{2}+|{\tilde{\bf d}}_{2,jk_{0}}|^{2}\right), \end{equation}

(19)

where N_h corresponds to the highest harmonics considered in the search, k₀ is the bin number for the binary orbital frequency. In the absence of GW signals, |$\mathcal {P}_{\rm {ecc}}$| follows a χ² distribution with 4N_h degrees of freedom. In practice N_h should be determined as the one that gives the lowest false alarm probability (FAP). When the orbital frequency is unknown one should search over all possible frequencies to find the maximum statistic.

For GW bursts with unknown waveforms, we adopt a time-frequency strategy in which the PTA data are first divided into small segments and then SVD is applied to each segment to output the two signal streams. Note that it is usually necessary to allow overlap between successive segments to catch signals occurring in the beginning or end of the segment. The detection of GW bursts of unknown waveforms generally involves searching for any ‘tracks’ or ‘clusters’ of excess power in the time-frequency space (e.g. Anderson & Balasubramanian 1999; Wen & Gair 2005). So the detection statistic is designed as accumulating signal power in (a ‘box’ of) the time-frequency domain and can be written as

\begin{eqnarray} \mathcal {P}_{\rm {GWb}}(i,k)= \sum ^{l/2}_{a=-l/2}\sum ^{m/2}_{b=-m/2} \left(|{\tilde{\bf d}}_{1,(i+a)(k+b)}|^{2}+|{\tilde{\bf d}}_{2,(i+a)(k+b)}|^{2}\right), \nonumber\\ \end{eqnarray}

(20)

for the ith segment and kth frequency bin. Here, l and m is the length of box in time and frequency, respectively. If the data consist of only Gaussian noise, |$\mathcal {P}_{\rm {GWb}}$| follows a χ² distribution with 4N_b degrees of freedom (where N_b = l × m).

The detection statistics proposed here all follow a non-central χ² distribution with their corresponding degrees of freedom when signals are present in the data. It is convenient to define the expected signal-to-noise ratio (ρ) as the non-centrality parameter

\begin{equation} \langle \mathcal {P}\rangle = N_{\rm {dof}}+\rho ^{2}, \end{equation}

(21)

where the brackets 〈...〉 denote the ensemble average of the random noise process, N_dof is the number of degrees of freedom for the ‘central’ distribution of |$\mathcal {P}$|⁠. We use ρ to quantify the signal strength in our simulations. Note that ρ² equals the detection statistic calculated for noiseless data. In a frequentist detection framework, we are interested in the FAP of a measured detection statistic |$\mathcal {P}_{0}$|⁠, i.e. the probability that |$\mathcal {P}$| exceeds the measured value for noise-only data. For the aforementioned methods the single-trial FAP is given by |$1-\mathrm{CDF}(\mathcal {P}_{0};\chi ^{2})$| where CDF( ; χ²) denotes the cumulative distribution function (CDF) for the central χ² distribution in question. When a search is performed over unknown source parameters, the total FAP is |$1-[\mathrm{CDF}(\mathcal {P}_{\rm {max}};\chi ^{2})]^{N_{\rm {trial}}}$| for the maximum detection statistic |$\mathcal {P}_{\rm {max}}$| found in the search with N_trial being the trials factor, which is defined as the number of independent cells in the searched parameter space for a grid-based search.

3.1 Relation to the ‘A₊A_×’ method

In the pulsar timing software package tempo2 (Hobbs, Edwards & Manchester 2006), there exists a functionality with which two time series |$A_{+,\times }^{{\rm {t2}}}(t)$| can be estimated for a given sky direction. These two time series correspond to two polarizations of the coherent Earth-term timing residuals. Here, we call it the ‘A₊A_×’ method. In this method |$A_{+,\times }^{{\rm {t2}}}(t)$| are modelled as time-varying parameters and treated just like any other normal parameters in a timing model so that they can be estimated through a global least-squares fitting routine. As observations for different pulsars are usually unevenly sampled and not at identical times, linear interpolation is used for such a global fit. The ‘A₊A_×’ method enables one to simultaneously fit for single-source GWs and normal pulsar timing parameters. To avoid the covariance between the global fit for |$A_{+,\times }^{{\rm {t2}}}(t)$| and the fit for timing model parameters of individual pulsars, some constraints must be set on the two time series. Currently three kinds of constraints are implemented in tempo2, namely, (1) quadratic constraints that correspond to pulsar spin parameters, (2) annual sinusoids for pulsar positions and proper motions, and (3) biannual sinusoids for pulsar parallax. These were first introduced and implemented in Keith et al. (2013) when correcting the dispersion measure variations for individual pulsars.

The ‘A₊A_×’ method was first illustrated for arbitrary GW bursts in fig. 5 of Hobbs (2013). It was used in Zhu et al. (2014) for an all-sky search for monochromatic GWs in the latest PPTA data set. A complete presentation on the ‘A₊A_×’ method and its applications to single-source GW detection will be provided in a forthcoming paper (Madison et al., in preparation). Here, we briefly discuss the relation between the ‘A₊A_×’ method and the SVD method proposed in this paper. First, the principle of both methods is identical since both use the response matrix |${\bf F}$| in a similar way and equation (16) essentially gives the least-squares solution to the two polarization waveforms, i.e. |$\hat{{\bf A}}$| in equation (16) is the frequency-domain equivalent of |$A_{+,\times }^{{\rm {t2}}}(t)$|⁠. The critical difference is in the implementation – the ‘A₊A_×’ method works in the time domain whereas the SVD method works in the frequency domain.

There are two advantages of the SVD method over the ‘A₊A_×’ method.

The SVD method is much faster when searching over the whole sky is required. This is because it works with post-fit residuals after data for each pulsar have been fitted for a timing model and a search over unknown source sky location only involves doing the SVD for the response matrix |${\bf F}$|⁠, while in the ‘A₊A_×’ method a global fit is done for each searched sky location (as in the current implementation).
A truncated SVD may be used to improve the detection sensitivity, which applies to the case when s₁ ≫ s₂, e.g. when the PTA has very low response to one of the two GW polarizations for some sky regions. This is possible especially for current PTAs whose sensitivities are dominated by a few best-timed pulsars.

The ‘A₊A_×’ method has been fully tested and implemented in real data for continuous GWs in Zhu et al. (2014). We will demonstrate the effectiveness of the SVD method with some examples using idealized data sets in Section 4 and then more realistic data sets in Section 5.

4 EXAMPLES USING IDEALIZED DATA SETS

For the purpose of illustration of our method, we consider an idealized PTA consisting of 20 PPTA pulsars. The simulated observations are evenly sampled once every two weeks with a time span of 10 yr. All observations contain stationary white Gaussian noise with an rms of 100 ns and equal error bars. The simulated data sets are produced as a combination of realizations of white Gaussian noise and Earth-term timing residuals. When we apply our detection statistics to noise-only data sets, it is confirmed that they follow a χ² distribution with their respective numbers of degrees of freedom.

In the following three subsections, we illustrate how our method works in terms of detection, localization and waveform estimation for the three types of signals considered in this work. The analysis is simplified as follows: (1) first the detection problem is demonstrated by evaluating the detection statistics at the injected source location; (2) then detection statistics are computed on a uniform grid of sky positions and for a range of frequencies to find the maximum; and (3) finally the frequency-domain waveforms A₊(f) and A_×(f) are estimated at the actual source sky location. To show the correctness of the estimation process, a number of noise realizations were performed in order to obtain average estimates of A_{+, ×}(f) which were then compared with the true waveforms. For simplicity, in the relevant figures we only plot the absolute values of A₊(f) and A_×(f), which are called spectral signatures. For the detection problem, we calculate the single-trial FAP to quantify the detection significance, which is only appropriate when source parameters (such as sky location and frequency) are known as we assume for the following examples. The simulated signals are weak to moderately strong with signal-to-noise ratios ranging from 5 to 10.

4.1 Monochromatic waves

For monochromatic waves, the simulated signal is characterized with the following parameters: h₀ = 1.16 × 10⁻¹⁵, f = 10 nHz, cos ι = 1, ψ = 0, ϕ₀ = 0, (α, δ) = (0, 0). The expected signal-to-noise ratio is ρ = 10. Fig. 1 shows the detection statistics (⁠|$\mathcal {P}_{\rm {mon}}$|⁠) as a function of frequency evaluated at the injected sky location. The maximum statistic, which gives extremely strong evidence of a detection, is found at a frequency of 9.94 nHz. To localize this source, we calculate |$\mathcal {P}_{\rm {mon}}$| for the same range of frequencies and for a grid of sky directions. At each sky direction, the maximum statistic over frequencies is recorded. Fig. 2 shows that the source is successfully localized to where the signal is generated.

$Detection statistics ($\mathcal {P}_{\rm {mon}}$) as a function of frequency for a simulated data set that includes a strong monochromatic signal (solid; solid black in online) and for the noise-only data set (dash; dashed red in online). The vertical line marks the injected frequency (10 nHz), while the horizontal line corresponds to a single-trial FAP of 10−4.$

Figure 1.

Detection statistics (⁠|$\mathcal {P}_{\rm {mon}}$|⁠) as a function of frequency for a simulated data set that includes a strong monochromatic signal (solid; solid black in online) and for the noise-only data set (dash; dashed red in online). The vertical line marks the injected frequency (10 nHz), while the horizontal line corresponds to a single-trial FAP of 10⁻⁴.

Open in new tab Download slide

Figure 2.

Sky map of detection statistics calculated for a simulated data set that includes a strong monochromatic signal. The signal is simulated at the centre of the map and the maximum detection statistic is found at ‘○’ (which is the same as the actual source location). Sky locations of the 20 PPTA pulsars are marked with ‘★’.

Open in new tab Download slide

Having already detected and localized the source, we use equations (16) and (17) to infer the GW waveforms. Fig. 3 shows the estimated spectral signatures along with the true spectra, indicating a very good estimation as one would expect since the injected signal is very strong. To check the correctness of our method, we overplot in Fig. 3 the average estimates taken over 100 noise realizations – they are almost identical to the true spectra.

Figure 3.

Estimated spectral signatures (thin solid; thin solid black in online) for the same data set (white Gaussian noise + a monochromatic GW) used in Figs 1 and 2, along with true spectra (solid; solid blue in online) and average estimates taken over 100 noise realizations (dash; dashed red in online).

Open in new tab Download slide

4.2 Eccentric binaries

The simulated source for this example is an eccentric SMBBH located in the Virgo cluster as described by the following parameters: m₁ = m₂ = 2 × 10⁸M_⊙, d_L = 16.5 Mpc, cos ι = 1, (α, δ) = (3.2594, 0.2219), an orbital frequency of 5 nHz, an initial orbital phase of 1 and an eccentricity of 0.5. The signal is simulated using equations give in section 2 of Jenet et al. (2004). In order to ‘detect’ this simulated signal, we first apply our analysis at the injected sky location and experiment on the number of harmonics that should be considered. It turns out that a harmonic summation up to the second harmonic gives the lowest FAP. We will further discuss in Section 6 the number of harmonics that should be included for binaries with different eccentricities. Fig. 4 shows the detection statistics (⁠|$\mathcal {P}_{\rm {ecc}}$|⁠) as a function of orbital frequency. The maximum statistic 59.85, corresponding to a FAP of 5 × 10⁻¹⁰, is found at 4.87 nHz. The expected signal-to-noise ratio for this injection is ρ = 7 when we perform a harmonic summing up to the second harmonic. Then we show in Fig. 5 the detection statistics calculated for the whole sky (after maximizing over orbital frequencies for each sky direction). The maximum statistic is found at a grid point close to the injected sky location.

$Detection statistics ($\mathcal {P}_{\rm {ecc}}$) as a function of orbital frequency for a simulated data set that includes a moderately strong signal (solid; solid black in online) and for the noise-only data set (dash; dashed red in online). The signal is produced by an eccentric SMBBH with an eccentricity of 0.5. Detection statistics are calculated as a harmonic summation up to the second harmonic. The vertical line marks the injected orbital frequency (5 nHz), while the horizontal line corresponds to a single-trial FAP of 10−4.$

Figure 4.

Detection statistics (⁠|$\mathcal {P}_{\rm {ecc}}$|⁠) as a function of orbital frequency for a simulated data set that includes a moderately strong signal (solid; solid black in online) and for the noise-only data set (dash; dashed red in online). The signal is produced by an eccentric SMBBH with an eccentricity of 0.5. Detection statistics are calculated as a harmonic summation up to the second harmonic. The vertical line marks the injected orbital frequency (5 nHz), while the horizontal line corresponds to a single-trial FAP of 10⁻⁴.

Open in new tab Download slide

Figure 5.

Sky map of detection statistics calculated for a simulated data set that includes a simulated signal produced by an eccentric SMBBH located in the Virgo cluster (marked by a white diamond). The maximum statistic is found at ‘○’. Sky locations of the 20 PPTA pulsars are marked with ‘★’.

Open in new tab Download slide

Fig. 6 shows the estimated spectral signatures assuming that we know the actual source location. Since the detection statistic varies only slightly within an area of tens of square degrees (as shown in Fig. 5), similar results should be obtained if we apply the spectral estimation analysis in the sky location where the maximum statistic is found. It is shown that reasonably good estimates are obtained only for the first two harmonics.

Figure 6.

Estimated spectral signatures (thin solid; thin solid black in online) for the same data set (white Gaussian noise + a signal expected from an eccentric SMBBH) used in Figs 4 and 5, along with true spectra (solid; solid blue in online) and average estimates taken over 100 noise realizations (dash; dashed red in online).

Open in new tab Download slide

4.3 GW bursts

In this example, the simulated GW burst signal is described by the following parameters: A = 100 ns, τ = 60 d, cos ι = 0.5, f₀=50 nHz, ψ = 1.57, ϕ₀ = 0, (α, δ) = (0.9512, − 0.6187), i.e. originating from the Fornax cluster. The signal occurs in the middle of our observations (MJD 55250). Waveforms of A_{+, ×}(t) are shown in the upper panel of Fig. 7. For this burst source the amplitudes of timing residuals are ≲ 100 ns, which means the signal should not be visibly apparent in individual pulsar data set.

$Upper panel: simulated timing signals from a GW burst. Lower panel: detection statistics $\mathcal {P}_{\rm {GWb}}$ (defined in equation 20) calculated in the time-frequency domain.$

Figure 7.

Upper panel: simulated timing signals from a GW burst. Lower panel: detection statistics |$\mathcal {P}_{\rm {GWb}}$| (defined in equation 20) calculated in the time-frequency domain.

Open in new tab Download slide

To dig out this burst signal from noisy data, we divide the 10-yr data set into segments of length of 300 d and calculate the detection statistic given by equation (20) in the time-frequency domain. The segment length is (roughly) chosen based on the knowledge of f₀, i.e. having ≳ one cycle per segment. In practice, since the time-frequency analysis for PTA data should not be limited by computational power, many trials of segment length can be performed to search for signals of different durations. The lower panel of Fig. 7 shows results of the time-frequency analysis assuming known source sky location. The most significant statistic 38.36, corresponding to a FAP of 9 × 10⁻⁸, is found at MJD 55257 and a frequency of 51.67 nHz. Therefore our analysis has clearly detected the injected signal and correctly identified its occurrence time and central frequency. For this example, the signal is very ‘isolated’ in the time-frequency space, so it is sufficient to only look at the maximum statistic in the time-frequency map (in this case, the signal-to-noise ratio is ρ = 5).

Fig. 8 shows the detection statistics evaluated at the whole sky (after maximizing over frequencies for each sky direction). The maximum statistic is found at a grid point close to the injected source location. Fig. 9 shows the inferred spectral signatures compared against the true spectra. As the injected signal is relatively weak for this example, the spectrum is not recovered as well as the previous two examples. For both plots, only the central segment that contains the majority of signal power (as illustrated in Fig. 7) was used.

Figure 8.

Sky map of detection statistics calculated for a simulated data set that includes a simulated GW burst originating from the Fornax cluster (marked by a white diamond). The maximum statistic is found at ‘○’. Sky locations of the 20 PPTA pulsars are marked with ‘★’.

Open in new tab Download slide

Figure 9.

Estimated spectral signatures (thin solid; thin solid black in online) for the same data set (white Gaussian noise + a GW burst) used in Figs 7 and 8, along with true spectra (solid; solid blue in online) and average estimates taken over 100 noise realizations (dash; dashed red in online).

Open in new tab Download slide

5 MORE REALISTIC DATA SETS

Examples given in the previous section have assumed idealized PTA observations that are evenly sampled³ and contain only stationary white Gaussian noise with equal error bars. Realistic features typical to real PTA data include: (1) observations are irregularly sampled and have varying ToA error bars. It is also common that the data span varies significantly among different pulsars; and (2) low-frequency (‘red’) timing noise may be present for some pulsars. To simulate realistic data sets, we make use of the actual data spans, sampling and error bars of the PPTA 6-yr Data Release 1 (DR1) data set that was published in Manchester et al. (2013) and used in Zhu et al. (2014) for an all-sky search for continuous GWs. The PPTA DR1 data set is publicly available online at a permanent link.⁴ As a check, in some simulations an uncorrelated red-noise process with a power-law spectrum is also included. This has no effect on the results since we assume the noise spectrum is known and thus include it explicitly in the noise covariance matrix. In actual analyses of real data, the noise estimation is a very important step and we do not attempt to address it in this work.

Fig. 10 shows the distribution of the detection statistics |$\mathcal {P}_{\rm {mon}}$| in the presence of signal along with the distribution of corresponding squared null streams (⁠|$|{\tilde{\bf d}}_{\rm {null}}|^2$|⁠) as defined in equation (15). The injected signal is characteristic of a circular SMBBH described by the following parameters: h₀ = 1 × 10⁻¹⁴, f = 20 nHz, cos ι = 0.6172, ψ = 4.1991, ϕ₀ = 1.9756, (α, δ) = (0, 0). The expected signal-to-noise ratio is ρ = 10, which is identical to the example shown in Section 4.1. To obtain the distribution of the detection statistic and null streams, we perform 10⁴ realizations of white Gaussian noise and for each noise realization keep the search parameters (i.e. source frequency and sky direction) fixed at their injected values. We can see that both match the expected distributions perfectly well. As mentioned in the previous section, null streams can be used as a consistency check in the case of a detected GW signal – if the detected signal is due to a GW, null streams should follow the expected noise distribution while otherwise false alarms caused by any noise processes are very unlikely to exhibit such a property.

$(a) Probability distribution of the detection statistic $\mathcal {P}_{\rm {mon}}$ (solid; solid blue in online) in the presence of a monochromatic signal compared to the expected non-central χ2 distribution (dash; dashed red in online) with 4 degrees of freedom and an non-centrality parameter ρ2 = 100. The simulation was performed with search parameters fixed at the injected values and 104 realizations of white Gaussian noise. (b) For the same simulation as panel (a) but instead showing the probability distribution of the squared null streams $|{\tilde{\bf d}}_{\rm {null}}|^2$ (solid; solid blue in online) compared to the expected χ2 distribution with 2 degrees of freedom (dash; dashed red in online).$

Figure 10.

(a) Probability distribution of the detection statistic |$\mathcal {P}_{\rm {mon}}$| (solid; solid blue in online) in the presence of a monochromatic signal compared to the expected non-central χ² distribution (dash; dashed red in online) with 4 degrees of freedom and an non-centrality parameter ρ² = 100. The simulation was performed with search parameters fixed at the injected values and 10⁴ realizations of white Gaussian noise. (b) For the same simulation as panel (a) but instead showing the probability distribution of the squared null streams |$|{\tilde{\bf d}}_{\rm {null}}|^2$| (solid; solid blue in online) compared to the expected χ² distribution with 2 degrees of freedom (dash; dashed red in online).

Open in new tab Download slide

Fig. 11 shows the detection statistics as a function of frequency calculated for realistic simulated data sets in the absence or presence of the same monochromatic signal as in Fig. 10. We have considered two statistics – |$\mathcal {P}_{\rm {mon}}$| proposed in this paper and the |$\mathcal {F}_{e}$|-statistic which was first proposed in Babak & Sesana (2012) and later strengthened in Ellis et al. (2012), both giving nearly identical results. In the signal-present case, the maximum value (129.4) of |$\mathcal {P}_{\rm {mon}}$| is found at 20.5 nHz, while the largest |$\mathcal {F}_{e}$|-statistic (125.4) is measured at 20.6 nHz. It is also obvious that at some high frequencies there is a spectral leakage problem that applies to both methods. This is not surprising given the similar principle of the two statistics: (1) the |$\mathcal {F}_{e}$|-statistic works fully in the time domain, but it also involves a process of maximum-likelihood estimation of the Fourier components in individual pulsar data set; (2) in the derivation of the |$\mathcal {F}_{e}$|-statistic, the four (time-dependent) basis functions of the monochromatic timing residuals are essentially equivalent to the two orthogonal sine-cosine pairs of A_{+, ×}(t).

$Detection statistics as a function of frequency for realistic simulated data sets in the absence (thin curves) and presence (thick curves) of a strong monochromatic signal. Here, we compare two statistics – $\mathcal {P}_{\rm {mon}}$ (solid; solid red in online) proposed in this paper and the $\mathcal {F}_{e}$-statistic (dash; black dash in online). All statistics are evaluated at the injected source sky direction. The vertical line marks the injected frequency (20 nHz), while the horizontal line corresponds to a single-trial FAP of 10−4.$

Figure 11.

Detection statistics as a function of frequency for realistic simulated data sets in the absence (thin curves) and presence (thick curves) of a strong monochromatic signal. Here, we compare two statistics – |$\mathcal {P}_{\rm {mon}}$| (solid; solid red in online) proposed in this paper and the |$\mathcal {F}_{e}$|-statistic (dash; black dash in online). All statistics are evaluated at the injected source sky direction. The vertical line marks the injected frequency (20 nHz), while the horizontal line corresponds to a single-trial FAP of 10⁻⁴.

Open in new tab Download slide

Fig. 12 shows the results of an all-sky search using the SVD method for the same signal-present data set as used in Fig. 11. The maximum detection statistic 130.1 is found at (α, δ) = (6.185, − 0.133) with a frequency of 20.5 nHz. The source is localized to a direction close to the injected location. However, compared with Fig. 2, we can see that the angular resolution for realistic PTAs is much worse because the sensitivity is dominated by a few good ‘timers’ in the array. It is worth pointing out that the SVD method is also advantageous in terms of computational efficiency over that of the |$\mathcal {F}_{e}$|-statistic, since for the latter the amount of computation is proportional to the total number of sky locations being searched.

Figure 12.

Sky map of detection statistics calculated for a simulated realistic data set that includes a strong monochromatic signal. The source is simulated at the centre of the map (indicated by a white diamond) and the maximum detection statistic is found at ‘○’. Sky locations of the 20 PPTA pulsars are marked with ‘★’.

Open in new tab Download slide

Similar results were obtained if the simulated data sets used in Figs 10–12 have gone through the tempo2 fitting process for a full timing model for each pulsar. In fact, it has been shown that the fitting process can be approximated by multiplying the timing residuals by a data-independent and non-invertible linear operator matrix (e.g. Demorest 2007). This matrix can be explicitly included in the noise covariance matrix and does not pose a problem for our analysis method. The timing model fit is known to absorb some GW power at low frequencies (close to 1/T_span with T_span being the data span, because of the fit for pulsar spin period and its first time derivative), and in some narrow bands centred around 1 yr^{− 1} (the fit for pulsar positions and proper motions), 2 yr^{− 1} (the fit for pulsar parallax) and high frequencies that correspond to binary periods for pulsars in binaries (e.g. fig. 1 in Cutler et al. 2014). These effects are similar to those caused by applying constraints on the |$A_{+,\times }^{{\rm {t2}}}(t)$| time series as discussed in Section 3.1.

6 SENSITIVITIES TO ECCENTRIC SMBBHs

In this section, we are interested in current PTA sensitivities to eccentric binaries. Given that previous published PTA searches for GWs from SMBBHs have assumed zero eccentricity, we attempt to answer the following question – when is a monochromatic search sufficient to detect eccentric binaries? For this purpose, we use simulated data sets that have the same sampling and error bars as the PPTA DR1 data set, and consider two detection statistics |$\mathcal {P}_{\rm {mon}}$| and |$\mathcal {P}_{\rm {ecc}}$|⁠. The sensitivities are parametrized by the luminosity distance within which 95 per cent of sources are detectable at a fixed FAP of 10⁻⁴.

To facilitate our calculations, we simulate signals due to eccentric binaries drawn from uniform distribution in cos δ and α and evaluate both |$\mathcal {P}_{\rm {mon}}$| and |$\mathcal {P}_{\rm {ecc}}$| for noiseless data at the injected sky location. For the demonstration here, we consider sources with a chirp mass of 10⁹M_⊙ and an orbital frequency of 10 nHz with other parameters such as cos ι, initial orbital phase and polarization angle randomized. Since here the data consist of signals only, the detection statistic equals the squared signal-to-noise ratio ρ² and scales inversely with |$d_{\rm L}^2$|⁠. We use this scaling relation to find the value of d_L that corresponds to the given detection threshold. We perform 10⁴ Monte Carlo simulations and choose the 95 per cent quantile as a point in the sensitivity plot.

Fig. 13 shows the sensitivity to GWs from eccentric SMBBHs using two detection strategies – a monochromatic search and a harmonic summing technique. It is clearly demonstrated that for high and low eccentricities a monochromatic search is better than the harmonic-summing search, while the latter is more sensitive in the moderate regime (0.3 ≲ e₀ ≲ 0.55). This is in good agreement with the fact that for high and low eccentricities the GW emission is dominated by the orbital frequency and its second harmonic, respectively. In the high- and low-eccentricity regimes, adding incoherently power from secondary harmonics is not beneficial because when a harmonic summing technique is adopted the detection threshold should be increased for a fixed FAP. We also find that the inclusion up to the third harmonic is sufficient for all possible eccentricities in this example.

$Luminosity distance (dL), within which 95 per cent of eccentric SMBBHs with a chirp mass of 109M⊙ and an orbital frequency of 10 nHz could be detectable with the current PTA sensitivity, as a function of orbital eccentricity (e0). Two detection strategies are considered: a monochromatic search (black dots) and a harmonic summation technique (blue open circles). Note that the sensitive distance scales as $M_{c}^{5/3}$.$

Figure 13.

Luminosity distance (d_L), within which 95 per cent of eccentric SMBBHs with a chirp mass of 10⁹M_⊙ and an orbital frequency of 10 nHz could be detectable with the current PTA sensitivity, as a function of orbital eccentricity (e₀). Two detection strategies are considered: a monochromatic search (black dots) and a harmonic summation technique (blue open circles). Note that the sensitive distance scales as |$M_{c}^{5/3}$|⁠.

Open in new tab Download slide

7 SUMMARY AND FUTURE WORK

PTA experiments have steadily improved their sensitivities over the past few years and produced very stringent constraints on the evolution of ensembles of SMBBHs. While the traditional focus was a stochastic background from SMBBHs throughout the Universe, possibilities of detecting and studying single-source GWs using PTAs have been explored in depth in recent years. Although it may be more challenging, detections of individual sources will provide rich information on the properties of the GW emitters such as the orbital eccentricities, masses and even spins of SMBBHs. This prospect becomes increasingly important as the Chinese 500-metre aperture spherical telescope (FAST) is expected to come online in 2016 (Nan et al. 2011; Hobbs et al. 2014) and the planned Square Kilometre Array (SKA) will be operational within a decade (Lazio 2013). Both FAST and SKA will provide a major step forward in terms of single-source GW detection with PTAs.

Based on the coherent network analysis method used in ground-based interferometers, we proposed a method that is capable of doing detection, localization and waveform estimation for various types of single-source GWs. We demonstrated the effectiveness of this method with proof-of-principle examples for GWs from SMBBHs in circular or eccentric orbits and GW bursts. We also demonstrated the implementation of this technique using realistic data sets that include effects such as uneven sampling and varying data spans and error bars. The new method is found to have the following features: (1) it is fast to run especially for all-sky blind searches; (2) it performs as well as published time-domain methods for realistic data sets; and (3) null streams can be constructed as a consistency check in the case of detected GW signals. Finally, we presented sensitivities to eccentric SMBBHs and found that (1) a monochromatic search that is designed for circular binaries can efficiently detect SMBBHs with both high and low eccentricities; and (2) a harmonic summing technique provides better detection sensitivities for moderate eccentricities.

Our future work will be:

to develop a fully functional data analysis pipeline that can be tested in future IPTA mock data challenges⁵ and applied to real data;
to make comprehensive comparisons (in terms of detection sensitivity, parameter estimation and speed) between the method proposed in this paper and published methods such as the ‘A₊A_×’ method and the |$\mathcal {F}_{e}$|-statistic. Again, this could be done along with a future IPTA data challenge;
to include pulsar terms in the detection framework for continuous GWs and investigate how the information of pulsar terms can be exploited to improve the PTA angular resolution.

We thank the referee for useful comments. This work was supported by iVEC through the use of advanced computing resources located at iVEC@UWA. We are grateful to Bill Coles and Yan Wang for useful discussions. X-JZ acknowledges the support of an University Postgraduate Award at UWA. LW acknowledges funding support from the Australian Research Council. GH is supported by an Australian Research Council Future Fellowship. J-BW is supported by West Light Foundation of CAS XBBS201322 and NSFC project no.11403086.

Here, it is assumed that noise for different pulsars is uncorrelated. If not the full covariance matrix should be used in a whitening process.

Here, the indices 1 and 2 refer to the first two terms in the new data streams |${\tilde{\bf d}}$|⁠. More indices are used for time and frequency when necessary.

Simultaneous to the estimation and correction of the dispersion measure variations for multiple-band pulsar timing data (Keith et al. 2013), one obtains estimation of the common-mode signal that is independent of radio wavelengths. These common-mode data sets are evenly sampled and can be searched for GWs.

http://dx.doi.org/10.4225/08/534CC21379C12

http://www.3ipta4gw.org/

REFERENCES

Aasi

, et al. ,

2013

preprint (arXiv:1304.0670)

Abadie

, et al. ,

Class. Quantum Gravity

2010

, vol.

pg.

173001

Month:	Total Views:
February 2017	2
March 2017	3
April 2017	3
May 2017	11
July 2017	3
August 2017	5
September 2017	11
October 2017	5
November 2017	5
December 2017	5
January 2018	2
February 2018	5
March 2018	5
April 2018	2
May 2018	9
June 2018	8
July 2018	20
August 2018	15
September 2018	10
October 2018	9
November 2018	11
December 2018	11
January 2019	7
February 2019	14
March 2019	18
April 2019	20
May 2019	17
June 2019	17
July 2019	13
August 2019	13
September 2019	13
October 2019	24
November 2019	20
December 2019	13
January 2020	11
February 2020	14
March 2020	21
April 2020	10
May 2020	3
June 2020	14

Article Contents

Abstract

1 INTRODUCTION

2 SINGLE-SOURCE GWs AND PULSAR TIMING

3 THE DATA ANALYSIS METHOD

3.1 Relation to the ‘A+A×’ method

4 EXAMPLES USING IDEALIZED DATA SETS

4.1 Monochromatic waves

4.2 Eccentric binaries

4.3 GW bursts

5 MORE REALISTIC DATA SETS

6 SENSITIVITIES TO ECCENTRIC SMBBHs

7 SUMMARY AND FUTURE WORK

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

3.1 Relation to the ‘A₊A_×’ method