Cosmic variance of the Hellings and Downs correlation for ensembles of universes having nonzero angular power spectra

Gravitational waves (GWs) induce correlated perturbations to the arrival times of pulses from an array of galactic millisecond pulsars. The expected correlations, obtained by averaging over many pairs of pulsars having the same angular separation (pulsar averaging) and over an ensemble of model universes (ensemble averaging), are described by the Hellings and Downs curve. As shown by Allen [Phys. Rev. D 107, 043018 (2023)], the pulsar-averaged correlation will not agree exactly with the expected Hellings and Downs prediction if the gravitational-wave sources interfere with one another, differing instead by a "cosmic variance" contribution. The precise shape and size of the cosmic variance depends on the statistical properties of the ensemble of universes used to model the background. Here, we extend the calculations of the cosmic variance for the standard Gaussian ensemble to an ensemble of model universes which collectively has rotationally invariant correlations in the GW power on different angular scales (described by an angular power spectrum, $C_\ell$ for $\ell=0,1,\cdots$.). We obtain an analytic form for the cosmic variance in terms of the $C_\ell$'s and show that for realistic values $C_{\ell}/C_0\lesssim 10^{-3}$, there is virtually no difference in the cosmic variance compared to that for the standard Gaussian ensemble (which has a zero angular power spectrum).