  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.
