Marco Crisostomi, Rutger van Haasteren, Patrick M. Meyers, Michele Vallisneri

Pulsar Timing Array (PTA) searches for nHz gravitational-wave backgrounds (GWBs) typically model time-correlated noise by assuming a diagonal covariance in Fourier space, neglecting inter-frequency correlations introduced by the finite observation window. We show that this diagonal approximation can lead to biased estimates of spectral parameters, especially for the common red process that represents the GWB. To address these limitations, we present a method that (i) computes the time-domain autocorrelation on a coarse grid using a fast Fourier transform (FFT), (ii) interpolates it accurately to the unevenly sampled observation times, and (iii) incorporates it into a low-rank likelihood via the Sherman--Morrison--Woodbury identity. Using both analytic covariance comparisons and end-to-end simulations inspired by the NANOGrav 15-year dataset, we demonstrate that our method captures frequency correlations faithfully, avoids Gibbs ringing, and recovers unbiased spectral parameters with modest computational cost. As PTA datasets increase in sensitivity and complexity, our approach offers a practical and scalable path to fully accurate covariance modeling for current and future analyses.