A finite window couples Fourier bins
Pulsar timing array (PTA) searches look for nanohertz gravitational-wave backgrounds by modeling time-correlated noise as a Gaussian process. In the Fourier domain, the convenient assumption is that each frequency bin is independent: the covariance becomes diagonal, with power S(f_j)/T on the diagonal and zeros elsewhere.
The paper derives why that cannot be exactly right. When a stationary process is expanded in a Fourier series over a finite interval of length T, the covariance between coefficients at frequencies f_j and f_k is an integral of the true spectrum weighted by two sinc functions, sinc(pi T (f - f_j)) and sinc(pi T (f - f_k)). Those sinc factors overlap across neighboring bins, so the Fourier coefficients are correlated even for a perfectly stationary process. The diagonal form in the standard analysis simply drops these off-diagonal terms.
Shallow red-noise spectra show the largest shifts
The authors do not claim the diagonal approximation is always wrong. Their analytic covariance comparisons and end-to-end simulations, inspired by the NANOGrav 15-year dataset, indicate that for smooth, low-frequency-dominated spectra the diagonal approximation is adequate, particularly after timing-model projection.
The failure mode is sharper. For spectra with sharp in-band features, the diagonal Fourier-domain prior cannot represent power falling between the bin frequencies k/T, and errors reach order unity even after projection. In the simulated dataset, the approximation also induces small errors in individual pulsars' timing-noise parameters for shallow-spectrum pulsars. Notably, in the realization studied it did not produce discernible method-dependent shifts in the common-process parameters, so the effect is method- and parameter-dependent rather than a uniform bias.
FFT on a coarse grid, then local interpolation
The proposed method replaces the diagonal Fourier prior with a covariance built in the time domain. It first evaluates the Wiener-Khinchin integral as a Fourier sum on a coarse, regularly spaced grid, using an inverse fast cosine transform. Two knobs control accuracy and cost: an oversampling factor omega that increases frequency resolution, and a factor zeta that integrates up to a multiple of the coarse grid's Nyquist frequency.
The resulting autocorrelation vector is then interpolated onto the unevenly sampled times of actual time-of-arrival measurements. Each measurement is expressed as a linear combination of only the two neighboring coarse-grid nodes that bracket it, so the interpolation stays local. Gaps in the data are handled naturally: the coarse grid spans the full observation window, and correlation across a gap is carried by the autocorrelation at that lag. The same construction extends to chromatic processes such as dispersion-measure variations, where the row weight picks up a (nu/nu_ref)^(-alpha) factor.
Local interpolation avoids Fourier edge ringing
A central comparison in the paper is between two rank-reduced representations of the same coarse signal. On the coarse node grid, the FFT-interpolation approach and a non-diagonal Fourier (SincPhi) approach agree exactly at the nodes; both are complete representations there.
They diverge when upsampled back to the full fine grid. The Fourier basis is non-local, so evaluating it on the full grid produces the Gibbs phenomenon, a ringing artifact concentrated at the edges of the observation window and spreading across the domain. The paper illustrates this with a signal defined on 10,000 fine points over x in [0, 8000], observed only within [2000, 6000], using 61 coarse nodes corresponding to 30 Fourier modes. The time-domain interpolation follows the original signal without artifacts, while the Fourier reconstruction rings. The authors note that this edge effect would affect the unevenly spaced real PTA data, even though it does not affect the coarse basis itself.
Spectral lines expose the diagonal prior
To probe the methods beyond a single realistic simulation, the authors use a family of Gaussian power spectral densities with analytic time-domain covariance, parameterized by amplitude, peak frequency mu_f, and width sigma_f. This gives an exact reference against which each approximation can be scored.
The results separate the methods by spectral shape. The diagonal prior breaks down for spectral lines and other in-band features narrower than the frequency-bin width 1/T, with errors of order unity for both tested node counts. The SincPhi method captures such features considerably better but remains limited in accuracy, especially at the smaller node count. The FFT-interpolation method is the one the authors present as capturing frequency correlations faithfully. A separate covariance comparison across Matérn 3/2 and broken power-law kernels, with constant and quadratic timing-model projection, shows the mean absolute element-wise error decreasing as the number of coarse timestamps and the oversampling factor increase.
Common-process posteriors stay stable in this realization
The end-to-end test uses a simulated NANOGrav 15-year-like dataset and compares posteriors for the common red-noise parameters gamma_gw, log10 A_gw, and log10 f_b. The diagonal runs use n = 30, 80, and 120 Fourier modes; the new method uses N-hat = 61, 161, and 241 coarse nodes, with a direct comparison at N-hat = 241 versus n = 120.
The honest reading is narrower than a blanket claim of bias. The paper reports that the diagonal approximation produces method-dependent shifts in spectral-parameter estimates especially for shallow intrinsic red-noise spectra, and small errors in individual pulsar timing-noise parameters in this simulation. But for the common-process parameters in this particular realization, it finds no discernible method-dependent shifts. The case for the new method therefore rests on the analytic failure for in-band features, the edge behavior, and the expectation that tolerance for mismodeling shrinks as datasets grow, not on a demonstrated common-process bias in this one dataset.
Accuracy matters more as PTA sensitivity grows
The paper's practical priority follows from its own scenario split. For smooth, low-frequency-dominated spectra, the diagonal approximation remains adequate, particularly after timing-model projection, so switching methods is not urgent on accuracy grounds alone. For spectra with sharp in-band features, or for shallow-spectrum pulsars whose individual noise parameters are being inferred, the diagonal prior's order-unity errors and parameter shifts matter now.
The forward-looking argument is about scale. As datasets grow in size and sensitivity, through consolidation across multiple PTAs and next-generation facilities, the tolerance for covariance mismodeling shrinks. The authors position their method as a computationally tractable route to fully accurate covariance modeling in that regime, with the cost controlled by the coarse-grid size and oversampling choices rather than by inverting a large dense covariance matrix. The decision threshold the paper supports is therefore spectral shape and target parameter: use the diagonal shortcut where the spectrum is smooth and low-frequency-dominated, and adopt the correlation-aware construction where in-band structure or shallow red noise is what the inference depends on.
Sources used in this answer
Beyond diagonal approximations: improved covariance modeling for pulsar timing array data analysis
A diagonal Fourier-domain covariance neglects inter-frequency correlations from the finite observation window, producing method-dependent spectral-parameter shifts especially for shallow red-noise spectra; an FFT-plus-interpolation low-rank method captures these correlations, avoids Gibbs ringing, and gives estimates closer to injected values in the studied realization at modest cost.