Pulsar Timing Residuals

A pulsar timing array does not compare pulse times with a perfect clock prediction fixed in advance. It fits a timing model—spin, astrometry, binary motion, dispersion, instrumental offsets, and more—to measured times of arrival. Nanohertz gravitational waves appear in the residuals left after that fit, so any signal component resembling a timing-model basis function is partly absorbed. The likelihood must fit deterministic timing and stochastic processes together.

Helpful background. The Gravitational-Wave Likelihood supplies Gaussian-process likelihoods and marginalization. Pulsar Timing as a Gravity and Matter Probe develops the pulsar timing model.

From pulse profiles to residuals

Radio data are calibrated, dedispersed, folded, and compared with a template to estimate topocentric times of arrival (TOAs). Clock corrections and an adopted Solar-System ephemeris map them to the solar-system barycenter. A model $t_{\rm model}(\boldsymbol\beta)$ predicts pulse phase from parameters $\boldsymbol\beta$.

Linearizing around a reference solution gives prefit residuals

$$\delta\boldsymbol t =\boldsymbol M\,\delta\boldsymbol\beta +\boldsymbol n+\boldsymbol s,$$

where $M_{ia}=\partial t_i/\partial\beta_a$ is the design matrix, $\boldsymbol n$ collects measurement and stochastic noise, and $\boldsymbol s$ is a gravitational-wave or other common signal.

Generalized least squares projects out the column space of $M$. If the columns of $G$ form an ordinary orthonormal basis for $\ker M^T$, so that $G^TM=0$, the marginalized likelihood can be written in terms of

$$\boldsymbol r=\boldsymbol G^{T}\delta\boldsymbol t, \qquad C_r=\boldsymbol G^{T}C\boldsymbol G.$$

For a fixed design matrix, the compressed data $G^T\delta\boldsymbol t$ do not themselves depend on $C$; their probability distribution does. This formula assumes an effectively flat prior over the linear timing offsets and a full-rank design matrix. If $N$ TOAs fit $m$ independent timing parameters, $G$ has size $N\times(N-m)$ and a positive-definite $C$ gives an invertible $(N-m)\times(N-m)$ covariance $C_r$ van Haasteren and Levin 2013, §3.2, Eq. (15).

An equivalent noise-weighted implementation uses the generalized-least-squares projector $R=I-M(M^TC^{-1}M)^{-1}M^TC^{-1}$. Unlike $G$, $R$ depends on $C$; its $N\times N$ projected covariance $RCR^T$ is singular because $m$ timing directions have been removed. It must be evaluated in the surviving subspace, not naively inverted as a full covariance. Fitting once and ignoring this projection underestimates covariance with red processes.

What the timing fit absorbs

Spin frequency and spin-down remove low-order polynomials in time. Astrometric position and proper motion introduce annual and annually growing signatures. Parallax has a semiannual component. Binary orbital parameters produce their orbital harmonics. Telescope/backend jumps absorb constant offsets. Dispersion-measure variations produce chromatic delays approximately proportional to radio frequency $\nu^{-2}$.

A very-low-frequency gravitational-wave process can resemble a polynomial over a finite baseline; power near one year can covary with astrometry. The array’s response therefore has a transmission function determined by cadence, gaps, TOA errors, and the fitted basis. Sensitivity is not simply the Fourier transform of an evenly sampled record.

The observation span $T$ sets a characteristic lowest resolvable Fourier scale $1/T$, but fitting spin parameters further suppresses the lowest bins. Irregular sampling and changing receivers mix spectral power. Time-domain covariance or carefully constructed Fourier bases must include these effects.

Noise in one pulsar

White-noise parameters rescale formal TOA uncertainties and represent additional uncorrelated or epoch-correlated scatter. Red spin noise, dispersion variations, scattering, profile evolution, and band- or system-dependent processes can each be modeled as Gaussian processes with distinct chromatic scaling and spectra. A simple power-law residual spectrum is often written

$$P_a(f) =\frac{A_a^2}{12\pi^2} \left(\frac{f}{f_{\rm yr}}\right)^{-\gamma_a} f_{\rm yr}^{-3},$$

Here $f_{\rm yr}=(1\,\mathrm{yr})^{-1}$, $A_a$ is dimensionless, and $P_a$ has units $\mathrm{s}^3=\mathrm{s}^2/\mathrm{Hz}$. For a gravitational-wave process with characteristic strain $h_c=A(f/f_{\rm yr})^\alpha$, the usual PTA normalization gives $P_{\rm gw}=h_c^2/(12\pi^2f^3)$ and therefore $\gamma=3-2\alpha$. The same amplitude convention is often used for intrinsic red noise, without interpreting it as gravitational strain. A steep power law needs a low-frequency prescription; the finite Fourier basis and timing fit are part of its statistical definition.

A common process can be mistaken for single-pulsar red noise if cross-correlations are omitted. Conversely, underfit pulsar noise can inflate apparent common power. Noise models should be selected or marginalized without tuning them to erase the spatial signature being tested.

Dataset and reproduction boundary

The public NANOGrav 15-year release provides TOAs, timing solutions, narrowband/wideband products, and analysis-related data. Its detector-characterization paper documents receiver history, noise processes, timing treatment, and validation Agazie et al. 2023, §§2–6. This is the concrete dataset interface through the evidence cutoff 2026-08-18 NANOGrav, 15-year data portal.

A reproducible analysis pins .tim and .par file versions, clock and Solar-System ephemeris, TOA flags, observing systems, frequency channels, timing software, noise-model configuration, prior ranges, Fourier basis, sampler, and seeds. It verifies residual and design-matrix dimensions, reproduces single-pulsar noise posteriors, and injects known common processes before a multi-pulsar claim.

Diagnostic checks

• Refit after adding or removing plausible chromatic terms; gravitational-wave cross-correlation should not hinge on one pulsar’s dispersion model.
• Plot timing-model transmission functions and identify annual and lowest-frequency suppression.
• Compare narrowband and wideband timing where available.
• Perform leave-one-pulsar and leave-one-backend-out analyses.
• Recover simulated red processes with the real cadence and TOA uncertainties.
• Check residual normality and nonstationary epochs without assuming that a Gaussian process captures every anomaly.

Exercises

1. Spin fit. A signal is exactly $s(t)=a+bt+ct^2$ over the data span. What happens if phase, spin frequency, and spin-down are freely fitted?

Solution

Those timing parameters span constant, linear, and quadratic residual trends, so the signal lies in the design-matrix column space and is projected out in the linear approximation. Only departures from that polynomial or external constraints could retain information.

2. Chromatic discrimination. Why can multifrequency observations distinguish dispersion variation from an achromatic gravitational-wave residual?

Solution

Cold-plasma dispersion scales approximately as $\nu^{-2}$ with observing radio frequency, while a gravitational-wave shift in pulse arrival time is achromatic. Simultaneous frequency coverage can fit the chromatic dependence, though profile evolution and scattering require their own models.

References

• Agazie, Gabriella, et al. “The NANOGrav 15-Year Data Set: Detector Characterization and Noise Budget.” Astrophysical Journal Letters 951 (2023): L10. DOI. Open PDF.
• NANOGrav Collaboration. “Publicly Released Data Sets.” Accessed 2026-08-18. Data portal.
• van Haasteren, Rutger, and Yuri Levin. “Understanding and Analysing Time-Correlated Stochastic Signals in Pulsar Timing.” Monthly Notices of the Royal Astronomical Society 428 (2013): 1147–1159. DOI. Open PDF.