On the overlap reduction function of pulsar timing array searches for gravitational waves in modified gravity

Pulsar Timing Array (PTA) searches for gravitational waves (GWs) aim to detect a characteristic correlation pattern in the timing residuals of galactic millisecond pulsars. This pattern is described by the PTA overlap reduction function (ORF) \Gamma_ab(\xi_ab), which is known as the Hellings--Downs (HD) curve in general relativity (GR). In theories of modified gravity, the HD curve often receives corrections. Assuming, e.g., a subluminal GW phase velocity, one finds a drastically enhanced ORF in the limit of small angular separations between pulsar a and pulsar b in the sky, \xi_ab --> 0. In particular, working in harmonic space and performing an approximate resummation of all multipole contributions, the auto correlation coefficient \Gamma_aa seems to diverge. In this paper, we confirm that this divergence is unphysical and provide an exact and analytical expression for \Gamma_aa in dependence of the pulsar distance L_a and the GW phase velocity v_ph. In the GR limit and assuming a large pulsar distance, our expression reduces to \Gamma_aa = 1. In the case of subluminal phase velocity, we show that the regularization of the naive divergent result is a finite-distance effect, meaning that \Gamma_aa scales linearly with fL_a, where f is the GW frequency. For superluminal phase velocity (subluminal group velocity), which is relevant in the case of massive gravity, we correct an earlier analytical result for \Gamma_ab. Our results pave the way for fitting modified-gravity theories with nonstandard phase velocity to PTA data, which requires a proper understanding of the auto correlation coefficient \Gamma_aa.