Limits to PTA sensitivity: spin stability and arrival time precision of millisecond pulsars

This paper assesses the quality of pulsars as astrophysical clocks for probing unmodeled processes, such as orbital motion and spacetime perturbations from gravitational waves (GWs). There are two broad questions that arise in this assessment:

• (i)  
  How stable are the spin rates of neutron stars (NSs)? and
• (ii)  
  How well can pulses be time tagged at an observatory?

The first issue relates to the quality of a NS as a clock while the second involves determination of when clock ticks (i.e. the pulses) are measured against an observatory clock. Once TOAs (times of arrival) are obtained, referencing them to the solar-system barycenter (SSBC) requires knowledge of astrometric parameters of the pulsar (location, proper motion, and parallax), which can be determined by nonlinear least-squares fitting to a multi-year sequence of TOAs or obtained through very-long-baseline interferometry. Pulsar timing relies on the following empirical foundation.

• (i)  
  Rotational stability allows individual turns of a NS to be predicted to high precision over time periods long enough (days to years) so that only infrequent monitoring is needed.
• (ii)  
  The emitted pulse is a cross-sectional slice through a radiation beam whose directionality is locked to the NS crust, at least in the mean.
• (iii)  
  Individual pulses vary in both phase (relative to a fiducial phase locked to the NS) and amplitude with stationary statistics, so that long time averages of N ≫ 1 pulses converge to a stable shape.
• (iv)  
  For objects with two or more stable shapes associated with metastable states of the magnetosphere, the shapes do not evolve with epoch.
• (v)  
  Frequency dependences of pulse shapes that derive from radiation physics and emission altitudes are regular and epoch independent.

Canonical pulsars (CPs; those with surface magnetic fields ∼10^{12 ± 1} G and periods P ≳ 10 ms) show large amounts of spin noise and metastable pulse shapes to varying degrees in different objects. However, millisecond pulsars (MSPs) are much more spin stable and have not as yet shown metastable profile changes. In particular, the rms perturbation from spin noise is ∼10⁻⁶ times smaller in the most spin-stable MSPs compared to typical CPs [17], allowing MSPs in pulsar timing arrays (PTAs) to probe processes well below the sub-μs level. Other effects intrinsic to emission from pulsar magnetospheres, from the interstellar medium (ISM) and other plasmas, and from instrumental effects provide a noise floor for timing precision that are discussed here. In the following, we use ν to represent electromagnetic frequency, f for GW frequency, and f_s for pulse rates of pulsars.

The topocentric TOA t of a pulse at radio frequency ν corresponds to a pulse phase

$$\begin{eqnarray} &&\phi (t, \nu ) = \phi _{\rm spin}(t) + \phi _{\rm mag}(t,\nu ) + \phi _{\rm orb}(t) + \phi _{\rm plasma}(t,\nu ) + \phi _{\rm ss}(t,\nu ) + \phi _{\rm tel}(t, \nu ). \end{eqnarray} \tag{ 1 }$$

Each term has a deterministic and a random component: ϕ_spin counts spin phase¹. ϕ_mag includes offsets from a fiducial point on the NS due to the variable locations and beaming directions of emission regions. Locations appear to be ν dependent and vary stochastically from motions within the overall particle beam volume. ϕ_orb includes motion of the NS from binary motion or from any orbital debris. ϕ_plasma includes ν-dependent propagation effects from dispersion, scattering, refraction, and birefringence in the ISM, the interplanetary medium (IPM), and ionosphere [8]. While the ISM produces the biggest effects, the IPM requires correction for objects viewed near the Sun. ϕ_ss describes orbital and gravitational propagation effects in the solar system. ϕ_tel includes contributions from radiometer noise and polarization of the antenna system. Terms in equation (1) are selectively discussed in the paper according to their relevance for determining current PTA sensitivity.

2.1. Minimum TOA error from matched filtering

The minimum TOA error is due to additive noise n_t that contributes to ϕ_tel in 1, from the receiver, ground pickup, and sky backgrounds (including the cosmic microwave background). Fitting a template U_t to a measured pulse $I_t = a + bU_{t-t_0} + n_t$, the TOA error is signal-to-noise ratio (S = b/σ_n) dependent. For U_t normalized to unit amplitude with N_ϕ samples at intervals Δt, the error from MF is

$$\begin{eqnarray} &&\sigma _{\rm t, MF} = \frac{\Delta t}{S \big[\sum\nolimits_{t=1}^{N_\phi -1} \left(U_t - U_{t-1}\right)^2\big]^{1/2}}. \end{eqnarray} \tag{ 2 }$$

We note that this is the minimum possible error for the TOA. For a Gaussian pulse with width W (FWHM) and Δt = P/N_ϕ (P = period),

$$\begin{eqnarray} &&\fl \sigma _{\rm t, MF} = \frac{\left(WP\right)^{1/2}}{(2\pi \ln 2)^{1/4} N_{\phi }^{1/2}S} \approx \frac{108 \;\rm ns}{S/100} \left[ \left( \frac{W} {100\ \rm \mu {\rm s}}\right) \left( \frac{P}{5\;\rm ms} \right) \left(\frac{2048}{N_{\phi }}\right)\right]^{1/2}, \end{eqnarray} \tag{ 3 }$$

where a fiducial pulse width of 100 μs and a period of 5 ms have been used. Typically, N ≫ 1 single pulses are averaged to increase S before TOAs are calculated. Overall, S∝(NΔν)^1/2/T_sys where Δν is the radiometer bandwidth and T_sys ≈ 20 to 50 K is the system temperature. For some MSPs, σ_{t, MF} as small as a few tens of ns have been obtained. Bandwidths of new spectrometers are about ten times larger than the previous generation, leading to MF fitting errors that are 1/3 of previous fitting errors. These do not necessarily mean that timing residuals will show this full benefit because other effects, such as pulse jitter and spin noise are unaffected by bandwidth improvements.

2.2. Neutron stars as clocks: spin noise

NSs spin down on time scales from ∼10³ to 10⁹ years due to radiation and particle emission from the magnetosphere. However, faster variations due to internal torques from crust-superfluid interactions and changes in magnetospheric torque are prominent in young, canonical pulsars and in magnetars, which have much higher fields. 'Glitches' are discontinuities in spin frequency and its time derivative (f_s or $\dot{f}_s$) with rise times ∼ seconds in many cases, and decay times from days to months, or longer. Stochastic noise in many CPs has nonstationary statistics associated with a 'red' power spectrum. Red noise is due to magnetospheric fluctuations that include jumps between discrete values of $\dot{f}_s$ [11] as well as unresolved variations in both f_s and $\dot{f}_s$. Fortunately, red noise has been identified in only a few MSPs and upper limits ≲ 100 ns over data spans T ≈ 5 yr have been placed on some objects. A global analysis of CPs and MSPs [17] shows broad consistency with rms TOAs (after removal of a second-order polynomial to account for f_s and $\dot{f}_s$) that scale as

$$\begin{eqnarray} &&\sigma _{\rm rn, 2} \approx 5^{+2.5}_{-1.7} \,\mu {\rm s} f_s^{-1.4\pm 0.1} \dot{f}_{s,-15}^{1.1\pm 0.1} T_{\rm yr}^{2\pm 0.2}, \end{eqnarray} \tag{ 4 }$$

for f_s in s⁻¹, $\dot{f}_s = 10^{-15}\,s^{-2}$ and T in yr. Object-to-object variations from this mean relation are large, about ±0.7 decade, and are not well understood. The corresponding spectrum is S_rn(f)∝f^{−5 ± 0.4}. This power law is problematic for detection of the stochastic background from supermassive black hole (SMBH) binaries which has a similar, nominal power-law spectrum ∝f^−13/3. The scaling law quantifies how MSPs are much better clocks owing to their very large spin rates and small $\dot{f}_s$. Current research aims to refine the characterization of red spin noise and to establish whether there are state changes in MSPs like those seen in CPs. If so, spin noise may be correctable [11].

2.3. Magnetospheric noise, pulse phase jitter, and beam evolution

Adding to torque variations that alter the spindown rate $\dot{f}_s$, individual pulses occur with stochastic offsets relative to a fiducial phase locked to the NS's spin. Offsets are comparable to the individual pulse widths and involve large (∼100%) amplitude variations causing $\sqrt{N}$ variations of a sample average of N pulses with respect to a long term average used as a template for MF. The strict conditions for MF are thus violated, leading to a TOA error that adds quadratically to the MF error given by

$$\begin{eqnarray} &&\sigma _J \approx \frac{f_JW}{2} \left(\frac{1+m_I^2}{2 N \ln 2}\right)^{1/2} \approx 20\ {\rm ns} \left( \frac{f_J}{1/3}\right) \left( \frac{W}{100\ \rm \mu {\rm s}}\right) \left( \frac{10^6}{N}\right)^{1/2} \left( \frac{1+m_I^2}{2}\right)^{1/2}, \end{eqnarray} \tag{ 5 }$$

where f_J is a dimensionless 'jitter' parameter (typically ∼1/3) equal to the rms time jitter of individual pulses in units of the width of the average pulse profile and m_I ∼ 1 is the modulation index (rms/mean) of pulse amplitudes. Pulse jitter differs qualitatively from torque fluctuations because it induces a TOA error even if there were no torque variations. The jitter produces a TOA error with white noise statistics for TOAs obtained from non-overlapping sets of pulses, in contrast to red-noise torque variations that have long correlation times.

A systematic TOA error results from changes in a pulsar's radiation beam with radio frequency, causing pulse shapes to vary in shape and in polarization. Emission varies across the beam owing to angular effects or to variations in emission altitude. Traditional MF using narrowband observations is now being replaced with wideband measurements that require attention to profile/beam evolution, as discussed in section Two-dimensional fitting in both time (or pulse phase) and frequency is necessary to avoid systematic errors and to achieve the theoretical minimum error.

2.4. Orbital perturbations

2.5. Interstellar medium effects: dispersion and scattering

Cold, magnetized plasma in the ISM has left-and-right hand circularly polarized normalized modes (LHCP, RHCP) with refractive indices $n_{l,r} \approx 1 - {\nu _p}^2/2\nu ^2 \mp {\nu _p}^2\nu _{\rm B\parallel }/2\nu ^3$, where n_e is the electron density, ν_p = (n_ee²/πm_e)^1/2 ≪ ν is the plasma frequency and ν_B∥ = eBcos θ/m_ec ≪ ν is the electron gyrofrequency, calculated for the magnetic field component along the line of sight (LOS); the ∓ cases apply for the two polarization modes. In the warm ionized ISM ν_p ≈ 1.56 kHz (n_e/0.03 cm⁻³)^1/2 and ν_B∥ ≈ 8.4 Hz(B/3 μG). The group velocity of a pulse yields frequency-dependent delays,

$$\begin{eqnarray} &&t = t_{\rm DM} \pm t_{\rm RM}, \end{eqnarray} \tag{ 6 }$$

where, the dispersive time delay and the small correction due to birefringence are

$$\begin{eqnarray} \fl t_{\rm DM} &=& \frac{{\rm e}^2}{2\pi m_e c} \frac{\int _0^D {\rm d}s\, n_e}{\nu ^2} = 4.15 \,{\rm ms}\, {\rm DM}\ \nu _{\rm GHz}^{-2}, \quad t_{\rm RM} = 2\,{\rm RM}c^2 \nu ^{-3} = 0.18 \,{\rm ns}\, {\rm RM}\nu _{\rm GHz}^{-3}, \end{eqnarray} \tag{ 7 }$$

where the dispersion and rotation measures and their standard units (for D in pc, n_e in cm⁻³, and B_∥ in μG) are

$$\begin{equation*} \fl {\rm DM}= \int _0^D {\rm d}s\, n_e(s) ({\rm pc\, cm^{-3}}),\quad {\rm RM}= \frac{{\rm e}^3}{2\pi m_e^2c^4} \int _0^D {\rm d}s \,n_e{\rm B_\parallel }= 0.81 \int _0^D {\rm d}s \,n_e{\rm B_\parallel }({\rm rad \,m^{-2}}). \end{equation*}$$

Phase perturbations from electron-density fluctuations on scales ≳ 10³ km are

$$\begin{eqnarray} &&\phi ({\bf x_{\perp }}, z) = \frac{2\pi }{\lambda }\int _0^D {\rm d}s\, n_{l,r} \approx -\lambda r_e \int _0^D {\rm d}z \,\delta n_e({\bf x_{\perp }}, z), \end{eqnarray} \tag{ 8 }$$

where r_e = e²/m_ec² is the classical electron radius and δn_e(x_⊥, z) is the electron-density variation versus perpendicular offset x_⊥ and line-of-sight coordinate z. Geometrical optics gives the refraction angle, θ_r = k⁻¹∇_⊥ϕ, while physical optics gives the diffraction angle, θ_d ∼ λ/ℓ_⊥, where ℓ_⊥ is the length scale on which the phase changes by ∼1 radian. For the ISM, θ_d ranges from <1 mas to ∼1 arc sec at 1 GHz for different lines of sight. Effects from multipath propagation include [15].

Variable dispersive delays. The dispersive timing delay is deterministic in frequency and is routinely removed to high precision using incoherent or coherent dedispersion algorithms [7], with DM typically known to one part in 10³. However, motions of the pulsar and ISM cause statistical variations in DM over days to years and therefore must be measured in order to dedisperse with high precision and to correct TOAs. Typical variations are σ_DM(T) ≈ 10⁻⁴ to 10⁻³ pc cm⁻³ over T = 1 yr [5, 10], scale as T^5/3, and have a power spectrum S_δDM∝f^−8/3 for a Kolmogorov spectrum (see below). The TOA variation has rms amplitude ∝σ_DM(T) and spectrum ∝S_δDM(f) as follows from equation (7).

Scattering delays. Pulse broadening from multipath propagation is described by convolution of an epoch and ν-dependent pulse broadening function (PBF) that is causal and whose width $\tau _d \sim D\theta _d^2/2c \propto \nu ^{-4}$ is approximately equal to the mean delay. τ_d is also strongly dependent on distance, so nearby MSPs (D ≲ 1 kpc) have small τ_d ≲ 10 ns at 1.4 GHz but can exceed 1 ms for MSPs with DM ≳ 300 pc cm⁻³ or at lower frequencies. Current PTAs therefore focus on low-DM MSPs but larger-MSP samples may include high-DM objects that will require careful removal of scattering delays. The wavenumber (q) spectrum for δn_e is broad, in many cases consistent with a Kolmogorov power law $S_{\delta n_e}(q) \propto q^{-11/3}$ [1]. In others, the spectrum has a different slope or corresponds to anisotropic irregularities. Larger scales that cause refraction also modulate the PBF, causing epoch dependence of τ_d. An active area of research involves algorithms for coherent descattering [4].

Refraction delays: Variable refraction causes two TOA variations. One is simply the geometric delay from angle-of-arrival (AOA) variations,

$$\begin{eqnarray} &&\sigma _{\rm AOA, geo} \sim \frac{D\theta _r^2(t)}{2c} \approx 1.2 \,\mu {\rm s}\, D(\rm kpc)\theta _r^2(\rm mas) \end{eqnarray} \tag{ 9 }$$

and another is from referencing TOAs to the SSBC,

$$\begin{eqnarray} &&\sigma _{\rm AOA, ssbc} \sim c^{-1} {\bf \hat{n}}\cdot {\bf r}_{\oplus } \cos \Phi (t) \approx 2.4 \,\mu {\rm s}\, \theta _r(\rm mas)\cos \Phi (t), \end{eqnarray} \tag{ 10 }$$

where r_⊕ = 1 AU and Φ(t) is the cyclical angle between the pulsar direction ${\bf \hat{n}}$ and the Earth-SSBC vector. Refraction angles ∝ν⁻² are expected to be ≲ 0.1 mas for nearby pulsars observed at 1.4 GHz, so the refractive TOA errors will be in the 10-100 ns range, but grow rapidly at lower frequencies and for more distant objects.

Intensity scintillations. Diffractive intensity scintillations (DISS) in time and frequency have characteristic scales Δt_d ∼ 100 sec and Δν_d ∼ 1 MHz that are highly dependent on ν, direction, source distance, and velocity. Patches of constructive interference ('scintles') in the ν–t plane have sizes ∼Δν_dΔt_d and a filling factor η_tη_ν ≈ 0.1. DISS is from diffraction of radiation into an angle θ_d ∼ λ/ℓ_d, where ℓ_d is the diffraction scale, the transverse length on which the phase perturbation from δn_e is one radian. Refractive intensity scintillation (RISS) is from large scale focusing and defocusing. The refraction scale ℓ_r ∼ Dθ_d, implying the well known relation, ℓ_rℓ_d ∼ λD ∼ (Fresnel  scale)² [15], and implies RISS time scales of days to years. RISS and DISS both modulate a pulsar's S/N, so the MF timing error can vary significantly. The finite number N_s of scintles across a bandwidth B and integration time T_int used to generate a TOA causes a $\sim \!\!1/\sqrt{N_s}$ variation in the shape of the PBF causing an additional TOA error $\sigma _{\rm t, scintles} \sim \tau _d / \sqrt{N_s}$. For nearby MSPs, N_s is small at 1.4 GHz but so too is τ_d, so this error is only a few ns but grows rapidly for higher DM pulsars or lower-frequency observations.

New spectrometers with nearly octave bandwidths will reduce radiometer noise significantly and allow ISM timing perturbations, all of which are strongly chromatic and with distinct wavelength dependences, to be more easily diagnosed and removed. Achieving rms timing residuals of tens of ns requires that scattering as well as dispersion effects be removed. The primary challenges of wideband spectrometers include the greater exposure to radio frequency interference and the intrinsic frequency dependence of pulse shapes that cause shifts between estimated pulse phases and the true spin phase of each NS. as discussed in section 2.3. We note that achromatic effects, including pulse jitter and spin noise, are not reduced with larger bandwidths; however lower noise will in fact reveal these achromatic effects in more MSPs.

GWs contribute two terms to timing residuals, the 'Earth' term that is correlated in the residuals for different pulsars and the 'pulsar' term that is uncorrelated along with other TOA perturbations. Stochastic backgrounds from SMBH binaries (e.g. [16]) are likely to be detected first in Earth terms by exploiting the angular correlation through a cross-correlation analysis. The same is true for continuous wave or chirped signals from individual SMBH binaries. GW bursts, however, may be detected first in a time-domain (rather than spectral-domain) analysis that reveals variance in one or more of the pulsar terms because the event rate of large bursts (such as bursts with memory, BWMs) is likely small and a PTA increases the time coverage by a factor equal to the number of MSPs, N_MSP. While BWMs may produce significant variance in the residuals of an MSP, it will be challenging to demonstrate that the variance is from a BWM [2].

Contaminating noise in timing residuals $\mathcal{R}(t)$ includes the pulsar term of the stochastic GW background along with other effects we have described. The composite spectrum for uncorrelated contributions to the residuals $\mathcal{R}(t)$ (i.e. excluding any GW signal terms) is

$$\begin{eqnarray*} &&S_{\mathcal{R}}(f) = S_{\rm gw}(f) + S_{\rm wn}(f) + S_{\rm rn}(f) + S_{\rm t, \delta {\rm DM}}(f) + S_{\rm AOA,geo}(f) + S_{\rm AOA, ssbc}(f) + S_{\rm \delta \tau _d}(f), \nonumber \end{eqnarray*}$$

that, respectively, are from the GW stochastic background (pulsar terms), white noise, red spin noise, DM variations, AOA variations (two terms), and variations in the DISS scattering delay. The last three terms are difficult to quantify but will have low-pass spectra. The gw, rn, and δDM terms are steep power laws that can be written in the form (s − 1)σ²(T)T^{−(s − 1)}f^−s with spectral indices s_rn = 5, s_gw = −13/3 for the stochastic background from SMBH binaries, and s_δDM = 8/3. Here the variance from a given spectral component in a time series of length T is σ²(T) and the expression explicitly shows the nonstationary statistics of the process. Variances for each term can be determined empirically or from models. For wn, the spectral form is $\sigma _{\rm wn}^2/(f_{\rm Ny}-1/T)$ where f_Ny is the Nyquist frequency for discrete sampling.

The residuals' noise spectrum $S_{\mathcal{R}}(f)$ is related to the commonly used h_c(f), the dimensionless strain per unit log frequency interval, and Ω_gw(f), the energy in GWs per unit log frequency relative to the cosmological closure density, according to [4, 9, 12, 14]

$$\begin{eqnarray} &&h_c^2(f) = 12\pi ^2 f^3 S_{\mathcal{R}}(f) = \frac{3H_0^2\Omega {\rm gw}(f)}{2\pi ^2 f^2}. \end{eqnarray} \tag{ 11 }$$

The GW contribution to the uncorrelated part of the noise amplitude (from the SMBH background contribution to the pulsar terms) is h_c(f) = Af^−2/3.

Table 1 summarizes the scaling laws for selected components of timing residuals and also gives some representative amplitudes and dependences on electromagnetic frequency ν for those effects that are chromatic.

Table 1. Scaling laws for contributions to RMS timing residuals of MSPs prominent contributions to timing residuals with typical RMS values corresponding to data spans T = 5 yr. Nonstationarities correspond to growth with increasing T. Spectral dependences are shown for the standard power spectrum $S_{\mathcal{R}}$ for the residuals $\mathcal{R}(t)$, for the dimensionless strain per unit log f, h_c(f), and for fractional cosmological energy density, Ω_gw(f) (see equation (11)).

	Residuals
Effect	Typical RMS	T Scaling	$S_{\mathcal{R}}(f)$	hc(f)	Ωgw(f)	Comment
White noisea	10 ns –10 μs	T0	f0	f3/2	f5	Smoothing:
						T−1/2
Spin noiseb	<100 ns–10 μs	T2 ± 0.2	f−5 ± 0.4	f−1 ± 0.2	f±0.4	CPs + MSPs
δDM(t)	1–10 μs	T5/6	f−8/3	f1/6	f7/3	ν−2, c
DISS: τd	1 ns–10 μs					ν−22/5, c
δAOA(t)
Geometric	≲ 1 μs	–	f−2/3	f7/6	f13/3	ν−4, c
SSBC	≲ 1 μs	–				ν−2, c
Stochastic GW	∼100 ns	T5/3	f−13/3	f−2/3	f2/3	SMBH
background						binaries

^aThree effects contribute to white noise: radiometer noise is uncorrelated in time and in radio frequency ν; pulse phase jitter that is uncorrelated in time but highly correlated across receiver bandwidths; and the finite scintle effect from DISS due characteristic correlation times and frequencies. ^bSpin noise is characterized from [17]. ^cPlasma effects are dependent on radio frequency, as indicated, and are correctable to varying degrees using multifrequency measurements. The quoted ranges of values in column 2 are given for typical frequencies (0.4–2 GHz).

Figure 1 shows the transmission function T(f) (left panel) that quantifies which frequencies are unaffected by the fit for spindown and astrometric parameters. The transmission function affects both the GW signal in the Earth terms and all sources of noise. A schematic form of the spectrum in the right-hand panel shows the noise contributions (i.e. no GW signal) to h_c(f) calculated from $S_{\mathcal{R}}(f)$ using only the wn, rn, and pulsar-GW terms. Values σ_wn = 30 ns, σ_rn = 20 ns(T/5 yr)² , and A = 10⁻¹⁵ at f = 1 yr⁻¹ have been used along with f_Ny = 20 yr⁻¹. Curves are shown for a single MSP and for a PTA comprising 50 MSPs, which will reduce the noise spectrum by a factor of $1/\sqrt{50}$. The curves represent the noise floor for detection of any correlated signals from the stochastic background or from continuous-wave signals. The apparent changes in slope at low frequencies are due to the effects of spindown fitting, as seen in the transmission function. Exclusion of the ISM terms is based on the following rationale. While DM variations contribute substantially with an h_c(f)∝f^1/6 scaling, current practice includes multiple radio frequency observations that allow estimation of DM(t) at each epoch and removal of the corresponding TOA variation. The refraction terms described earlier have ∼ν⁻² and ∼ν⁻⁴ dependences. The first of these is nearly degenerate with the δDM-induced variation while the second is nearly degenerate with pulse-broadening from multipath propagation. For detection of GWs at the tens of ns level, these degeneracies are probably adequate but for weaker signals, the specific differences between the refraction terms and the other terms may need to be addressed, particularly if future PTAs include MSPs that have higher DMs than in the current sample.

Zoom In Zoom Out Reset image size

4.1. Other effects

Current work includes efforts to make higher-order corrections of TOAs and torque variations by investigating pulse shape variations. 'Descattering' of pulse shapes aims to coherently deconvolve the interstellar PBF from the measured signal [4], which will remove the mean delay ∼τ_d and any epoch-to-epoch variations in τ_d. Some canonical pulsars show variations in torque that correlate with pulse shape changes on time scales of weeks to months. To varying degrees of precision, the torque changes can be removed from timing data [11]. Similar phenomena have not yet been seen in MSPs, either because the underlying physics is different or because MSPs have not been investigated on pulse-to-pulse time scales to the extent that CPs have. Lastly, amplitude and phase jitter of single pulses also is amenable to correction by exploiting any correlation between TOA variations and jitter-induced stochastic shape changes [6, unpublished notes]. To achieve significant improvement in TOA precision, the correlation coefficient ρ between a pulse shape parameter and TOA perturbations needs to be high. The reduction factor in rms TOA is $\sqrt{1-\rho ^2}$, so a 50% reduction requires $\rho = \sqrt{3}/2 = 0.87$ [6, unpublished notes]. Use of multiple shape parameters alters the expression but the basic point remains the same. Recent work [13] indicates ∼20% improvement and has been extended to ∼40% [20]. The degree of improvement depends on the statistics of jitter and seems to rely on an asymmetric phase-jitter distribution (unpublished work).

Current PTA limits on stochastic GW backgrounds approach predictions based on mergers of SMBH binaries yet are based largely on the previous generation of digital spectrometers. New, larger-bandwidth spectrometers combined with new algorithms will yield significantly better sensitivities and will provide constraining limits or detections of the stochastic background, of continuous wave sources (individual SMBH binaries or cosmic strings), or of bursts. The inherent nonstationarity of some of the GW signatures in timing residuals implies growth in signal levels as ∼T^5/3 (stochastic background) and ∼T for bursts with memory [2]. However, red spin noise also increases with data span (rms ∝T², cf equation (4)), so long-term projections depend on the uncertainty in the population average scaling law for spin noise and also on object-to-object differences. Nonetheless, analyses that combine measurements of multiple MSPs to detect the Earth term of the GW perturbation benefit from the $1/\sqrt{N_{\rm MSP}}$ reduction of noise levels. Based on several studies of detectability (Siemens, this volume) and [3], some general statements can be made about prospects for detection of a stochastic background. In most analyses, the detection significance is defined as $S_d = \hat{C} /\sigma _C$ where C is a cross correlation parameter that exploits the correlated Earth term and σ_C is its error. When σ_C is dominated by non-GW noise, S_d is said to be noise limited. In the opposite case, it is GW dominated.

• (i)  
  Increasing the number of MSPs (N_MSP) helps in any regime, but has greater impact in the noise limited case (rn or wn) where the detection significance S_d is linear in N_MSP. In the extreme case where the residuals are GW-dominated, S_d increases only as the $\sqrt{N_{\rm MSP}}$ [3]. Detection of GWs almost certainly will occur in the regime where S_d has a shallower than linear dependence on N_MSP.
• (ii)  
  In the wn-limited regime larger S_d can be obtained through a combination of more pulsars, as above, greater timing throughput (higher cadence and longer integration times), smoothing, and a decrease in timing error per TOA.
• (iii)  
  Because red noise is spectrally similar to the GW contribution and does not average out significantly in temporal correlation functions because of its long correlation time, the primary means for increasing S_d is to sum detection statistics over many pulsars.
• (iv)  
  In the GW-dominated regime increasing the number of pulsars N_MSP is the only option.

I thank D Madison for providing some of the code used to make figure 1 and to my NANOGrav and IPTA colleagues for many fruitful discussions. Comments from the referees improved the presentation. This work was supported by an NSF PIRE grant to Cornell University.