On the calculation of p-values for quadratic statistics in Pulsar Timing Arrays

Pulsar Timing Array (PTA) projects have reported various lines of evidence suggesting the presence of a stochastic gravitational wave (GW) background in their data. One key line of evidence involves a detection statistic sensitive to inter-pulsar correlations, such as those induced by GWs. A $p$-value is then calculated to assess how unlikely it is for the observed signal to arise under the null hypothesis $H_0$, purely by chance. However, PTAs cannot empirically draw samples from $H_0$. As a workaround, various techniques are used in the literature to approximate $p$-values under $H_0$. One such technique, which has been heralded as a model-independent method, is the use of “scrambling” transformations that modify the data to cancel out pulsar correlations, thereby simulating realizations from $H_0$. In this work, scrambling methods and the detection statistic are investigated from first principles. The $p$-value methodology that is discussed is general, but the discussions regarding a specific detection statistic apply to the detection of a stochastic background of gravitational waves with PTAs. All methods in the literature to calculate $p$-values for such a detection statistic are rigorously analyzed, and many analytical expressions are derived. All this leads to the conclusion that scrambling methods are not model-independent and thus not completely empirical. Rigorous Bayesian and Frequentist $p$-value calculation methods are advocated, the evaluation of which depend on the generalized ${\chi }^2$ distribution. This view is consistent with the posterior predictive $p$-value approach that is already in the literature. Efficient expressions are derived to evaluate the generalized ${\chi }^2$ distribution of the detection statistic on real data. It is highlighted that no Frequentist $p$-values have been calculated correctly in the PTA literature to date.