Cosmic Variance in Anisotropy Searches at Pulsar Timing Arrays

TL;DR

This paper investigates how cosmic variance affects anisotropy measurements in pulsar timing array data, demonstrating that statistical inference remains robust and can distinguish between astrophysical and cosmological gravitational wave backgrounds.

Contribution

It provides a detailed analysis of cosmic variance impacts on PTA anisotropy searches and confirms the robustness of statistical inference in this context.

Findings

Anisotropy measurements are robust despite cosmic variance.

Maximum multipole scales linearly with the number of pulsars.

PTA data can effectively discriminate between astrophysical and cosmological GWBs.

Abstract

Recent pulsar timing array (PTA) analyses show evidence for a gravitational wave background (GWB) with angular correlations consistent with the Hellings-Downs curve. Anisotropies are a key discriminator of the origin of this GWB, as they are expected to be at 1--20\% for astrophysical sources, but suppressed for cosmological GWBs. However, contrary to gravitational wave detectors at higher frequencies, PTAs only take a few independent measurements of a GWB and consequently are highly sensitive to cosmic variance, which induces apparent anisotropies in individual realizations of an isotropic GWB. We demonstrate explicitly that statistical inference nevertheless remains robust, i.e., measurements are consistent with the underlying assumption of isotropy. This confirms that searches for anisotropies will be able to robustly discriminate astrophysical from cosmological GWBs. En route, we demonstrate that the maximum multipole constrained by a PTA dataset scales linearly with the number of pulsars $\ell_{\rm max} \sim N_p$.

Equations

$$\delta t_I=\frac{{\hat{p}}_I^a{\hat{p}}_I^b}{2}{\int }_0^{D_I}\mathrm{d}s{h}_{ab}\left(t(s),\vec{x}(s)\right),$$

$$h_{ab}(t,\vec{x})$$

$$\times \left[{\tilde{h}}_{+}(f,\hat{\Omega })e_{ab}^{+}(\hat{\Omega })+{\tilde{h}}_{\times }(f,\hat{\Omega })e_{ab}^{\times }(\hat{\Omega })\right],$$

$$⟨{\tilde{h}}_P(f,\hat{\Omega }){\tilde{h}}_{P^{\prime }}^{\ast }(f^{\prime },{\hat{\Omega }}^{\prime })⟩=\frac{\mathcal{P}(f,\hat{\Omega })}{4}{\delta }_{P{P}^{\prime }}{\delta }^2(\hat{\Omega },{\hat{\Omega }}^{\prime })\delta (f-f^{\prime }),$$

$$\mathcal{P}(f,\hat{\Omega })\equiv S_h(f)P(\hat{\Omega }),$$

$${\Omega }_{\mathrm{GW}}{h}^2=\frac{h^2}{{\rho }_c}\frac{\mathrm{d}{\rho }_{\mathrm{GW}}}{\mathrm{d}\log f}=\frac{2{\pi }^2{f}^3}{3H_0^2/h^2}{S}_h,$$

$$P(\hat{\Omega })=\sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{c}_{\ell m}{Y}_{\ell m}(\theta ,\varphi ),$$

$$D_{IJ,ij}={\widetilde{\delta t}}_{I,i}{\widetilde{\delta t}}_{J,j}^{\ast },$$

$${\widetilde{\delta t}}_{I,i}\equiv \frac{1}{i2\pi }\sum _P{\int }_{S^2}{\mathrm{d}}^2\hat{\Omega }{\int }_{-\infty }^{\infty }\frac{\mathrm{d}f}{f}{\tilde{h}}_P(f,\hat{\Omega })R_I^P(f,\hat{\Omega })$$

$$\times \text{sinc}[\pi (f_i-f)T_{\text{obs}}],$$

$$R_I^P(f,\hat{\Omega })$$

$$C_{IJ,ij}$$

$$\times \text{sinc}[\pi (f_i-f)T_{\text{obs}}]\text{sinc}[\pi (f_j-f)T_{\text{obs}}],$$

$${\Gamma }_{IJ}(f)\equiv {\int }_{S^2}{\mathrm{d}}^2\hat{\Omega }{\kappa }_{IJ}(f,\hat{\Omega }){\gamma }_{IJ}(\hat{\Omega })P(\hat{\Omega }),$$

$${\gamma }_{IJ}(\hat{\Omega })\equiv$$

$$-\frac{3}{8}(1-{\hat{p}}_I\cdot \hat{\Omega })(1-{\hat{p}}_J\cdot \hat{\Omega }),$$

$${\kappa }_{IJ}\equiv \left(1-e^{-2\pi if{D}_I(1+{\hat{p}}_I\cdot \hat{\Omega })}\right)\left(1-e^{2\pi if{D}_J(1+{\hat{p}}_j\cdot \hat{\Omega })}\right).$$

$${\kappa }_{IJ}(f,\hat{\Omega })\simeq 1+{\delta }_{IJ}.$$

$$C_{IJ,ij}=C_{IJ}^p\otimes C_{ij}^f.$$

$${\Gamma }_{IJ}=\sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{c}_{\ell m}{\Gamma }_{IJ,\ell m}.$$

$${\widetilde{\delta t}}_{I,i}$$

$$s_{ik}^{\pm }\equiv \mathrm{sinc}(\pi (f_i\pm f_k)T_{\mathrm{obs}}),$$

$$D_{IJ,ij>0}=2\mathrm{Re}\left[{\widetilde{\delta t}}_{I,i}{\widetilde{\delta t}}_{J,j}^{\ast }\right].$$

$$⟨{\tilde{h}}_{P,n,k}^{\ast }{\tilde{h}}_{P^{\prime },n^{\prime },k^{\prime }}⟩$$

$$⟨{\tilde{h}}_{P,n,k}{\tilde{h}}_{P^{\prime },n^{\prime },k^{\prime }}⟩$$

$$C_{IJ,ij>0}=\frac{4}{3}{\Gamma }_{IJ}$$

$${\Gamma }_{IJ}=C_{IJ}^p=\frac{3}{2}\Delta \hat{\Omega }\sum _{P,n}{H}_n{R}_{I,n,k}^P{R}_{J,n,k}^{P\ast ,},$$

$$-\ln \mathcal{L}(\widetilde{\delta t}|\mathbf{\theta })=\sum _{IJ,ij>0}\ln [\det C_{IJ,ij}]+{\widetilde{\delta t}}_{I,i}{C}_{IJ,ij}^{-1}{\widetilde{\delta t}}_{J,j}^{\ast },$$

$$\displaystyle\delta{\hat{\bm{\theta}}}_{\alpha}^{(i)}=\sum_{\beta}[F^{-1}(\bm{\theta})]_{\alpha\beta}\frac{\partial\ln{\cal L}(\tilde{\delta}t|\bm{\theta})}{\partial\bm{\theta}_{\beta}}\Bigg{|}_{\bm{\theta}={\hat{\bm{\theta}}}^{(i)}}\,,$$

$$\frac{\partial \ln \mathcal{L}(\tilde{\delta }t|\mathbf{\theta })}{\partial {\mathbf{\theta }}_{\alpha }}=\mathrm{Tr}\left[C^{-1}\frac{\partial C}{\partial {\mathbf{\theta }}_{\alpha }}{C}^{-1}(D-C)\right],$$

$$F_{\alpha \beta }=\mathrm{Tr}\left[C^{-1}\frac{\partial C}{\partial {\mathbf{\theta }}_{\alpha }}{C}^{-1}\frac{\partial C}{\partial {\mathbf{\theta }}_{\beta }}\right].$$

$$C_{\ell }\equiv \sum _m\frac{c_{\ell m}^2}{2\ell +1}.$$

$$c(f)\sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{c}_{\ell m}{\Gamma }_{IJ,\ell m}=⟨{\tilde{d}}_I{\tilde{d}}_J⟩,$$

Full text

††preprint: CERN-TH-2025-159

**Introduction. ** Pulsar Timing Arrays (PTAs) monitor the arrival times of pulses from millisecond pulsars, whose intrinsic stability allows for precise predictions of their time-of-arrival at Earth. Gravitational waves (GWs) passing between Earth and these pulsars induce characteristic deviations, known as timing residuals, which are predictably correlated across pulsars. For an isotropic, gravitational wave background (GWB), the ensemble-averaged correlations trace the Hellings-Downs (HD) curve Hellings:1983fr , providing a distinctive signature for GW detection.

Recent results from worldwide PTA collaborations have reported evidence for a stochastic signal consistent with nanohertz GWs, displaying angular correlations in line with the HD expectation NANOGrav:2023gor ; Antoniadis:2023ott ; Reardon:2023gzh ; Xu:2023wog . The precise origin of this background remains uncertain. While a GWB from supermassive black hole binaries (SMBHBs) is the leading interpretation, alternative scenarios involving a cosmological GWB from the early Universe remain compatible with current data Madge:2023cak ; NANOGrav:2023hvm ; Antoniadis:2023xlr ; Figueroa:2023zhu ; Ellis:2023oxs .

A promising avenue to distinguish these possibilities lies in the search for anisotropies in the GWB. Most cosmological sources predict highly isotropic backgrounds, whereas SMBHBs are expected to generate anisotropies at the $1-20\%$ level at low multipoles Mingarelli:2013dsa ; Taylor:2013esa ; Mingarelli:2017fbe ; NANOGrav:2023tcn ; Sah:2024oyg ; Sah:2024etc . Thus, the detection of such anisotropies would provide strong evidence for an astrophysical origin of the signal, motivating the need to assess the sensitivity of current and future PTA datasets to these features.

In this work, we clarify the role of cosmic variance in PTA anisotropy measurements. Given their intrinsic frequency and angular coverage compared to the correlation length and time scales of the GWBs probed, PTAs are only sensitive to a few effective realizations of the GWB, unlike GW observatories at higher frequencies. Even in the limit of very many isotropically distributed pulsars and a GWB arising as the superposition of very many independent sources (e.g., of cosmological origin), this implies that the output of PTA measurements depends on the realization of the GWB in the local Universe, referred to as cosmic variance.111 Our terminology of cosmic variance is in analogy to CMB measurements and is referred to as sample variance in Nay:2023pwu ; Pitrou:2024scp . It refers not only to the variance of the HD curve, but to all PTA observables. Additional scatter in the measurements arises from noise, from the finite number of pulsars in a PTA (also referred to as pulsar variance Allen:2022dzg ) and, for a GWB due to SMBHBs, from their distribution in the sky (also referred to as shot noise Grimm:2024lfj ; Cusin:2025xle ). These effects arise in addition to the irreducible cosmic variance discussed here and are not the subject of this work. In particular, this implies that locally measured anisotropies may significantly differ from ensemble averages. Additionally, cosmic variance induces mild variations in the HD correlation itself, potentially influencing the sensitivity to anisotropies Agarwal:2024hlj .

Using realistic data simulations, we explicitly show that although maximum likelihood estimators obtained from individual realizations of an isotropic GWB can exhibit deviations due to cosmic variance, the recovered uncertainties render the inference statistically consistent with the underlying isotropic nature of the background, even in the signal-dominated regime. This in particular implies that cosmic variance cannot ‘fake’ anisotropies in a statistically significant way, and that the discrimination of astrophysical and cosmological GWBs based on anisotropies is robust. These results can be reproduced with fastPTA, a publicly available Python code to forecast the constraining power of PTA configurations.

**GWB anisotropies in PTA data. ** PTAs measure the timing residual $\delta t_{I}$ induced by GWs along the line-of-sight between a pulsar $I$ (located in direction $\hat{p}_{I}$ at distance $D_{I}$) and Earth (located here at position $\vec{0}$). The time shift induced at time $t$ then reads

[TABLE]

with $s$ denoting the affine parameter parameterising the geodesic, $t(s)=t-(D_{I}-s)$, $\vec{x}=(D_{I}-s)\hat{p}_{I}$, and $a,b$ the spatial components of the vectors. In the Fourier domain, the spin-2 perturbations of the metric read

[TABLE]

where we introduced the polarization tensors $e^{+,\times}_{ab}(\hat{k})$ with $P=\{+,\times\}$. Reality of $h_{ab}$ imposes $\tilde{h}_{P}(f,\hat{\Omega})=\tilde{h}_{P}^{*}(-f,\hat{\Omega})$, with $\hat{\Omega}$ denoting the normalized gravitational wave vector. For a stationary and unpolarized GWB, we can express the ensemble average as

[TABLE]

where we take ${\cal P}(f,\hat{\Omega})$ to be a deterministic function of the model. In this work, we assume the power spectral density (PSD) to be factorizable as

[TABLE]

with normalization condition $\int{\rm d}^{2}\hat{\Omega}P(\hat{\Omega})=1$. In general, the factorization in Eq. (4) may not hold, as for the case of GWB produced by SMBHBs, for which anisotropies become larger at large frequencies Becsy:2022pnr ; Sato-Polito:2023spo ; Lemke:2024cdu ; Gersbach:2024hcc ; Raidal:2024tui ; Moreschi:2025qtm . However, this assumption does not affect our argument, and we will come back to this point in the conclusions.

The GW abundance today reads

[TABLE]

where $\rho_{c}$ is the critical energy density of the Universe and $H_{0}/h=1/(9.78\,{\rm Gyr})$ is the Hubble parameter today. Building on the framework developed in previous studies NANOGrav:2023tcn ; Allen:1996gp ; Anholm:2008wy ; Mingarelli:2013dsa ; Taylor:2013esa ; Taylor:2015udp ; Taylor:2020zpk ; Ali-Haimoud:2020ozu ; Ali-Haimoud:2020iyz ; Pol:2022sjn , we expand the GWB power spectrum in spherical harmonics and combine it with the geometric configuration of PTAs to predict the induced timing residual correlations. To account for an anisotropic GWB, we expand $P(\hat{\Omega})$ as

[TABLE]

where $Y_{\ell m}(\theta,\phi)$ are the real-valued spherical harmonics. We adopt the same notation used in Depta:2024ykq . The normalization condition we use imposes $c_{00}=1/\sqrt{4\pi}$.

Given $N_{p}$ pulsars, the maximum number of multiples that can be constrained scales as $\ell_{\rm max}\propto N_{p}$ (not its square root, as reported in Romano:2016dpx ; NANOGrav:2023tcn ). The argument proposed in Romano:2016dpx is based on the assumption that from each pulsar data stream one can extract at most one complex number. This idea resonates with CMB temperature maps, where each pixel brings information about temperature or intensity of that portion of the sky alone (neglecting polarization, for simplicity). However, in the PTA case, each additional pulsar adds $N_{p}-1$ pulsar pairs, or new detectors, which are sensitive to the full sky through a unique pattern (given by the response function) that is not fully degenerate with the remaining pulsar pairs. We elaborate on this point in App. A.

**Frequency domain correlation function. ** The information on the GWB (anisotropies) is encoded in the correlations in the time delays measured between pulsars $I$ and $J$ at frequencies $i,j$, and can be expressed as

[TABLE]

where we introduced the Fourier space time delays

[TABLE]

and the response function

[TABLE]

The ensemble average of $D_{IJ,ij}$ can be derived using the statistical properties of the GWB (2), which gives

[TABLE]

in terms of the overlap reduction function (ORF)

[TABLE]

with the pairwise timing response function $\gamma_{IJ}(\hat{\Omega})$ Ali-Haimoud:2020ozu

[TABLE]

and

[TABLE]

For PTAs, the minimum frequency is limited by the observation time $f\sim 1/T_{\rm obs}\sim\text{nHz}\sim 0.1/\text{pc}$, while the typical distance to the pulsars is $D_{I}\sim~\text{kpc}$. Therefore, $fD_{I}\gg 1$, so that the rapidly oscillating pieces in Eq. (13) are negligible when integrated over the directions $\hat{\Omega}$ (see also Mingarelli:2013dsa ), and $\kappa_{IJ}$ reduces to

[TABLE]

This simplification also implicitly assumes that fitting to the timing models reduces power below $1/T_{\rm obs}$, see, e.g., Refs. Blandford:1984xwb ; Hazboun:2019vhv ; Babak:2024yhu ; Pitrou:2024scp , for related discussion. As a further test for the validity of this approximation, we will not adopt it for data generation, but only in the likelihood evaluation. As a consequence of this and the initial assumption (4), the correlation matrix is factorizable

[TABLE]

Substituting (6) into (11), and using the orthogonality of $Y_{\ell m}(\theta,\phi)$, the ORF reads

[TABLE]

The quantities $\Gamma_{IJ,\ell m}$, referred to as generalized ORFs Mingarelli:2013dsa , encode the geometric structure of the PTA experiment. They can be viewed as a generalization of the HD function $\chi_{IJ}\equiv c_{00}\Gamma_{IJ,00}$ Hellings:1983fr , see Fig. 1 (which is only a function of the angle $\cos\zeta_{IJ}=\hat{p}_{I}\hat{p}_{J}$ between pulsars $I$ and $J$ as dictated by symmetries, see, e.g., Kehagias:2024plp ), to higher multipoles.

**Data simulation. ** We start from Eq. (8) and turn the Fourier space integration into a sum over a discrete set of frequencies $N_{f}^{\rm int}$, as well as the angular integration over a pixelized sky with $N_{\rm pix}$ pixels to obtain

[TABLE]

where we introduced

[TABLE]

and the pixel index $n$ labels directions on the sky using an equal-area tessellation with angular size $\Delta\hat{\Omega}=4\pi/N_{\rm pix}$ obtained with the HEALPix package Gorski:2004by , while the index k labels positive only frequencies $f_{k}=f_{\rm min}+k\Delta f$, with $k=0,1,...,N_{f}^{\rm int}$. We take $f_{\rm min}=1/2T_{\rm obs}$ and $\Delta f=1/10\,T_{\rm obs}$222The low frequency integration cut-off is motivated by the effective sharp reduction of power at low frequencies induced by fitting the timing model  Blandford:1984xwb ; Hazboun:2019vhv ; Babak:2024yhu ; Pitrou:2024scp ; Allen:2024bnk . We tested that our choice does not significantly affect our results. Choosing a lower $f_{\rm min}$ and very red spectra would require going beyond the approximation (14) in the likelihood, as corrections proportional to $f_{\rm min}D_{I}$ would become relevant for large $N_{p}$ networks.. We perform the analysis for a base of positive frequencies $f_{i}=i/T_{\rm obs}$, with $i>0$, therefore

[TABLE]

To generate GWB realizations, we simulate the complex coefficients $\tilde{h}_{P,n,k}$ appearing in Eq. (17) as Gaussian random variables with zero mean and two-point function

[TABLE]

This equation is the discretized version of Eq. (3), with $H_{n,k}=P(\hat{\Omega}_{n})S_{h}(f_{k})/4$, and $H_{k}=\sum_{n}H_{n,k}=S_{h}(f_{k})/4$ is directly related to the sky-integrated strain PSD in the $k$-th frequency bin333Notice that $\Delta\hat{\Omega}$ in Eq. (20) ensures that the total power in a given sky patch is pixelization-independent.. When assuming anisotropic injections with $c_{\ell>0,m}\neq 0$, we distribute power across pixels according to (6). Because these variables are complex and Gaussian, we take the real and imaginary parts of $\tilde{h}_{P,n,k}$ as $\tilde{h}_{P,n,k}=(h^{r}_{P,n,k}+ih^{i}_{P,n,k})/\sqrt{2}$, where $h_{r}$ and $h_{i}$ are independent real Gaussian random variables with zero mean and variance $\sigma_{n,k}^{2}=H_{n,k}/{\Delta\hat{\Omega}\Delta f}$.

The expectation value of $D_{IJ,ij}$ (10) becomes

[TABLE]

with

[TABLE]

where, assuming the factorization (4), i.e. $H_{n,k}=H_{n}H_{k}$, holds, (14) implies the l.h.s. to be frequency independent.

In our forecasts, for each realization, we generate pulsar positions randomly on the sky, assuming a uniform distribution over the sphere. We illustrate these realizations in Fig. 1, showing the normalised $D_{IJ,ij}$ as a function of the pulsar-pulsar angular separation. On the left, we display the total variance computed over all pulsar pairs and frequencies. On the right, we perform the pulsar average by binning according to angular distances and averaging within each bin. The mean closely follows the HD prediction, while the scatter agrees with the expected cosmic variance Allen:2022dzg . The latter is invariant under the increase of the number of pulsars, reflecting the irreducible variance arising from the GWB realization in our Universe.

**Inference. ** In this work, we focus on the strong-signal regime and neglect the role of instrumental and intrinsic noise. While technically straightforward, incorporating noise in the analysis would primarily reduce the effective signal-to-noise ratio without qualitatively altering our conclusions (see Ref. Babak:2024yhu for an in-depth discussion). In practice, we neglect the noise contribution, $C_{n,IJ,ij}$, to the total cross-correlation introduced in Babak:2024yhu , and the full ensemble average covariance matrix for the timing residuals from pulsars $I,J$ and frequencies $i,j$ takes the compact form (21). This covariance matrix extends the one considered in Ref. Babak:2024yhu as it accounts for frequency-frequency correlations and an arbitrary GWB angular dependence, i.e., it does not rely on the assumption of an isotropic GWB. We stress again that we rely on the factorization in (15). Under this assumption, since frequency and pulsar indices do not mix, we can neglect the transmission function from the analysis without affecting the conclusions regarding anisotropies.

Assuming the data $\tilde{\delta t}_{I,i}$ to be Gaussian variables with zero mean and described only by their variance $C_{IJ,ij}(\bm{\theta})$, the log-likelihood can be written as Moran1951HypothesisTI ; Contaldi:2020rht ; Bond:1998zw ; Babak:2024yhu

[TABLE]

where we have ignored additive terms that are independent of model parameters $\bm{\theta}$, which includes only the coefficients $c_{\ell m}$ parametrizing GWB anisotropies. We do not infer the GWB spectral properties, as, thanks to the factorizability (15), they are irrelevant for analyses of anisotropies in the strong-signal regime. To speed up the inference, we adopt an iterative procedure (e.g., Franciolini:2025leq ). We find the maximum likelihood (ML) estimators for $\hat{\bm{\theta}}\equiv\hat{c}_{\ell m}$ in (6) by iteratively performing a gradient descent with a Fisher Information Matrix (FIM)-informed step size, i.e. $\hat{\bm{\theta}}_{\alpha}^{(i+1)}=\hat{\bm{\theta}}_{\alpha}^{(i)}+\delta\hat{\bm{\theta}}_{\alpha}^{(i)}$ with

[TABLE]

where the partial derivative of the log-likelihood reads

[TABLE]

with matrix multiplication and trace over on both pulsar $IJ$ and frequency $ij$ indeces (omitted to avoid clutter). These steps are repeated until the changes $|\delta\hat{\bm{\theta}}_{\alpha}^{(i)}|$ become smaller than a predefined convergence threshold, which we set to be $1\%$ of the current parameter uncertainty estimated from the FIM defined as (e.g., Bond:1998zw )

[TABLE]

Assuming the posterior distribution is nearly Gaussian, one can obtain the measurement uncertainties, as well as the model parameter correlations, by inverting the FIM (see, e.g., Coe:2009xf ). We validate the minimization procedure by also performing likelihood minimization with other numerical routines, confirming our findings. Also, we validated the FIM procedure by running an MCMC to confirm both ML and FIM estimates presented below. Finally, the accuracy of the inference is demonstrated by the P-P plots we report in App. B.

**Results. ** We analyse two distinct scenarios: isotropic and anisotropic GWB simulations, with various $N_{p}$. We will focus on the widely adopted angular power spectra

[TABLE]

Our results are $N_{\rm pix}$ (or $\Delta\hat{\Omega}$)-independent, as we adopt values well below the angular resolution of our PTA configurations. We use $S_{h}=f^{-7/3}$ as predicted by GW driven SMBHB inspirals Phinney:2001di and $N_{f}=15$.

Isotropic simulations. For isotropic GWB realizations, we inject signals with $c_{00}=1/\sqrt{4\pi}$ and all higher-order multipole coefficients set to zero. Fig. 2 (leftmost column) shows the distribution of reconstructed multipole power spectra normalized by their associated uncertainties $\sigma_{C_{\ell}}$, for various PTA configurations.

In the top panel, we see that with larger $N_{p}$, the monopole $C_{0}=1/4\pi$ is measured with higher and higher accuracy, following the expected scaling $\sigma_{C_{0}}\sim 1/\sqrt{N_{p}}$ Babak:2024yhu as we include autocorrelations in the analysis. Moving to higher order multipoles, we see that the inferred $\hat{C}_{\ell}$ are always positive, as in each realization, due to the stochastic nature of the signal, the inferred $\hat{c}_{\ell m}$ are non-zero. In other words, the distribution of recovered $\hat{C}_{\ell}$ exhibits a characteristic positive expectation value inherent to power spectrum estimators: although the underlying $c_{\ell m}$ coefficients have zero mean, they possess non-zero variance, and their quadratic combination $C_{\ell}$, being the sum of squares, yields strictly positive values. However, the key finding is that the ratio $\hat{C}_{\ell}/\sigma_{C_{\ell}}$ consistently peaks around ${\cal O}(1)$ across all higher multipoles, confirming that the induced and reconstructed anisotropies remain small compared to the measurement uncertainties, and one infers the signal to be statistically consistent with isotropy. This analysis demonstrates that cosmic variance does not induce spurious anisotropy detections.

Our conclusion remains valid across different simulation parameters. The recovery of isotropy is independent of the sky pixelization scheme, provided the strong-signal approximation remains valid. Moreover, changing the number of pixels does not affect the results since it only redistributes the same total variance (see footnote 3). Similarly, the result remains unaffected by varying the number of pulsars $N_{p}$ from 70 to 500. This shows that, while having larger arrays reduces measurement uncertainties (as established in previous work Depta:2024ykq ), a larger $N_{p}$ also decreases the cosmic variance because one takes more independent measurements of the sky. All in all, we find that isotropic backgrounds are correctly identified as such, even when accounting for cosmic variance in the simulation of single data realizations.

We expect these results to extend to more realistic observational scenarios. Including instrumental noise would only increase the measurement uncertainties, making the results even more consistent with isotropy. For sparse PTA networks where pulsar-limited rather than cosmic variance dominates the error budget, the statistical inference remains equally reliable.

Anisotropic simulations. To validate the ability of our framework to detect genuine anisotropies, we perform injections with non-zero multipole coefficients up to $\ell_{\text{max}}=2$, i.e., a maximally dipolar signal $P(\hat{\Omega})=(1+\cos\theta)/4\pi$, corresponding to an injection with only non-zero $\{c_{00}=1/(2\sqrt{\pi}),c_{10}=1/(2\sqrt{3\pi})\}$, and a maximally quadrupolar signal $P(\hat{\Omega})=3\cos^{2}\theta/4\pi$, corresponding to $\{c_{00}=1/(2\sqrt{\pi}),c_{20}=1/(\sqrt{5\pi})\}$. The corresponding results are shown in the center and right columns of Fig. 2.

A few comments are necessary. First, we observe that the monopole is recovered with the same significance as in the isotropic scenario. This is because different multipoles are independent, and the total (integrated) power remains unchanged under anisotropic injections, i.e., $\int{\rm d}^{2}\hat{\Omega}P(\hat{\Omega})=1$. Therefore, in the strong-signal regime, the uncertainty is driven solely by $N_{f}$ and $N_{p}$. Furthermore, for the multipoles with no injected signal, the distribution of $\hat{C}_{\ell}/\sigma_{C_{\ell}}$ remains of order ${\cal O}(1)$, i.e., the multipole is consistent with zero, as in the isotropic case. Finally, for the multipoles corresponding to the injection (dipolar: $\ell=1$, quadrupolar: $\ell=2$), the $C_{\ell}$ are measured with increasing significance as $N_{p}$ grows.

Conclusions. In this work, we have performed simulations of realistic PTA datasets by injecting GWB signals with different angular power spectra. Adopting an inference framework based on the Whittle likelihood (23), we have analyzed the case of signal-dominated observations while varying the number of pulsars. Our results confirm the fundamental limitations of PTA experiments to access the underlying $C_{\ell}$ distribution. Nevertheless, the inferred parameters remain robustly consistent with the assumptions underlying each injection, and the behaviour of anisotropies due to cosmic variance follows the expected trends. Our simulations confirm that cosmic variance cannot induce significant false anisotropic signals, provided the correct likelihood is adopted444As shown in Konstandin:2024fyo , modeling the time-domain cross-correlations $C_{IJ}$ as Gaussian variables, rather than the square of Gaussian variables, neglects higher order correlations affecting the analysis. . Even when individual realizations of an isotropic GWB exhibit apparent directional variations due to cosmic variance, the statistical uncertainties properly account for this effect. For an isotropic GWB, the recovered confidence intervals remain consistent with the null hypothesis of isotropy, providing a robust foundation for distinguishing astrophysical from cosmological GWB origins based on the search for anisotropies.

The cosmic variance in PTAs is analogous to the well-known CMB case after taking into account the more complicated structure of the response function arising from the three-dimensional network geometry. In CMB temperature measurements, cosmic variance at multipole $\ell$ is fundamentally limited by the $\sim 2\ell+1$ independent modes available on the celestial sphere. This results in large error bars, often shown as the cosmic variance ‘band’. For the CMB, all information we have is a 2D map of the sky, so assuming there are no limitations due to the instrument resolution, the number of measurements we have at multipole $\ell$ is simply $\sim 2\ell+1$. However, PTAs probe the GWB through a network of pulsar pairs, each acting as an independent detector with its unique directional sensitivity encoded in the generalized overlap reduction functions $\Gamma_{IJ,\ell m}$. This geometric advantage allows PTAs to potentially overcome some cosmic variance limitations by increasing the number of pulsars. Each additional pulsar contributes $N_{p}-1$ new correlations, providing additional independent measurements of the same underlying multipole structure. Therefore, we confirm that future PTA datasets could detect anisotropies as small as $C_{\ell=1}=9\times 10^{-3}(1/N_{f})(100/N_{p})$ at $3\sigma$ C.L., where $N_{f}$ indicates the number of signal dominated observed frequencies, as reported in Depta:2024ykq .

Appendix A Maximum Number of Multipoles Constrained with $N_{p}$ Pulsars

We estimate the number of independent multipole moments $c_{\ell m}$ that can be constrained from GWB anisotropies using $N_{p}$ pulsars. The estimator for the multipole coefficients can be written as

[TABLE]

where we have factored out a frequency-dependent amplitude $c(f)$. A formal solution of this linear system involving $(\ell_{\rm max}+1)^{2}$ unknowns can be easily obtained. We rewrite this system of equations indicating as $\alpha$ and $\beta$ the index of pulsar pairs and $c_{\ell m}$ coefficients, to get

[TABLE]

The matrix inversion requires a square matrix of full rank, indicating the maximal number of ${c_{\ell m}}^{\prime}s$ which can be inferred is given by the length of $\{\alpha\}$, i.e. the number of pulsar pairs $N_{\rm pairs}$. Given $N_{p}$ pulsars, the number of independent cross-correlation pairs is $N_{\rm pairs}={N_{p}(N_{p}-1)}/{2}$. Assuming that only cross-correlations contribute significantly (neglecting autocorrelations relevant for low $\ell$, see Depta:2024ykq ), and considering a single frequency bin for simplicity, the estimation of $c_{\ell m}$ can be schematically written as

[TABLE]

We have verified explicitly that the rank of $\Gamma_{\alpha\beta}^{-1}$ is generically equal to the number of pulsar pairs, indicating that the number of multipole coefficients that can be independently constrained scales as

[TABLE]

This scaling differs from the $\ell_{\rm max}\sim\sqrt{N_{p}}$ estimate quoted in NANOGrav:2023tcn ; Romano:2016dpx . We provide a physical interpretation for this scaling argument in the main text.

Appendix B Model validation through probability - probability (P-P) plots

We validate our inference procedure, as well as the accuracy of the uncertainty estimates obtained using the FIM, by presenting the probability-probability (P-P) plot shown in Fig. 3.

The P–P plot compares the cumulative distribution function (CDF) of the inferred quantile for the injected $c_{\ell m}$ values in the posterior, with the uniform distribution expected under correct calibration, providing a diagnostic of the statistical consistency of the inferred posteriors. To construct it, we generate $1000$ independent data realizations. For each realization, we perform the full inference and determine the quantile at which the injected $c_{\ell m}$ value lies within the corresponding posterior distribution. The resulting set of quantiles is then sorted in ascending order and plotted against the theoretical expectation for ordered quantiles, which follows an equally spaced sequence between [math] and $1$. The $1\sigma$ sampling uncertainty on the P–P plot is computed assuming binomial statistics: for a given theoretical cumulative probability $p$, the standard deviation is

[TABLE]

where $N$ is the total number of realizations. This uncertainty band, shown as a gray shaded region around the diagonal, quantifies the expected statistical scatter for a finite sample of realizations. See App. A of Franciolini:2025leq for more explanations and a dictionary on how to interpret P-P plots.

In Fig. 3, from left to right, we display the P-P plot for each $c_{\ell m}$ coefficient up to $\ell_{\text{max}}=3$, corresponding respectively to injections with: monopole only, monopole plus maximal dipole, and monopole plus maximal quadrupole. All results are based on $10^{3}$ realizations, each using $N_{p}=120$ and $N_{p}=500$ pulsars. We find excellent agreement with the theoretical expectation in all cases, confirming both the validity of the inference method and the reliability of the associated uncertainty estimates.

Appendix C Distribution of $\hat{C}_{\ell}$

In Fig. 4 we show the distribution of the ML estimators $\hat{C}_{\ell}$ obtained performing 1000 realizations. We assume the same parameters as in Fig. 2. Each color (row) corresponds to a different number of pulsars, and each column corresponds to a different multiple index $\ell$. The simulations assume an isotropic GWB. We remove the mean from the $\ell=0$ result. Superimposed to the histograms, we show in black dashed line the expected distribution of $C_{\ell}$, built assuming a Gaussian $d_{I,i}$ for an isotropic signal. This reads Franciolini:2025leq

[TABLE]

where $D=\sum_{\ell}|\hat{c}_{\ell m}|^{2}/(2\ell+1)$, the variance is $P=\text{mean}(\hat{C}_{\ell})$ across all the realizations, and $N_{s}=(2\ell+1)/2$ is the number of effective degrees of freedom in the sum at each $\ell$. This distribution exactly matches a $\chi^{2}$ distribution with the appropriate degrees of freedom $k=2\ell+1$, as expected for a variable built as the sum of $2\ell+1$ Gaussian variables.