Optimal robust detection statistics for pulsar timing arrays

Pulsar timing arrays (PTAs) seek to detect a nano-Hz stochastic gravitational-wave background (GWB) by searching for the characteristic Hellings and Downs angular pattern of timing residual correlations. So far, the evidence remains below the conventional $5$-$\sigma$ threshold, as assessed using the literature-standard ``optimal cross-correlation detection statistic''. While this quadratic combination of cross-correlated data maximizes the {\em deflection} (signal-to-noise ratio), it does not maximize the detection probability at fixed false-alarm probability (FAP), and therefore is not Neyman-Pearson (NP) optimal for the assumed noise and signal models. The NP-optimal detection statistic is a different quadratic form, but is not used because it also incorporates autocorrelations, making it more susceptible to uncertainties in the modeling of pulsar timing noise. Here, we derive the best compromise: a quadratic detection statistic which is as close as possible to the NP-optimal detection statistic (minimizing the variance of its difference with the NP statistic) subject to the constraint that it only uses cross-correlations, so that it is less affected by pulsar noise modeling errors. We study the performance of this new $\NPMV$ statistic for a simulated PTA whose noise and (putative) signal match those of the NANOGrav 15-year data release: GWB amplitude $A_{\rm gw}=2.1\times 10^{-15}$ and spectral index $\gamma=13/3$. Compared to the literature-standard ``optimal" cross-correlation detection statistic, the $\NPMV$ statistic increases the detection probability by $47\%$ when operating at a $5$-$\sigma$ FAP of $\alpha = 2.9 \times 10^{-7}$.