Extending the Frequency Reach of Pulsar Timing Array-based Gravitational Wave Search without High-cadence Observations

The ever-growing trove of gravitational wave (GW) signals from compact binary coalescences (Abbott et al. 2016, 2019) collected by the Laser Interferometer Gravitational-Wave Observatory (LIGO; LIGO Scientific Collaboration et al. 2015) and Virgo (Acernese et al. 2015) detectors is revealing the GW universe in the ∼10 to 10³ Hz band. At lower frequencies, the space-based Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) mission will target the millihertz band from 10⁻⁴ to 10⁻¹ Hz while pulsar timing arrays (PTAs) are already putting meaningful constraints in the sub-μHz band on the stochastic GW background from an unresolved population of supermassive black hole binaries (SMBHBs; Lentati et al. 2015; Shannon et al. 2015; Arzoumanian et al. 2016, 2018), continuous waves from resolvable SMBHBs (Zhu et al. 2014; Babak et al. 2016; Aggarwal et al. 2019), and bursts (Wang et al. 2015a; Aggarwal et al. 2020).

The numbers of millisecond pulsars currently being timed by PTA consortia are 47 (Alam et al. 2020, NANOGrav), 26 (Kerr et al. 2020, PPTA), 42 (Desvignes et al. 2016, EPTA), and 65 (Perera et al. 2019, IPTA). Next-generation radio telescopes, namely the Five-hundred meter Aperture Spherical Telescope (Nan et al. 2011; Hobbs et al. 2019) and the Square Kilometer Array (SKA; Smits et al. 2009; Janssen et al. 2015) will grow the number of well-timed pulsars (noise rms ≲100 ns) to O(10³). Along with a more uniform sky coverage and standardized data spans, this will improve the sensitivity to GWs from resolvable sources by two orders of magnitude (Wang & Mohanty 2017, 2020a).

The high-frequency limit of the sensitive band for PTA-based GW searches is widely assumed (Zhu et al. 2014; Lentati et al. 2015; Shannon et al. 2015; Arzoumanian et al. 2016, 2018; Babak et al. 2016; Aggarwal et al. 2019) to be ≈4 × 10⁻⁷ Hz, corresponding to the Nyquist rate (Bracewell 2000) associated with the average cadence of timing observations, typically 1/(2 weeks), for individual pulsars in an array. Therefore, attempts at extending the high-frequency limit for resolvable GW sources, to frequency >1 μHz, are all based on high-cadence observations of single pulsars (Yardley et al. 2010; Yi et al. 2014; Dolch et al. 2016; Perera et al. 2018). It should be noted here that the actual cadences for pulsars in current PTAs vary over a large range and 1/(2 weeks) is more representative of its higher end.

In this Letter, we show that the high-frequency reach of PTAs for resolvable sources is not limited by the Nyquist rate, f_SP, of single-pulsar observations and that the limiting frequency, f_PTA, can be much higher than assumed so far. The key here is that a higher frequency signal is preserved due to aliasing in the sequence of timing observations from each pulsar and can be unscrambled using asynchronous observations (Bretthorst 2001; Wong et al. 2006) from multiple pulsars.

The lack of synchronicity, an inherent feature of PTA data, can be turned into an observational strategy, which we call staggered sampling, to boost the high-frequency reach of PTAs. Staggered sampling simply requires the introduction, by design, of relative time shifts between the sequences of timing observations without requiring a change in the individual observational cadence for any pulsar. Unlike high-cadence observations, this approach is scalable to future large-scale PTAs with O(10³) pulsars because it does not increase the total telescope time consumed by PTA observations.

While the sensitivity of PTA-based GW searches falls with increase in GW signal frequency, an increase in the number of pulsars enhances it. This motivates a first exploration in this Letter of the astrophysical implications of high-frequency searches with an SKA-era PTA, where staggered sampling could increase the frequency reach to ≈4 × 10⁻⁴ Hz and bridge the gap in coverage of the GW spectrum between PTAs and LISA.

In a PTA with N_p pulsars, the timing residual of the I-th pulsar after subtracting a best-fit model (excluding GWs) of the pulse time of arrival is given by d^I(t) = s^I(t) + n^I(t), where s^I(t) is the GW-induced signal and n^I(t) is noise.

At the high signal frequencies of interest to us, the samples of n^I(t) can be assumed to be drawn from an independent and identically distributed Normal random process with zero mean and constant variance ${\left({\sigma }^{I}\right)}^{2}$ (i.e., white Gaussian noise). The contribution to n^I(t) from errors in fitting the timing model are negligible at higher frequencies except at very specific ones (Kopeikin & Potapov 2004; Cutler et al. 2014) such as 1 yr⁻¹ and harmonics. The latter are ignored in our analysis due to the extremely narrow bands that are affected.

With ν^I(t) and ${\nu }_{0}^{I}(t)$ denoting the pulsar rotation frequencies observed at the solar system Barycenter and at the pulsar, respectively, s^I(t) is given by Estabrook & Wahlquist (1975), Sesana & Vecchio (2010)

$$\begin{eqnarray}&&{s}^{I}(t)={\int }_{0}^{t}{dt}^{\prime} {z}^{I}(t^{\prime} ),\end{eqnarray} \tag{ 1 }$$

where ${z}^{I}(t)\equiv ({\nu }^{I}(t)-{\nu }_{0}^{I}(t))/{\nu }_{0}^{I}(t)$ is the GW-induced Doppler shift. For a plane GW arriving from R.A. (α) and decl. (δ), with polarizations h_+,×(t; θ) parametrized by source parameters θ,

$$\begin{eqnarray}&&{z}^{I}(t)=\displaystyle \sum _{A=+,\times }{F}_{A}^{I}(\alpha ,\delta ){\rm{\Delta }}{h}_{A}(t;\theta ),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{h}_{+,\times }(t;\theta )={h}_{+,\times }(t;\theta )-{h}_{+,\times }(t-{\kappa }^{I};\theta ),\end{eqnarray} \tag{ 3 }$$

where, ${F}_{+,\times }^{I}(\alpha ,\delta )$ are the antenna pattern functions (Lee et al. 2011) and Δh_+,×, for the two-pulse response (Estabrook & Wahlquist 1975), contain the so-called Earth and pulsar terms that arise from the action of the GW on pulses at the time, t, of their reception and at the time, t − κ^I, of their emission, respectively.

For a non-evolving circular binary emitting a monochromatic signal, θ includes the overall amplitude (ζ), GW frequency (f_gw), inclination angle of the binary orbital angular momentum relative to the line of sight (ι), GW polarization angle (ψ), and initial orbital phase (φ₀) (Wang et al. 2014). The time delay κ^I appears as an unknown constant phase offset, called the pulsar phase parameter ϕ_I, for such a source.

To search for resolvable GW sources, we use the likelihood-based detection and parameter estimation method described in Wang & Mohanty (2015b), Wang et al. (2017), and Wang & Mohanty (2020b) that takes both the Earth and pulsar terms into account. The method partitions the estimation of parameters such that the pulsar phases are either maximized (Wang & Mohanty 2015b) or, as chosen here, marginalized (Wang et al. 2017) semi-analytically, allowing the method to scale to an arbitrarily large N_p. The remaining parameters are estimated numerically by maximizing the (marginalized) likelihood using Particle Swarm Optimization (Kennedy & Eberhart 1995; Zhu et al. 2016; Mohanty 2018). The maximum value serves as the detection statistic for deciding between the null (H₀) and alternative (H₁) hypotheses about given data that a signal is absent or present, respectively.

Let ${t}^{I}=\{{t}_{i}^{I}\}$, i = 1, 2,..., N^I, denote the times at which the residual d^I(t) is sampled and let their spacing, ${t}_{i+1}^{I}-{t}_{i}^{I}$, be Δ > 0 on the average (e.g., Δ ≥ 2 weeks for IPTA pulsars Verbiest et al. 2016). Consider the set ${ \mathcal T }={\cup }_{I=1}^{{N}_{{\rm{p}}}}{t}^{I}$ of sample times in ascending order from all the array pulsars and let x_k, $k=1,2,\ldots ,\,{\sum }_{I=1}^{{N}_{{\rm{p}}}}{N}^{I}$, denote an element of this set. We consider two specific schemes for staggered sampling in our analysis. In both of them, we set N^I = N to be the same for all the pulsars. This is mainly for reducing the complexity of our codes and not an essential requirement for staggered sampling.

The most straightforward scheme, called uniform staggered sampling, is one where x_i+1 − x_i is a constant. This implies that the samples of d^I(t) are uniformly spaced and the sequence of samples from one pulsar has a constant time shift relative to those of others: ${t}_{i}^{I}=(i-1){\rm{\Delta }}+{\delta }^{I}$, with δ^I = (I − 1)Δ/N_p.

In the second scheme, called randomized staggered sampling, x_i+1 − x_i is a random variable. This is a more realistic situation given the uncertainties inherent in planning astronomical observations. However, as with current PTAs, any real observation strategy would have a target that it seeks to approximate, which is assumed here to be uniform staggered sampling. Therefore, we adopt a reasonable model for randomized staggered sampling in which t^I_i is replaced by ${t}_{i}^{I}+{c}_{i}^{I}$, where c^I_i is a random variable. In our analysis, we will assume that c^I_i is drawn from a truncated Cauchy probability density function (PDF) (Papoulis 1984) with a location parameter set to zero, scale factor of 1/3 day, and $| {c}_{i}^{I}| \leqslant 7$ days. The heavy tails of this PDF allow large excursions—the 99% inter-percentile range is ≃10 days—from the planned observation times of uniform staggered sampling.

The search for individual GW sources is carried out on staggered sampling data as described earlier—no changes are required to the detection and estimation algorithm as it works entirely in the time domain. That this leads to a higher frequency reach is validated directly in this Letter using simulated data. While a rigorous analytic treatment of an arbitrary staggered sampling scheme is left to future work, the following argument indicates what the maximum detectable signal frequency should be. Take the trivial case of identical GW-induced residuals, s^I(t) = s^J(t), observed with uniform staggered sampling. Because the identical residuals lead to a common signal, pooling all data into a single time series will yield samples of the same signal but with a smaller spacing of Δ/N_p. Hence, the maximum detectable frequency is f_PTA = N_pf_SP. Simply pooling the data does not work for the real case of a heterogeneous, s^I(t) ≠ s^J(t) for I ≠ J, set of GW-induced residuals but one expects the same limit to hold.

Note that a higher cadence observational strategy to achieve the same high-frequency limit as staggered sampling, namely N_pf_SP, would increase the total telescope time occupied in timing observations by a factor N_p as each pulsar must be timed with a cadence of Δ/N_p. This makes the high-cadence strategy extremely costly and unviable for the large N_p of O(10³) in an SKA-era PTA.

We use the following simulation setup to show that the staggered sampling schemes described above increase the frequency reach of a PTA. We consider a PTA with N_p = 50 nearest pulsars chosen from the simulated catalog in Smits et al. (2009). The total observation period is set at T = 5 yr with observations spaced Δ = 2 weeks apart in the case of uniform staggered sampling: this results in f_SP ≈ 4 × 10⁻⁷ Hz. (T is set lower than the typical value of 10 yr or more for PTA data to keep computational costs of the simulation in check.)We consider non-evolving sources with four angular frequencies ω_gw: 64 rad yr⁻¹ (3.23 × 10⁻⁷ Hz), 256 rad yr⁻¹ (1.29 × 10⁻⁶ Hz), 1024 rad yr⁻¹ (5.16 × 10⁻⁶ Hz), and 4096 rad yr⁻¹ (2.0656 × 10⁻⁵ Hz). Note that the last three sources have frequencies > f_SP and the highest one is very close to the staggered sampling limit of f_PTA = N_pf_SP = 4098.09 rad yr⁻¹. The sources are located at α = 3.5 rad and δ = 0.3 rad in equatorial coordinates. This location corresponds to the lowest degree of ill-posedness in parameter estimation for the SKA-era PTA used in Wang & Mohanty (2017). The inclination angle and the GW polarization angle are given by ι = 0.5 rad and ψ = 0.5 rad, respectively. The initial orbital phase is set at φ₀ = 2.89 rad.

Following the noise model described earlier, the standard deviation of the noise n^I(t) is set at σ^I = 100 ns. The overall amplitude, ζ, of the GW signal, which depends on the distance to the source, its chirp mass, and GW frequency, is determined by the specified network signal-to-noise ratio (S/N) ρ. Here ${\rho }^{2}={\sum }_{I=1}^{{N}_{{\rm{p}}}}{\left({\rho }^{I}\right)}^{2}$, where

$$\begin{eqnarray}&&{\left({\rho }^{I}\right)}^{2}=\displaystyle \frac{1}{{\left({\sigma }^{I}\right)}^{2}{\rm{\Delta }}}{\int }_{0}^{T}{\left({s}^{I}(t)\right)}^{2}{dt}\end{eqnarray} \tag{ 4 }$$

is the squared signal-to-noise ratio (S/N) of the GW-induced timing residual for the I-th pulsar.

Figure 1 shows the distributions of the detection statistic under the H₀ and H₁ hypotheses for both the uniform and the randomized staggered sampling strategies. From these, we estimate the detection probability of a ρ = 10 signal to be ≳90% at a false alarm probability of ≃1/500. The latter corresponds to setting the detection threshold at the largest value of the detection statistic obtained from H₀ data realizations. Within the precision of our simulation, these numbers are fairly independent of the staggered sampling strategy and the GW signal frequency. While a two-sample Kolmogorov–Smirnoff test on the H₀ distributions does show their apparent relative shift to be statistically significant, this has no noticeable effect on detection probability at the above S/N.

Zoom In Zoom Out Reset image size

Given that a higher frequency signal has a larger number of cycles in a given observation period, it is natural to ask if it can be detected over shorter observation periods using staggered sampling. We verified this by repeating the above simulations with T = 1 yr and ω_gw = 512 rad yr⁻¹ and 1024 rad yr⁻¹, keeping all else fixed. The detection probabilities had insignificant changes, suggesting that the performance of a staggered sampling-based search depends primarily on the S/N of a signal.

Figure 2 shows the estimated sky locations of the source for uniform staggered sampling (T = 5 yr) and a moderately strong S/N of ρ = 20. Within the precision of our simulations, the error in localizing a GW source does not show a clear trend with the frequency of the signal. Resolving a trend, if it exists, would require a computationally much more expensive simulation that we leave for future work. The typical localization error in Figure 2 of O(100) deg² makes searches for optical counterparts of GW sources feasible with the Rubin observatory (Ivezić et al. 2019) across the entire range, [T⁻¹, f_PTA], of GW frequencies (Liu et al. 2015; Wang & Mohanty 2017).

Zoom In Zoom Out Reset image size

For an SKA-era PTA with N_p = 10³ pulsars and per-pulsar observational cadence of 1/(2 weeks), our results show that staggered sampling can increase f_PTA = N_pf_SP to as high as 4 × 10⁻⁴ Hz. To estimate the achievable sensitivity, we use the same simulated PTA as in Wang & Mohanty (2017), comprising millisecond pulsars within 3 kpc taken from the synthetic catalog in Smits et al. (2009). From an analysis similar to Figure 1, we find that the detection probability of a monochromatic signal for the SKA-era PTA is ≃60% at ρ = 10 for a false alarm probability of ≃1/50. We adopt this as the fiducial value for the minimum detectable S/N averaged over the sky angles α and δ. (The resulting geometrical factor is ≃1 for this PTA.) Non-detection of a signal at this S/N will result in a sky-averaged upper limit on monochromatic GW strain amplitude,

$$\begin{eqnarray}&&h=8.89\times {10}^{-15}\left(\displaystyle \frac{{f}_{\mathrm{gw}}}{{10}^{-6}\,\mathrm{Hz}}\right){\left(\displaystyle \frac{T}{5\mathrm{yr}}\right)}^{-\tfrac{1}{2}}\left(\displaystyle \frac{\sigma }{100\,\mathrm{ns}}\right),\end{eqnarray} \tag{ 5 }$$

that is about three orders of magnitude lower in the [10, 400] μHz band than the current one from high-cadence observation of millisecond pulsar J1713 + 0747 (Dolch et al. 2016).

In the extended frequency range, not only would the inspiral phase of an SMBHB signal be observable but also the merger and ringdown phases. To quantify the sensitivity to each of these phases, we use the luminosity distance, D_L, for a sky-averaged S/N = 10. (The inclination and polarization angles are also averaged over in the case of inspirals). Because searches for the pulsar and Earth term can be decoupled for a strongly evolving signal (Finn & Lommen 2010), we make the conservative choice of using the S/N of only the Earth term.

The inspiral signal is calculated in the Newtonian approximation (Peters & Mathews 1963) over an observation period $\min \{20\,\mathrm{yr},\tau \}$, where τ is the lifetime for the signal frequency to evolve from an initial value f_i to f_ISCO, the frequency at the innermost stable circular orbit (ISCO). The merger and ringdown phases are obtained from waveforms computed in the spin-aligned effective one body numerical relativity for eccentric binary (SEOBNRE; Cao & Han 2017; Liu et al. 2020) formalism: the part of the waveform between the instantaneous frequency exceeding f_ISCO and the instantaneous amplitude attaining its maximum value is defined as the merger, with the subsequent phase being the ringdown. For the latter, only the dominant l = 2, ∣m∣ = 2 mode, with its corresponding frequency f_2,2, is used. We consider only circular binaries with zero spin and equal mass components for which the defining parameters are only the chirp mass ${{ \mathcal M }}_{c}=0.435\,M$ (M is the total mass) and, for the inspiral, the chosen τ.

Figure 3 shows D_L for all the different phases as a function of ${{ \mathcal M }}_{c}$ and τ along with their characteristic frequencies. We see that for ${{ \mathcal M }}_{c}\geqslant 4.5\times {10}^{9}$ M_⊙, the inspiral signal always stays below f_SP, irrespective of τ. On the other hand, the inspiral signal for ${{ \mathcal M }}_{c}\lt 4.5\times {10}^{9}$ M_⊙ would cross f_SP even if f_i < f_SP. Table 1 summarizes the distance reach, with the distance (D_L = 20 Mpc) to the Virgo cluster as a baseline, for different signal phases that require the extended frequency range (≥f_SP) of staggered sampling to be observable. Further applications of Figure 3 are considered below.

Zoom In Zoom Out Reset image size

The observation of ringdown signals by an SKA-era PTA with staggered sampling could extend the test of the no-hair theorem to the extremely large mass range of SMBHB remnants. For this we consider the test in Isi et al. (2019) that achieves an ≈10% level as defined by the fractional difference in the estimated values of a particular combination of mass and spin parameters measured from the late ringdown and the (S/N ≈ 14) post peak-amplitude waveforms of GW150914.

For the above S/N of the post peak-amplitude (our ringdown) waveform, staggered sampling will allow the ringdown from a ${{ \mathcal M }}_{c}\lesssim 2\times {10}^{10}$ M_⊙ system, for which f_2,2 ≳ f_SP, to be detected out to D_L ≲ 1.32 Gpc. Given that the corresponding inspiral signal would be extremely loud, S/N = 620 for τ > 2 yr, the source would be localized well in advance of the ringdown. This would allow a subset of favorably located pulsars to be targeted for significantly better timing over the duration of the ringdown. Assuming the timing residual noise is reduced from 100 ns to ≈20 ns (Feng et al. 2020), the observed ringdown S/N would increase to ≈70, leading to a test of the no-hair theorem at the ≈2% level.

Considering a lower mass system such as ${{ \mathcal M }}_{c}=5\times {10}^{8}$ M_⊙, Figure 3 shows that a τ = 5 yr inspiral signal, with f_i ≳ f_SP, will be detectable out to D_L ≈ 100 Mpc. While the corresponding ringdown signal (f_2,2 = 0.02 mHz) would be too weak for the SKA-era PTA, it would be extremely loud, with S/N ≈ 220, for a concurrently operating LISA. Compared to the fiducial ringdown S/N above, the relative measurement accuracy of all the ringdown parameters—inversely proportional to S/N from Fisher information analysis (Berti et al. 2006)—would improve by a factor of ≈16. In combination with the PTA-detected inspiral, this would again lead to a stringent test of the no-hair theorem.

In recent work (D'Orazio & Loeb 2020), a scheme for measuring the Hubble constant, H₀, has been proposed that uses only GW observations by an SKA-era PTA without requiring an electromagnetic counterpart. It relies on measuring both D_L and the co-moving distance D_c = D_L/(1 + z), z is the cosmological redshift, of a GW source through the effect of GW wave front curvature on pulsar timing residuals. The governing condition for the measurability of this effect is (D'Orazio & Loeb 2020) γ = πf_gwL²/D_c ≳ 0.1, where L is the Earth-pulsar distance. For the assumed GW frequency f_gw = 10⁻⁷ Hz < f_SP in D'Orazio & Loeb (2020), achieving a precision of δH₀/H₀ ≲ 30% in this scheme puts a rather stringent observational requirement on the error, δL, in L of δL/L ∼ 1% at L > 10 kpc. However, the feasibility of this scheme is improved if staggered sampling is used to reach higher f_gw. For example, L reduces to 3 kpc for the same relative error if f_gw = 10⁻⁶ Hz. This could happen if f_gw = f_ISCO for an ${{ \mathcal M }}_{c}\,=2\times {10}^{9}$ M_⊙ system, which would be detectable out to D_L = 2 Gpc (4 Gpc) for τ = 5 yr (20 yr), yielding γ = 0.21 (0.125).

The impact of the extended frequency reach from staggered sampling on the detectability of a wider range of signals than considered here needs further study. Among these are higher signal harmonics (Peters & Mathews 1963) from unequal-mass SMBHBs that, orbital evolution studies indicate (Sesana 2010), could be driven to high eccentricities (∼0.3) by interactions with the stellar environment. Independent evidence comes from observations (Dey et al. 2018) of the SMBHB candidate OJ 287 that suggest a binary mass ratio of ≃122 and eccentricity 0.657. Besides SMBHBs, oscillation of a network of cosmic strings (Burke-Spolaor et al. 2019), superradiance from axion clouds around isolated black holes (Cardoso et al. 2017), near-zone waves induced by turbulence of solar convection (Bennett & Melatos 2014), and solar oscillation modes (Cutler & Lindblom 1996; Polnarev et al. 2009) could be potential targets for a staggered sampling-based search.

Y.W. gratefully acknowledges support from the National Natural Science Foundation of China (NSFC) under grants No. 11973024 and No. 11690021, and Guangdong Major Project of Basic and Applied Basic Research (grant No. 2019B030302001). The contribution of S.D.M. to this paper is supported by NSF grant No. PHY-1505861. Z.C. gratefully acknowledges support from NSFC under grant No. 11690023. We thank Wen-Fan Feng for discussions on LISA. We acknowledge the Texas Advanced Computing Center (TACC) at the University of Texas at Austin (www.tacc.utexas.edu) for providing high-performance computing resources. We thank the anonymous referee for helpful comments and suggestions.