Harmonic spectrum of pulsar timing array angular correlations

Pulsar timing arrays (PTAs) detect gravitational waves (GWs) via the correlations they create in the arrival times of pulses from different pulsars. The mean correlation, a function of the angle $\gamma$ between the directions to two pulsars, was predicted in 1983 by Hellings and Downs (HD). Observation of this angular pattern is crucial evidence that GWs are present, so PTAs "reconstruct the HD curve" by estimating the correlation using pulsar pairs separated by similar angles. The angular pattern may be also expressed as a "harmonic sum" of Legendre polynomials ${\rm P}_l(\cos \gamma)$, with coefficients $c_l$. Here, assuming that the GWs and pulsar noise are described by a Gaussian ensemble, we derive optimal estimators for the $c_l$ and compute their variance. We consider two choices for "optimal". The first minimizes the variance of each $c_l$, independent of the values of the others. The second finds the set of $c_l$ which minimizes the (squared) deviation of the reconstructed correlation curve from its mean. These are analogous to the so-called "dirty" and "clean" maps of the electromagnetic and (audio-band) GW backgrounds.