The maximum frequency of gravitational waves (GWs) detectable with traditional pulsar timing methods is set by the Nyquist frequency () of the observation. Beyond this frequency, GWs leave no temporal-correlated signals; instead, they appear as white noise in the timing residuals. The variance of the GW-induced white noise is a function of the position of the pulsars relative to the GW source. By observing this unique functional form in the timing data, we propose that we can detect GWs of frequency (super-Nyquist frequency GWs; SNFGWs). We demonstrate the feasibility of the proposed method with simulated timing data. Using a selected dataset from the Parkes Pulsar Timing Array data release 1 and the North American Nanohertz Observatory for Gravitational Waves publicly available datasets, we try to detect the signals from single SNFGW sources. The result is consistent with no GW detection with 65.5% probability. An all-sky map of the sensitivity of the selected pulsar timing array to single SNFGW sources is generated, and the position of the GW source where the selected pulsar timing array is most sensitive to is , (rad); the corresponding minimum GW strain is at Hz.
Detecting super-Nyquist-frequency gravitational waves using a pulsar timing array]Detecting super-Nyquist-frequency gravitational waves using a pulsar timing array 1,2]Yi Shu-Xu 1*,3]Zhang Shuang-Nan
|\hb@xt@1ex1Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China;|
|\hb@xt@1ex2University of Chinese Academy of Sciences, Beijing 100049, China;|
|\hb@xt@1ex3Space Science Division, National Astronomical Observatories of China, Chinese Academy of Sciences, Beijing 100012, China;|
Received ; accepted ; published
gravitational wave, pulsar, black hole
PACS number(s): 25.60.Je, 21.10.Jx, 25.40.Lw, 26.20.+f
S.-X Yi & S.-N. Zhang \@titlehead. Sci China-Phys Mech Astron, 2016, 58: ??, doi: ??
Introduction The recent direct detection of gravitational waves (GWs)  marks the beginning of GW astronomy era, after about five decades of effort on GW detection [2, 3, 4, 6, 7, 5, 8, 9]. Among the various proposed methods, the pulsar timing array (PTA) method shows promise in identifying GW-imprinted structure in the timing residuals of a number of pulsars [10, 11]. Although no detection has been made, more and more stringent upper limits of both individual GWs and GW background have been set using this method [12, 13, 14, 15]. The traditional pulsar timing method has an upper limit on the frequency of detectable GWs, which is known as the Nyquist frequency (). In the case of even sampling, , where is the observation time span, and is the number of the pulse time of arrival (TOA). In general, is set by the sampling rate of the TOA. When , where is the frequency of GWs, GWs leave no temporal-correlated structures in the timing residuals of the pulsar and are therefore undetectable using the traditional pulsar timing method; instead, additional white noise will be left in the timing residuals, as the TOAs are not coherent in phase with the GWs (see the illustration in Figure 1). For the most typical biweekly observation scheme, the frequency upper bound is Hz. Although, in principle, GWs can be searched for at arbitrarily high frequency in the TOA using a Bayesian method , the parameters of white timing noise will be completely correlated with the amplitude of the GWs, once , i.e., for super-Nyquist-frequency GWs (SNFGWs). However, an SNFGW will indicate itself by increasing the total white noise level in the timing residuals. The amplitude of additional GW-induced-white noise is a function of the coordinates of the pulsars relative to the GW source. In this paper, we propose that, by observing this unique relation between pulsar position and white noise variance, we can detect an SNFGW using a PTA.
The main purpose of this paper is to present the theory of the proposed method and also to demonstrate its feasibility with available data. The paper is organized as follows: The theory of our method is described at Section 2. We test the feasibility of our method on simulated timing data in Section 3. A subset of pulsars is selected from the Parkes Pulsar Timing Array data release 1 (PPTA DR1)  and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav dfg+12) publicly available datasets . The selected PTA is used to search for SNFGWs using the proposed method in Section 4. In Section 5, we study the sensitivity of our method and the data to single SNFGW sources. We summarize our conclusions and discuss potential targets and the advantages and shortcomings of the proposed method in the last section.
1 Relation between the coordinates of pulsars and the SNFGW-induced timing noise
The total white timing noise consists of three parts: 1. white noise resulting from the intrinsic properties of pulses, e.g., single-pulse shape variability , or from propagation through the interstellar medium (ISM) (); 2. measurement uncertainties of the TOAs (); and white noise induced by SNFGWs (). Since each is an intrinsic property of a pulsar, we expect that there is no correlation between from pulsar to pulsar. Now we consider how the SNFGW determines of a pulsar. Suppose that the right ascension and the declination of the pulsar are and , respectively, and that those of the GW source are and , respectively. The timing residuals brought by the GW are
where are the geometric factors 
where is the angle between the GW source and the pulsar, and are
In Equation (3), and are, respectively, the strain and the angular frequency of the GW, is the polarization angle of the GW, and is the inclination angle of the orbital plane of the GW source with respect to the line of sight. The initial phase of the GW is set to zero in Equation (3). We defer the treatment of the pulsar term to the discussion section; for now, we only consider the Earth term in Equation (3) for simplicity. Equation (3) can be simplified as follows: We denote , , and ; therefore, Equation (1) can be rewritten as
We further denote
then Equation (4) becomes
Since the sampling frequency of TOAs is less than the frequency of the sinusoid in Equation (7), and the TOAs are random in the phase of the sinusoid, the resulting extra variance in the timing residuals is
Among the three components of white noise, and are not correlated with the position of the pulsars, i.e., the parameter ; meanwhile, is proportional to . Therefore, as long as the strain of the GW is strong enough to make dominate, , where is the rms of the timing residuals. Note that is in fact also a function of the positions of the pulsars; therefore, the correlation between and will deviate from proportionality. as a function of is plotted in Figure 2: 1,000 pulsars were generated randomly in a uniform distribution in the sky, and is uniformly random from 0 to . The position of the GW source is assigned to . The corresponding is calculated using Equation (5). In Figure 2, the vertical axis indicates the mean values of , and the error bars are the standard deviation of .
When the inclination angle is far from 90, can be treated as a function of , whereas when , the scatter of around its mean value becomes more and more significant, which is the systematic scatter of Equation (8).
2 Testing the method with simulated data
To test the feasibility of the above-mentioned method, we generate simulated pulsar timing residuals for 1,000 pulsars. The pulsars are uniformly distributed at random over the celestial sphere. The time spans, the total number of observations, the intrinsic white noise levels, and the TOA uncertainties are assigned according to the data of PSR J0437-4715 in PPTA DR1, while the dates of observations are randomly assigned for each pulsar. The coordinates of the GW source are set at , . Since the pulsars are uniformly distributed, the location of the GW source does not matter. The inclination angle of the source is set to , and the polarization angle is , Hz, and Hz; therefore, is super-Nyquist. The timing residual at each observation is assigned as
where is a Gaussian with a mean and a standard deviation ; is the intrinsic white noise level, which is set to ns; is the TOA uncertainty at each . Figure 3 shows the relationship between the variance of timing residuals and the position parameter of the simulated timing data of the pulsars. The correlation becomes more and more significant with increasing GW strain.
The slope of the above-mentioned -variance relationship is equal to the combination of the GW, where the factor is included to take account of for different pulsars. We inject different GW strain into the timing residuals and perform linear fit to the resulting -variance pairs, and we plot the fitted slope, i.e., the estimated as a function of the injected value of in Figure 4, where is the parameter when is fixed to in Equation (4). The error bars show the error of the slope given by the fitting process. When the injected GW strain is small, the fitted slopes have large uncertainties and deviate from the injected values. As the injected GW strain increases enough, the slopes fall onto the red dashed lines, where the estimated equals the injected . By changing the intrinsic white noise level , we find that higher intrinsic white noise decreases the sensitivity of the PTA to the GW; this conclusion is in accordance with traditional pulsar timing methods.
3 Detecting SNFGWs using the selected PTA from PPTA DR1 and NANOG1mmr av dfg+12
In this section, we demonstrate how the proposed method can be used with real pulsar timing data by replacing the simulated pulsar timing data with real data selected from PPTA DR1 and NANOGrav dfg+12. As mentioned above, , and all contribute to the total residuals budget. The contribution of can be calculated from the uncertainty of the TOA and can be removed from the total residuals rms. We denote the rms of the remaining residuals as . is composed of and . In the ideal case, values of these pulsars are intrinsic to the individual pulsars and they do not correlate with the PTA. Therefore the only component that makes correlate with the PTA is . However, in practice, owing to the complexity of the observing systems (i.e., front end/back end combinations), the apparent TOA uncertainty cannot faithfully reflect . As a result, the rms of some instrument-related timing residuals enters , which are also correlated in complicated ways but beyond the scope of this current work. Therefore, the estimated significance of GW detection in this section should be considered only as an upper limit or optimistic evaluation. The timing residuals are obtained by fitting the TOA using TEMPO2  and the variance of the timing residuals are thus calculated. The timing residuals of PSR J1939+2134 and J1824-2452A are polynomial whitened, while the data of other pulsars are fitted using the downloaded ephemeris. The bias and scaling factors EQUAD and EFAC are all set to zero in the fitting. The resulting timing residuals of pulsars of PPTA DR1 and NANOGrav dfg+12 are plotted in Figures 5 and 6, respectively. of each pulsar is estimated as follows: We generate a new series of timing residuals for each pulsar such that
from a small starting value; we increase until the variance of the simulated timing residuals (Var) equals the real variance (Var). In practice, we consider these two quantities to be identical when the relative difference (. We list the total rms (), , and the average of the TOA uncertainties in Tables 1 and 2.
|name||(ns)||(ns)||Ave. TOA (ns)|
|name||(ns)||(ns)||Ave. TOA (ns)|
After is known for the th pulsar (where ranges from 1 to , where is the number of pulsars in the PTA), we need to obtain using Equation (2) to test the correlation described in Equation (8). Since the location of the GW source is unknown, we divide the celestial sphere into equal-area grids. For each grid we suppose that the GW source is locates within it and we calculate . We want to test the correlations between and . Owing to the nature of Pearson correlation coefficient (PCC), data with greater distance to the barycenter contribute more to the PCC. As a result, if we use the PCC to study the correlation, the minority of the pulsars that have the largest will dominate the PCC of -. To avoid this problem, we study the correlation on a logarithm scale, in which the scatter of data decreases. In this way, we also assign less weight to the outliers. We therefore use the weighted correlation coefficient (WCC) between and . The WCC between two lists of data and is defined as
where are the normalized weights. The weight that we assign to each data point is its distance to the barycenter. A loop of all the sky grids gives the all-sky map of --WCC, which is plotted in Figure 7; the color scale indicates the WCC. The celestial grid where --WCC is maximum is indicated by the green circle. The coordinates of this point are , (rad), and the corresponding log-log-WCC is 0.31. We need to know the probability that the above - correlation is due to the intrinsic white noise of the PTA. Therefore, we randomly shuffle the values of the pulsars 1,000 times, and we calculate the all-sky map of log-log-WCC and the largest WCC for each permutation. We then get the distribution of the maximum WCC, which is plotted in Figure 8. We notice that 65.5% of the realizations have a maximum WCC larger than the observed value 0.31; thus the probability that the observed - correlation is the consequence of intrinsic white noise of the pulsars is 65.5%, and the result is consistent with a nondetection.
4 Sensitivity to single SNFGW sources
In the above section we found a nondetection result and we want to study the sensitivity of this method to single SNFGW sources using PPTA DR1 and NANOGrav dfg+12 data. The procedure is outlined as follows:
Divide the sky sphere uniformly into grids. In each grid, put an SNFGW source. The inclination angle is set to optimal , and the frequency of the GW is Hz.
Starting from a small GW strain value and a random polarization angle , generate a series of timing residuals based on Equation (9).
Follow the SNFGW source-detecting procedure described in the above section. Increase and return to step 2, until the detection significance reaches 99%.
Record the current value of as the minimum GW strain that the dataset is sensitive to. Move to the next grid point of the sky sphere.
The resulting all-sky map of sensitivity is presented in Figure 9. The position of the GW source where the selected PTA is most sensitive to is , (rad), which is indicated with a green circle in Figure 9; the corresponding minimum is at Hz. According to Equation (3), the sensitive scales with the frequency . We present our sensitivity results in the super-Nyquist band in Figure 10, compared with limits from previous PTA works (the sensitivity curves of the LIGO and the proposed LISA).
5 Summary, Conclusions, and Discussion
The theory of the method for detecting SNFGWs is presented in Section 1.
We tested the theory in Section 2 using simulated PTA data and showed that studying the correlation between and can serve as a method for detecting SNFGWs.
The all-sky map of the WCC between and is shown in Figure 7.
The all-sky map of the sensitivity of the selected PTA to single SNFGW sources is shown in Figure 9.
We summarize our conclusions as follows: Theory:
SNFGWs leave additional white noise in the timing residuals.
of the GW-induced white noise is proportional to the parameter , and the - proportional relationship is scaled by . For definitions of the other symbols see Section 1.
When is given, the stronger GW strain will give a more significant - correlation.
The combination of the GW source can be estimated by fitting the - relationship.
The GW strain needs to be higher than a lower limit to make the - correlation unambiguous. The lower limit of GW strain increases with the intrinsic white noise level of the PTA.
The coordinates of the GW source where the WCC is optimized are , (rad), and the corresponding WCC is 0.31.
Monte Carlo simulation indicates that the probability that the observed - correlation is the consequence of the intrinsic white noise of the pulsars is 65.5%.
Sensitivity to single SNFGW sources:
The position of the GW source where the selected PTA is most sensitive to is , (rad); the corresponding minimum GW strain is at Hz.
5.3.1 Target sources
The SNFGW sources that we aim to study are the merging supermassive black hole binaries (SMBHBs). The frequency of the GW is higher than the typical of a PTA, i.e., Hz. Therefore, we can use filtering techniques to remove red noise from other origins. We treat the strain and as steady throughout the paper (stationary assumption); however, they are both evolving as the SMBHB merges. The demand that the amplitude of GW-induced white noise be stationary during the TOA time span set an upper limit on the GW frequency obtained from this method. We estimate the upper frequency limit as follows. The GW strain is related to its frequency by 
where is a constant scaling factor determined by the chirp mass () and distance of the GW source. The GW frequency at the observer is related to the time before final coalescence of the binary by 
where is the redshift of the GW source (). We want the relative change of in Equation (7) to be less 10% during the time span (so that the change of the variance of the GW-induced timing residuals ). From Equations (12) and (13) we know that
Using Equation (13) we get the upper frequency limit
Inserting , yr, and into Equation (16) gives the upper frequency limit of Hz. Therefore, this method increases the GW frequency upper limit by an order of magnitude without increasing the cadence of observations. Longer decreases ; however, we can divide the whole data span into small segments and apply the method on each segment. If the chirp mass increases to then Hz, which is only a small extension toward the high-frequency end of the detectable GW range by traditional pulsar timing methods. The impact of relaxing the stationary-amplitude condition will be studied in the future.
5.3.2 Pulsar term
When the GW passes the pulsar, a sinusoidal structure similar to that in Equations (1) and (3) will be left in the TOA, which is known as the pulsar term. We denote the frequency of the pulsar term and the Earth term as and , respectively. can be related to by
where is the distance of the pulsar, and is the time to coalescence, and is the angle between the pulsar and the GW source. When , we can use a high-pass filter to remove the contribution of the pulsar term and then process as described above. When , the power contributed by the pulsar term cannot be separated from the Earth term. We denote the amplitude of the signal of the pulsar term and the Earth term as and , respectively, which are related by
The total rms contributed by both the pulsar term and the Earth term can also be written as
When we try to detect the GW source with expected frequency , and , we can calculate its value via Equation (13), and then get . For each pulsar’s , we compare it with . If , we apply a high-pass filter to the timing residuals and thus remove the pulsar term. If , we refine by Equation (20) and process the data as described above. In the second case, the uncertainty of the distance of the pulsar will affect the determination of and thus distort the expected linear correlation in Equation (19). Since we need only the variance of white noise, we use the whitened timing residuals. Therefore a large number of pulsars that are not usable in traditional pulsar timing methods because of their red noise can be included in our treatment, including some normal pulsars that have low noise at super-Nyquist Fourier frequencies. With more pulsars and a wider distribution on the celestial sphere, we have a larger range of and therefore a benefit to GW detection using this method. However, the diversity of the intrinsic white noise level of pulsars also decreases the sensitivity to GW signals.
SNZ acknowledges partial funding support from the National Basic Research Program (“973” Program) of China (Grant Nos. 2014CB845802 and 2012CB821801), the National Natural Science Foundation of China (Grant Nos. 11103019, 11133002, 11103022, and 11373036), the Qianren start-up grant 292012312D1117210, and the Strategic Priority Research Program “The Emergence of Cosmological Structures” (Grant No. XDB09000000) of the Chinese Academy of Sciences.
- 1 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters, 116, 061102
- 2 Weber, J. 1977, Nature, 266, 243
- 3 Abbott, B. P., Abbott, R., Adhikari, R., et al. 2009, Reports on Progress in Physics, 72, 076901
- 4 Shaddock, D. A. 2009, PASA, 26, 128
- 5 Lee H M, Le Bigot E, Du Z H, et al. Gravitational wave astrophysics, data analysis and multimessenger astronomy. Sci China-Phys Mech Astron, 2015, 58: 120403
- 6 Manchester, R. N., & IPTA 2013, Classical and Quantum Gravity, 30, 224010
- 7 Blair D, Ju L, Zhao C N, et al. Gravitational wave astronomy: the current status. Sci China-Phys Mech Astron, 2015, 58: 120402
- 8 Mitrofanov V P, Chao S, Pan H-W, et al. Technology for the next gravitational wave detectors. Sci China-Phys Mech Astron, 2015, 58: 120404
- 9 Blair D, Ju L, Zhao C N, et al. The next detectors for gravitational wave astronomy. Sci China-Phys Mech Astron, 2015, 58: 120405
- 10 Sazhin, M. V. 1978, Soviet Astronomy, 22, 36
- 11 Detweiler, S. 1979, ApJ, 234, 1100
- 12 Demorest, P. B., Ferdman, R. D., Gonzalez, M. E., et al. 2013, ApJ, 762, 94
- 13 Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2014, Apj, 794, 141
- 14 Lentati, L., Taylor, S. R., Mingarelli, C. M. F., et al. 2016, MNRAS, 453, 2576
- 15 Zhu, X.-J., Wen, L., Hobbs, G., et al. 2015, MNRAS, 449, 1650
- 16 van Haasteren, R. 2011, Ph.D. Thesis, 176
- 17 Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017
- 18 Lee, K. J., Wex, N., Kramer, M., et al. 2011, MNRAS, 414, 3251
- 19 Shannon, R. M., Osłowski, S., Dai, S., et al. 2014, MNRAS, 443, 1463
- 20 Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006, MNRAS, 369, 655
- 21 Thrane, E., & Romano, J. D. 2013, Phys. Rev. D., 88, 124032
- 22 Hughes, S. A. 2009, ARA&A, 47, 107