[    [

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.

.Article .

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) [1] 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 [16], 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.

Figure 1: Upper panel: When , the timing residuals (blue line) show a temporal-correlated structure following the wave form of the GW (red line). Middle panel: When , the timing residuals are white noise. Bottom panel: The weaker strain of SNFGWs leaves lower level of white noise in the timing residuals, compared to that in the middle panel.

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) [17] and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav dfg+12) publicly available datasets [12]. 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 [19], 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 [18]


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 .

Figure 2: Mean values of corresponding to 1,000 randomly distributed pulsars and random as a function of . 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.

Figure 3: Relationship between the timing residuals variance and of simulated timing residuals of pulsars with different injected GW strain . The values of are indicated in each panel.

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.

Figure 4: Fitted slopes of the -variance relationship, i.e., the estimated as functions of the injected . Different intrinsic white noise levels are indicated in each panel.

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 [20] 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)
J0437-4715 69 42 41
J0613-0200 1301 19 1042
J0711-6830 4405 3 3357
J1022+1001 2315 625 1327
J1024-0719 2981 7 2279
J1045-4509 3230 15 2596
J1600-3053 758 9 540
J1603-7202 2207 24 1283
J1643-1224 2722 3 2022
J1713+0747 424 9 269
J1730-2304 2296 4 1677
J1732-5049 3224 2 2585
J1744-1134 920 6 573
J1824-2452A 2337 13 1687
J1857+0943 1384 8 1292
J1909-3744 255 8 232
J1939+2134 402 295 142
J2124-3358 3641 8 2602
J2129-5721 3703 5 3017
J2145-0750 3532 3 2175
Table 1: , and the average of TOA uncertainties (Ave. TOA) of PPTA DR1; the unit is nanosecond.
Figure 5: Timing residuals of PPTA DR1 pulsars. The names of the pulsars are indicated at the top of each panel. The axes are hidden for clarity.
Figure 6: Timing residuals of NANOGrav dfg+12 pulsars. The names of the pulsars are indicated at the top of each panel. The axes are hidden for clarity.
name (ns) (ns) Ave. TOA (ns)
J1857+0943 3492 11 2050
J0613-0200 2884 3 2058
J1600-3053 1549 87 1303
J1713+0747 1842 143 817
J1909-3744 1654 317 811
J2145-0750 7173 8 4836
J1955+2908 6104 5001 7107
J1012+5307 7493 11 4457
J1640+2224 6353 10 2753
J1744-1134 4530 2 2341
J1910+1256 3308 17 2233
J2317+1439 1031 31 629
J0030+0451 2646 3 2080
J1455-3330 12747 5001 13868
J1643-1224 3571 2039 2243
J1853+1308 4208 2 3431
J1918-0642 7172 9 5180
Table 2: , and the average of TOA uncertainties (Ave. TOA) of NANOGrav dfg+12; the unit is nanosecond.

From Tables 1 and 2 we chose pulsars with ns, and, for pulsars shared by both tables, the smaller values are selected. The resulting pulsars are listed in Table 3.

name (ns) name (ns)
J0437-4715 42 J1857+0943 8
J0613-0200 19 J1909-3744 8
J0711-6830 3 J2124-3358 8
J1024-0719 7 J2129-5721 5
J1045-4509 15 J2145-0750 3
J1600-3053 9 J1012+5307 11
J1603-7202 24 J1640+2224 10
J1643-1224 3 J1910+1256 17
J1713+0747 9 J2317+1439 31
J1730-2304 4 J0030+0451 3
J1732-5049 2 J1853+1308 2
J1744-1134 6 J1918-0642 9
J1824-2452A 13
Table 3: Selected 25 pulsars in this work.

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.

Figure 7: All-sky map of the WCC between and . The green stars indicate the locations of the pulsars used; the green circle is the position where the WCC reaches its maximum.

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:

  1. 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.

  2. Starting from a small GW strain value and a random polarization angle , generate a series of timing residuals based on Equation (9).

  3. Follow the SNFGW source-detecting procedure described in the above section. Increase and return to step 2, until the detection significance reaches 99%.

  4. 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.

Figure 8: Normalized probability density distribution (PDF) of the maximum WCC between - of 1,000 permutations of . The vertical dashed line indicates the observed WCC 0.31; the probability that any permutation has the maximum WCC larger than the observed one is 65.5%. Therefore, the probability that the observed - correlation is the consequence of intrinsic white noise of the pulsars is 65.5%, and this result indicates a nondetection.
Figure 9: All-sky map of the sensitivity of our data to single SNFGW sources. The green stars indicate the locations of the pulsars used; the green circle is the position where our PTA is most sensitive to the GW source.

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).

Figure 10: The sensitivity of the method described in this work to single gravitational wave sources for the optimal binary orientation, compared with the sensitivity or upper limits of other methods. The red star marks the peak signal of the first GW event GW150914 [1]. The figure is edited from [21].

5 Summary, Conclusions, and Discussion

5.1 Summary

  1. The theory of the method for detecting SNFGWs is presented in Section 1.

  2. 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.

  3. The all-sky map of the WCC between and is shown in Figure 7.

  4. The all-sky map of the sensitivity of the selected PTA to single SNFGW sources is shown in Figure 9.

5.2 Conclusions

We summarize our conclusions as follows: Theory:

  1. SNFGWs leave additional white noise in the timing residuals.

  2. 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.


  1. When is given, the stronger GW strain will give a more significant - correlation.

  2. The combination of the GW source can be estimated by fitting the - relationship.

  3. 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.


  1. The coordinates of the GW source where the WCC is optimized are , (rad), and the corresponding WCC is 0.31.

  2. 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:

  1. 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 Discussion

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 [18]


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 [22]


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


However, in Equation (19) is defined differently compared with that in Equation (8):


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.


  • Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters, 116, 061102
  • Weber, J. 1977, Nature, 266, 243
  • Abbott, B. P., Abbott, R., Adhikari, R., et al. 2009, Reports on Progress in Physics, 72, 076901
  • Shaddock, D. A. 2009, PASA, 26, 128
  • 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
  • Manchester, R. N., & IPTA 2013, Classical and Quantum Gravity, 30, 224010
  • Blair D, Ju L, Zhao C N, et al. Gravitational wave astronomy: the current status. Sci China-Phys Mech Astron, 2015, 58: 120402
  • 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
  • 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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description