[
Abstract
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 temporalcorrelated signals; instead, they appear as white noise in the timing residuals. The variance of the GWinduced 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 (superNyquist 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 allsky 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 superNyquistfrequency gravitational waves using a pulsar timing array]Detecting superNyquistfrequency gravitational waves using a pulsar timing array 1,2]Yi ShuXu 1*,3]Zhang ShuangNan
[
[ 
\hb@xt@1ex^{1}Key Laboratory of Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences, Beijing 100049, China; 
\hb@xt@1ex^{2}University of Chinese Academy of Sciences, Beijing 100049, China; 
\hb@xt@1ex^{3}Space 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
Citation:
S.X Yi & S.N. Zhang \@titlehead. Sci ChinaPhys 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 GWimprinted 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 temporalcorrelated 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 superNyquistfrequency GWs (SNFGWs). However, an SNFGW will indicate itself by increasing the total white noise level in the timing residuals. The amplitude of additional GWinducedwhite 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) [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 SNFGWinduced timing noise
The total white timing noise consists of three parts: 1. white noise resulting from the intrinsic properties of pulses, e.g., singlepulse 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
(1) 
where are the geometric factors [18]
(2)  
where is the angle between the GW source and the pulsar, and are
(3)  
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
(4) 
We further denote
(5) 
and
(6) 
then Equation (4) becomes
(7) 
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
(8) 
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 abovementioned 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 J04374715 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 superNyquist. The timing residual at each observation is assigned as
(9) 
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 abovementioned 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 instrumentrelated 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 J18242452A 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
(10) 
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) 

J04374715  69  42  41 
J06130200  1301  19  1042 
J07116830  4405  3  3357 
J1022+1001  2315  625  1327 
J10240719  2981  7  2279 
J10454509  3230  15  2596 
J16003053  758  9  540 
J16037202  2207  24  1283 
J16431224  2722  3  2022 
J1713+0747  424  9  269 
J17302304  2296  4  1677 
J17325049  3224  2  2585 
J17441134  920  6  573 
J18242452A  2337  13  1687 
J1857+0943  1384  8  1292 
J19093744  255  8  232 
J1939+2134  402  295  142 
J21243358  3641  8  2602 
J21295721  3703  5  3017 
J21450750  3532  3  2175 
name  (ns)  (ns)  Ave. TOA (ns) 

J1857+0943  3492  11  2050 
J06130200  2884  3  2058 
J16003053  1549  87  1303 
J1713+0747  1842  143  817 
J19093744  1654  317  811 
J21450750  7173  8  4836 
J1955+2908  6104  5001  7107 
J1012+5307  7493  11  4457 
J1640+2224  6353  10  2753 
J17441134  4530  2  2341 
J1910+1256  3308  17  2233 
J2317+1439  1031  31  629 
J0030+0451  2646  3  2080 
J14553330  12747  5001  13868 
J16431224  3571  2039  2243 
J1853+1308  4208  2  3431 
J19180642  7172  9  5180 
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) 

J04374715  42  J1857+0943  8 
J06130200  19  J19093744  8 
J07116830  3  J21243358  8 
J10240719  7  J21295721  5 
J10454509  15  J21450750  3 
J16003053  9  J1012+5307  11 
J16037202  24  J1640+2224  10 
J16431224  3  J1910+1256  17 
J1713+0747  9  J2317+1439  31 
J17302304  4  J0030+0451  3 
J17325049  2  J1853+1308  2 
J17441134  6  J19180642  9 
J18242452A  13 
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 equalarea 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
(11) 
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 allsky 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 loglogWCC 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 allsky map of loglogWCC 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 sourcedetecting 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 allsky 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 superNyquist 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
5.1 Summary

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 allsky map of the WCC between and is shown in Figure 7.

The allsky 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:

SNFGWs leave additional white noise in the timing residuals.

of the GWinduced white noise is proportional to the parameter , and the  proportional relationship is scaled by . For definitions of the other symbols see Section 1.
Simulation:

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

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 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 GWinduced 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]
(12) 
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]
(13) 
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 GWinduced timing residuals ). From Equations (12) and (13) we know that
(14) 
Therefore,
(15) 
Using Equation (13) we get the upper frequency limit
(16) 
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 highfrequency end of the detectable GW range by traditional pulsar timing methods. The impact of relaxing the stationaryamplitude 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
(17) 
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 highpass 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
(18) 
The total rms contributed by both the pulsar term and the Earth term can also be written as
(19) 
However, in Equation (19) is defined differently compared with that in Equation (8):
(20) 
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 highpass 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 superNyquist 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 startup grant 292012312D1117210, and the Strategic Priority Research Program “The Emergence of Cosmological Structures” (Grant No. XDB09000000) of the Chinese Academy of Sciences.
References
 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 ChinaPhys 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 ChinaPhys Mech Astron, 2015, 58: 120402
 8 Mitrofanov V P, Chao S, Pan HW, et al. Technology for the next gravitational wave detectors. Sci ChinaPhys Mech Astron, 2015, 58: 120404
 9 Blair D, Ju L, Zhao C N, et al. The next detectors for gravitational wave astronomy. Sci ChinaPhys 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., BurkeSpolaor, 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