# Measuring anisotropies in the PTA band with cross-correlations

Giulia Cusin, Cyril Pitrou, Martin Pijnenburg, Alberto Sesana

► To cite this version:

Giulia Cusin, Cyril Pitrou, Martin Pijnenburg, Alberto Sesana. Measuring anisotropies in the PTA band with cross-correlations. Monthly Notices of the Royal Astronomical Society, 2025, 2884, pp.2896. 10.1093/mnras/staf1074. hal-04986176

HAL Id: hal-04986176

https://hal.science/hal-04986176v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Measuring anisotropies in the PTA band with cross-correlations

Giulia Cusin,<fnref n="1" />,<fnref n="2" />* Cyril Pitrou,<fnref n="1" /> Martin Pijnenburg<fnref n="2" /> and Alberto Sesana<fnref n="3" />,<fnref n="4" />,<fnref n="5" />

Accepted 2025 June 27. Received 2025 June 5; in original form 2025 April 15

## ABSTRACT

The astrophysical gravitational wave background in the nanohertz (nHz) band is expected to be primarily composed of the superposition of signals from binaries of supermassive black holes. The spatial discreteness of these sources introduces shot noise, which, in certain regimes, would overwhelm efforts to measure the anisotropy of the gravitational wave background. Moreover, the fact the time-residual map has a finite angular resolution and the presence of pulsar noise, affects our ability to construct the angular power spectrum of the anisotropy from a time-residual map (finite resolution noise). In this work, we explicitly demonstrate, starting from first principles, that cross-correlating a gravitational wave background map with a sufficiently dense galaxy survey can mitigate this issue. This approach could potentially reveal underlying properties of the gravitational wave background that would otherwise remain obscured. We quantify both the shot noise and the finite resolution noise level and show that cross-correlating the gravitational wave background with a galaxy catalogue improves by more than one order of magnitude the prospects for a first detection of the background anisotropy by a gravitational wave observatory operating in the nHz frequency range. In particular, we find that with a futuristic scenario with an effective number of frequencies equal to $N_f = 10$, the detection of the spectral amplitude can be achieved combining the first 20 multipoles, with a threshold to resolve single events that corresponds to the Square Kilometer Array (SKA) sensitivity. Increasing observation time, pulsar number or reducing the pulsar white noise considerably improves the detection significance.

**Key words:** gravitational waves – large-scale structure of Universe.

## 1 INTRODUCTION

Pulsar timing arrays (PTAs) were used to provide the first evidence of a stochastic gravitational wave (GW) background in the nHz band (Antoniadis et al. 2023; Reardon et al. 2023; Agazie et al. 2023a; Xu et al. 2023), with an integrated energy density of $\Omega_{GW} = 9.3_{-4.0}^{+5.8} \times 10^{-9}$ (Agazie et al. 2023a) and evidence of spatial correlation among different pulsar redshifts following the Hellings–Downs (HD) function. Various detection methods are reviewed in Romano & Cornish (2017), and the physical interpretation of the HD correlation is discussed in Jenet & Romano (2015), Romano & Allen (2024), Wu, Bi & Huang (2024), and Grimm et al. (2025a,b). Methods for mapping the background have also been developed (Anholm et al. 2009; Mingarelli et al. 2013; Semenzato et al. 2024; Depta et al. 2025; see also Cusin, Durrer & Ferreira 2019a; Bernardo, Liu & Ng 2024, for a discussion of polarization with Stokes parameters), and it has recently been proposed that a loud source may account for the evidence of anisotropy reported in Grunthal et al. (2024). The most recent analysis of the HD signature by NANOGrav (Agazie et al. 2025) is expressed in multipole space (see Gair et al.

2014; Roebber & Holder 2017; Hotinli, Kamionkowski & Jaffe 2019; Qin et al. 2019; Bernardo & Ng 2022, 2023, 2025; Allen 2024; Nay et al. 2024; Pitrou & Cusin 2025 for a harmonic formulation of the HD correlation), with a clear detection of the quadrupole and only marginal evidence for the octupole. The significance of background detection is expected to improve in the future, as the signal-to-noise ratio (SNR) of PTA observables increases with observation time (Siemens et al. 2013; Pol, Taylor & Romano 2022; Nay et al. 2024).

The origin of the signal is still uncertain, with anisotropies in energy density providing a key observable to pin down its nature. In fact, while the overall GW signal generated by an astrophysical population of sources will inevitably show some degree of anisotropy, a cosmological background sourced in the early Universe is expected to be isotropically distributed (at least at angular scales of tens of degrees, probed by PTAs).¹Current sky bounds are set by NANOGrav (Agazie et al. 2023c), while Depta et al. (2025) uses Fisher Information Matrix forecasts to show that sensitivity to angular power spectra

<footnote n="1">¹Institut d'Astrophysique de Paris, UMR-7095 du CNRS et de Sorbonne Université, 75014 Paris, France</footnote>

<footnote n="2">²Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, Quai E. Ansermet 24, CH-1211 Genève 4, Switzerland</footnote>

<footnote n="3">³Dipartimento di Fisica 'G. Occhiali', Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="4">⁴INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy</footnote>

<footnote n="5">⁵INAF – Osservatorio Astronomico di Brera, via Brera 20, I-20121 Milano, Italy</footnote>

<footnote n="1">¹Cosmological backgrounds can in fact exhibit large kinematic anisotropies in the presence of spectral features (Cusin & Tasinato 2022; Cruz et al. 2024b, see also Cruz et al. 2024a) for detection prospects exploiting synergies with astrometry.</footnote>

<footnote>* E-mail: giulia.cusin@iap.fr</footnote>
improves rapidly with the number of pulsars monitored, expected to
increase by 20–40 per year.

The astrophysical gravitational wave background (AGWB) in the frequency band targeted by PTA experiments is expected to be primarily dominated by signals from binary systems of massive black holes (BBHs) in the inspiraling phase (Rajagopal & Romani 1995; Jaffe & Backer 2003; Sesana, Vecchio & Colacino 2008). Although the duration of these signals is exceedingly long compared to typical observation times, making the resulting background effectively continuous and stationary over the observation period (but see Falxa, Leclere & Sesana 2025), the sources themselves have a discrete spatial distribution. The number of sources fluctuates according to a Poisson distribution, introducing shot noise into the angular power spectrum. Additionally, sources are located in galaxies, hence they are expected to be a biased tracer of the large-scale structure distribution. In certain situations, shot noise can obscure or dominate the clustering component, making it challenging to isolate the underlying anisotropic structure of the AGWB spectrum. For example, Semenzato et al. (2024) recently showed using numerical simulations that the autocorrelation map is Poisson-noise dominated, for the black hole models used in the analysis.

Cross-correlating a shot noise dominated background map with a dense galaxy map has been proved to be an effective method to alleviate shot noise limitation in the context of Earth-based detectors, where the Poissonian nature of sources is both spatial and temporal, due to the very short emission duration of mergers in band (Jenkins & Sakellariadou 2019; Jenkins, Romano & Sakellariadou 2019; Cusin et al. 2019b; Alonso et al. 2020, 2024; Yang et al. 2023). In particular, in Alonso et al. (2020) it is shown, in the context of ground-based detectors and considering an ideal scenario where shot noise is the only noise component, that the SNR of the cross-correlation outperforms that of the autocorrelation by several orders of magnitude. See also Tian, Ding & Kou (2024) for a general analysis of cross-correlation methods applicable to PTA.

In the context of PTAs, in addition to shot noise, there exists an additional noise component stemming from the finite angular resolution of the time-residual map. This resolution is determined by the number of pulsars in the network, characterized by $L_{max} \sim \sqrt{N_p}$. Even in a perfect experiment with zero pulsar noise, this limitation impacts our ability to constrain the anisotropy of the time-residual map, as discussed in Hotinli et al. (2019). Furthermore, it propagates into the reconstruction of the angular power spectrum of the background energy density, which is the primary observable used to quantify deviations from isotropy, see Pitrou & Cusin (2025). We will refer to this additional noise component as finite resolution noise. The size of the finite resolution noise is set by the number of effective frequencies which for a perfect experiment (where pulsars are noiseless) is extremely large ($N_f = 2Tf_{max} \gg 1$) while it is much suppressed in the presence of pulsar noise (Pitrou & Cusin 2025). Our aim here is to explore the issues of shot noise and finite resolution noise in the context of PTA observations, and to show that cross-correlating a background map with a dense galaxy catalogue can mitigate both noise contributions.

After having identified the observables we are interested in Section 2.1, we present a first principle derivation of the shot noise contribution of auto and cross-correlation maps for an AGWB in the nHz band. We need to account for two distinct layers of stochasticity describing the distribution of massive black hole binaries: a given galaxy may or not contain a massive black hole binary, and galaxies are themselves a discrete and random biased sampling of the underlying matter density. Using properties of compound statistics, we show that while the shot noise of the autocorrelation map is

proportional to the inverse number of sources (black hole binaries), the cross-correlation shot noise is proportional to the inverse galaxy number, hence it is much suppressed. In Section 2.3, we turn into the computation of the finite resolution noise, for both auto and cross-correlation angular spectra, adapting results from Hotinli et al. (2019) and Pitrou & Cusin (2025). In Section 3, we present numerical results: we consider the theoretical framework of Cusin, Pitrou & Uzan (2017, 2018) and Pitrou, Cusin & Uzan (2020) to describe anisotropies and we consider the astrophysical models of Sesana (2013) and Rosado, Sesana & Gair (2015) to describe the underlying massive black hole population. The AGWB strain map is dominated by a few very bright sources which obscure a background from an unresolved population, whose anisotropy map is a representative sample of large-scale structure anisotropies. We therefore consider the effect of filtering out bright sources, for three different threshold values of source strain: a threshold corresponding to the European Pulsar Timing Array (EPTA) strain curve ('no cut' situation as effectively there are no resolvable sources resulting from the filtering procedure), a threshold corresponding to the projected SKA strain sensitivity, and an intermediate value. For each situation we compute the auto and cross-correlation spectrum with galaxy distribution, and the corresponding shot noise and finite resolution noise. We estimate the SNR of auto and cross-correlation maps and we show that the detectability depends on two factors: on the threshold for filtering out resolvable events, which depends on the sensitivity of the instrument, and on the number of effective frequencies one observatory has access to $N_f$, which depends on the number of pulsars monitored, on the time of observation and on the pulsar noise level (both white noise and red noise), see section VII C 2 of Pitrou & Cusin (2025). In that reference, it was actually found that the number of effective frequency bins for fixed observational setting, is a decreasing function of the multipole $L$ of the time residual map. Here, we consider an effective flat number $N_f \sim 10$, which has to be understood as the average $N_f$ over the number of multipoles one has access to with a given map resolution. The number $N_f = 10$ corresponds to a futuristic scenario,<fnref n="2" /> which can only be achieved by considerably increasing observation time and the cadence of pulsar monitoring. We also consider the case of a perfect experiment with $N_f = 2Tf_{max} \gg 1$. In this latter case, the finite resolution noise is very small (pulsar white noise is negligible) and shot noise dominates. For $N_f = 10$, we find that detection can be achieved with cross-correlation, integrating the background signal over frequency and combining the first multipoles up to $\ell = 20$ with a threshold for resolvable events at $10^{-16}$. For $N_f \gg 1$ detectability is obtained combining the first five multipoles for threshold $10^{-16}$ and the first 13 multipoles for $10^{-15}$. In this article, technical derivations are presented in a series of appendices.

**Notation.** The index $L$ denotes the angular decomposition of a time-residual map. The angular resolution of the map is related to the number of pulsars in the network via $L_{max} \sim \sqrt{N_p}$. The lower case $\ell$ is used for the multipole decomposition of the energy density correlation function, and $\ell_{max} = 2L_{max}$. We call perfect experiment an observational setting where the pulsar noise is negligible. Hence a perfect experiment still is characterized by a finite pulsar number, observation time, and cadence of observation. In our analysis, the entire redshift distribution of galaxies is utilized, employing optimal Wiener weights, hence assuming that the redshift distribution of GW sources is known.

<footnote n="2">²Assuming the pulsar noise characteristics of Agazie et al. (2023b), getting $N_f(L = 5) \sim 1$ would require 450 pulsars monitored during 25 yr.</footnote>
## 2 GENERAL CONCEPTS AND DERIVATIONS

### 2.1 Our observables

The observed GW energy density parameter, $\Omega_{GW}$ is defined as the background energy density $\rho_{GW}$ per units of logarithmic frequency $f$ and solid angle $n$, normalized by the critical density of the Universe today $\rho_c$. It can be divided into an isotropic background contribution $\bar{\Omega}_{GW}$ and a contribution from anisotropic perturbations $\delta\Omega_{GW}$ (Cusin et al. 2017, 2018, 2019b):

$$ \Omega_{GW}(\mathbf{n}, f) = \frac{f}{\rho_c} \frac{d^3 \rho_{GW}}{d^2 \mathbf{n} \, d f}(\mathbf{n}, f) = \frac{\bar{\Omega}_{GW}(f)}{4\pi} + \delta\Omega_{GW}(\mathbf{n}, f), \quad (1) $$

where the isotropic background spectrum can be written as the integral over conformal distance $r$ (where we set the speed of light equal to unity, $c = 1$):

$$ \bar{\Omega}_{GW}(f) = \int dr \, \partial_r \bar{\Omega}_{GW}(f, r), \quad (2a) $$

where the integrand plays the role of an astrophysical kernel that contains information on the local production of GWs at galaxy scales. This kernel can be built out of a model of black holes formation, or a catalogue implementing such model, as detailed in Appendix D. Different astrophysical models give quite different predictions for this kernel, see e.g. Cusin et al. (2019b) for an explorative approach in the frequency band of Earth-based and space-based detectors.

We now consider observables integrated over the whole frequency range a given observatory is sensitive to. Explicitly, we introduce

$$ \bar{\Omega}_{GW} = \int d \log f \, \Omega_{GW}(f), \quad (3a) $$

$$ \Omega_{GW}(\mathbf{n}) = \int d \log f \, \Omega_{GW}(\mathbf{n}, f), \quad (3b) $$

$$ \partial_r \bar{\Omega}_{GW}(r) = \int d \log f \, \partial_r \bar{\Omega}_{GW}(f, r), \quad (3c) $$

$$ \delta\Omega_{GW}(\mathbf{n}) = \int d \log f \, \delta\Omega_{GW}(\mathbf{n}, f). \quad (3d) $$

We observe that since $\delta\Omega_{GW}$ is a stochastic quantity, it can correlate with other cosmological stochastic observables. An interesting observable to look at is the cross-correlation of the AGWB with the distribution of galaxies, i.e. with the galaxy number counts $\Delta$ defined as the overdensity of the number of galaxies per unit of redshift and solid angle

$$ \Delta(\mathbf{n}, z) = \frac{N(z, \mathbf{n}) - \bar{N}(z)}{\bar{N}(z)}. \quad (4) $$

First, if astrophysical GW sources are located in galaxies, we would expect the SGWB and the galaxy distribution to have a high correlation level. Second, cross-correlating with galaxies helps to mitigate the problem of shot noise and to possibly extract the clustering information out of the shot noise threshold. This has been shown in the context of ground-based detectors, where shot noise is dominated by its popcorn component due to the transient emission of sources, see Jenkins & Sakellariadou (2019), Jenkins et al. (2019), Cusin et al. (2019b), and Alonso et al. (2020).

We want to show that this is also the case when the emission is continuous in time, like in the PTA band. Finally, by cross-correlating with the galaxy distribution at different redshifts, one could try to get a tomographic reconstruction of the redshift distribution of sources. In this work, our goal is to maximize the detection chances hence we do not follow this approach.

The angular power spectrum of the GW and galaxy counts autocorrelations and for their cross-correlations are defined as

$$ (2\ell + 1) C_{\ell}^{GW} \equiv \sum_{m=-\ell}^{\ell} \langle a_{\ell m} a_{\ell m}^* \rangle, \quad (5a) $$

$$ (2\ell + 1) C_{\ell}^{\Delta} \equiv \sum_{m=-\ell}^{\ell} \langle b_{\ell m} b_{\ell m}^* \rangle, \quad (5b) $$

$$ (2\ell + 1) C_{\ell}^{GW\Delta} \equiv \sum_{m=-\ell}^{\ell} \langle a_{\ell m} b_{\ell m}^* \rangle, \quad (5c) $$

where the bracket denotes an ensemble average and $a_{\ell m}$ and $b_{\ell m}$ are the coefficients of the spherical harmonics decomposition of the AGWB energy density and galaxy number counts, respectively. Explicitly

$$ \delta\Omega_{GW}(\mathbf{n}) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} a_{\ell m} Y_{\ell m}(\mathbf{n}), \quad (6a) $$

$$ \Delta(\mathbf{n}) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} b_{\ell m} Y_{\ell m}(\mathbf{n}). \quad (6b) $$

It can be shown that the angular power spectra of the auto and cross-correlation are given by (Cusin et al. 2017):

$$ C_{\ell}^{GW} = \frac{2}{\pi} \int dk \, k^2 |\delta\Omega_{GW,\ell}(k)|^2, \quad (7a) $$

$$ C_{\ell}^{\Delta} = \frac{2}{\pi} \int dk \, k^2 |\Delta_{\ell}(k)|^2, \quad (7b) $$

$$ C_{\ell}^{GW\Delta} = \frac{2}{\pi} \int dk \, k^2 \delta\Omega_{GW,\ell}^*(k) \Delta_{\ell}(k), \quad (7c) $$

where $k$ is the Fourier mode norm. Keeping only the leading-order contribution to the anisotropy given by clustering, which is expressed by equation (D10), we have

$$ \delta\Omega_{GW,\ell}(k) = \int d \log f \, \delta\Omega_{GW,\ell}(k, f), \quad (8) $$

$$ \delta\Omega_{GW,\ell}(k, f) = \int \frac{dr}{4\pi} \partial_r \bar{\Omega}_{GW}(f, r) [b(r) \delta_{m,k}(r) j_\ell(kr)], \quad (9) $$

where $j_\ell$ are spherical Bessel functions, while $\delta_m$ is the dark matter overdensity, related to galaxy overdensity via the bias factor $b(r)$. The corresponding contribution from galaxy overdensities (retaining only the dominant clustering contribution) reads

$$ \Delta_{\ell}(k) = \int dr \, W(r) [b(r) \delta_{m,k}(r) j_{\ell}(kr)], \quad (10) $$

where $W(r)$ is a window function normalized to one which selects the redshift bin in the galaxy catalogue we want to consider in the cross-correlations.

The framework presented so far is general and can be applied to any astrophysical background component across any frequency band. We will now focus on the specific case of PTAs, where the sources are BBHs residing in galaxies. In the next section, we will estimate the contribution of shot noise in this context, both for the autocorrelation spectrum and for the cross-correlation with the galaxy distribution.

Recently, the variance on the background energy density due to the discrete distribution of sources has been studied in Lamb & Taylor (2024), while Sato-Polito & Kamionkowski (2024) looks at anisotropies due to discrete source distribution, for different astrophysical models.

### 2.2 Shot noise

Our goal is to demonstrate from first principles that if the source distribution is a sub-sample of the underlying galaxy distribution, cross-correlations are less affected by shot noise than autocorrelations. We will also derive explicit expressions for both contributions.
## 2.2.1 Compound statistics

When counting the number of massive BBHs, we have two levels of stochasticity: one galaxy may or may not contain a black hole binary, and the distribution of galaxies is stochastic and modeled as Poisson-distributed. We write the number of black hole binaries as

$$N_{\text{BBH}} = \sum_{i=1}^{N_G} X_i, \qquad (11)$$

where $X_i$ follows a binomial distribution $X_i \sim B(1, \beta_X)$, hence $\langle X_i \rangle = \beta_X$ and $\text{Var}(X_i) = \beta_X(1 - \beta_X)$. In practice $\beta_X \ll 1$, so that most galaxies do not contain a binary emitting in the relevant frequency band.

Using properties of compound statistics derived in Appendix A and recalling that the covariance of two variables Y and Z is defined as $\text{Var}(Y, Z) = \langle YZ \rangle - \langle Y \rangle \langle Z \rangle$, one finds

$$\langle N_{\text{BBH}} \rangle = \beta_X \langle N_G \rangle, \qquad (12a)$$

$$\text{Var}(N_{\text{BBH}}) = \beta_X \langle N_G \rangle, \qquad (12b)$$

$$\text{Var}(N_{\text{BBH}}, N_G) = \beta_X \langle N_G \rangle. \qquad (12c)$$

We observe that $N_{\text{BBH}}$ is also a Poissonian variable.

## 2.2.2 A heuristic derivation

Let us try to have an intuition of why cross-correlating the stochastic background map to a galaxy map helps with the shot noise problem. To present a heuristic derivation, we simplify our treatment and notation: we fix a frequency bin, we consider a set of radial bins and we assume that all sources have the same luminosity L. Recalling that only objects in the same distance bin correlate, we can schematically write

$$\Omega_{\text{GW}} = \sum_{r} \Omega_{\text{GW}}(r) = \sum_{r} N_{\text{BBH}}(r) \frac{L}{\rho_c 4\pi r^2}, \qquad (13)$$

where $L/(4\pi r^2)$ is the energy flux associated to the GW event. For the galaxy distribution, we write

$$N_G = \sum_r N_G(r), \qquad (14a)$$

$$\Delta = \sum_r \frac{N_G(r) - \langle N_G \rangle}{\langle N_G \rangle}. \qquad (14b)$$

Then using results from previous section, one has

$$\text{Var}(\Delta) = \frac{1}{\langle N_G \rangle^2} \sum_r \text{Var}(N_G(r)) = \frac{1}{\langle N_G \rangle}, \qquad (15a)$$

$$\begin{aligned} \text{Var}(\Omega_{\text{GW}}) &= \sum_r \left( \frac{L}{\rho_c 4\pi r^2} \right)^2 \text{Var}(N_{\text{BBH}}(r)) \\ &= \sum_r \frac{\langle \Omega_{\text{GW}}(r) \rangle^2}{\langle N_{\text{BBH}}(r) \rangle}, \end{aligned} \qquad (15b)$$

$$\begin{aligned} \text{Var}(\Omega_{\text{GW}}, \Delta) &= \frac{1}{\langle N_G \rangle} \sum_r \frac{L}{\rho_c 4\pi r^2} \text{Var}(N_{\text{BBH}}(r), N_G(r)) \\ &= \frac{\langle \Omega_{\text{GW}} \rangle}{\langle N_G \rangle}. \end{aligned} \qquad (15c)$$

Note that $\text{Var}(\Omega_{\text{GW}})$ is suppressed by the number of massive black hole binaries (at fixed $\Omega_{\text{GW}}$). Since $\langle N_G \rangle \gg \langle N_{\text{BBH}} \rangle$ (for $\beta_X \ll 1$), the relative cross-correlation variance is expected to be smaller than the autocorrelation.

## 2.2.3 A more formal derivation of shot noise

Let us refine the derivation above, giving up the assumption that sources have all the same luminosity. We consider a binned catalogue of sources (in our case massive black holes in binaries) and look at their correlation function. Define $\Delta_{i,f,M}$ to be the source overdensity in a spatial bin $i$, frequency bin $f$, and chirp mass bin $M$. The two-point function of this observable $\langle \Delta_{i,f,M} \Delta_{j,f',M'} \rangle$ is given by two terms

$$\langle \Delta_{i,f,M} \Delta_{j,f',M'} \rangle = \frac{1}{V_c \Delta \log f \Delta \log M a^3 \bar{n}_s} \delta_{ff'} \delta_{MM'} \delta_{ij} + \xi_{ij}^s, \qquad (16)$$

where $\xi_{ij}^s$ is the binned correlation function of the underlying, smooth, density field, $\delta_s$, and $V_c \Delta \log f \Delta \log M a^3 \bar{n}_s$ is the mean number of sources in a spatial pixel of comoving volume $V_c$, frequency pixel of volume $\Delta \log f$, and chirp mass pixel $\Delta \log M$. Finally, $\delta_{ij}$ is the Kronecker symbol which arises from the Poisson sampling of the underlying, smooth density field, while the Kronecker in frequency and chirp mass space comes from the fact that only sources emitting at the same frequency are correlated. If we take the continuous limit of a binned survey, we find

$$\begin{aligned} \Delta_{i,f,M} &\rightarrow \delta_s(\mathbf{x}, \log f, \log M), \\ \delta_{ij} &\rightarrow V_c \delta^{(3)}(\mathbf{x} - \mathbf{x}'), \\ \delta_{ff'} &\rightarrow \Delta \log f \delta(\log f - \log f'), \\ \delta_{MM'} &\rightarrow \Delta \log M \delta(\log M - \log M'). \end{aligned} \qquad (17)$$

It is convenient to split the observed source overdensity into a smooth clustering part $\delta_s(\mathbf{x}, \log f, \log M)$ and Poissonian contribution

$$\hat{\delta}_s = \delta_s(\mathbf{x}, \log f, \log M) + \delta_s^{\text{shot}}(\mathbf{x}, \log f, \log M), \qquad (18)$$

where with a hat we denote observable quantities and the arguments on the left-hand side are understood for shortness. The first contribution on the right-hand side is due to clustering while for the Poissonian contribution one has

$$\begin{aligned} \langle \delta_s^{\text{shot}}(\mathbf{x}, \log f, \log M) \delta_s^{\text{shot}}(\mathbf{x}', \log f', \log M') \rangle \\ = \frac{1}{a^3 \bar{n}_s} \delta^{(3)}(\mathbf{x} - \mathbf{x}') \delta(\log f - \log f') \delta(\log M - \log M'), \end{aligned} \qquad (19)$$

where the quantity at the denominator $n_s = n_s(z, \log f, \log M)$ is the density of binaries per units $\log f$ and $\log M$. We use this result to compute the correlation function of the background anisotropies accounting for both clustering and Poissonian components. We assume in first approximation that the clustering of sources follows the clustering of galaxies i.e. $\delta_s(\mathbf{x}, \log f, \log M) = \delta_G(\mathbf{x})$.<fnref n="3" /> We get

$$\begin{aligned} \langle \delta \hat{\Omega}_{\text{GW}}(\mathbf{n}) \delta \hat{\Omega}_{\text{GW}}(\mathbf{n}') \rangle \\ = \frac{1}{(4\pi)^2} \int d\mathbf{r} \int d\log f \int d\log M \frac{\partial \bar{\Omega}_{\text{GW}}}{\partial \mathbf{r} \partial \log M} \\ \cdot \int d\mathbf{r}' \int d\log f' \int d\log M' \frac{\partial \bar{\Omega}_{\text{GW}}}{d\mathbf{r}' d\log M'} \\ \cdot \left[ \langle \delta_G(\mathbf{x} = \mathbf{n}\mathbf{r}) \delta_G(\mathbf{x}' = \mathbf{n}'\mathbf{r}') \rangle \right. \\ \left. + \delta^{(3)}(\mathbf{x} - \mathbf{x}') \delta(\log f - \log f') \delta(\log M - \log M') \frac{1}{a^3 \bar{n}_s} \right], \end{aligned} \qquad (20)$$

<footnote n="3">³Notice that this is just an approximation and a more refined model should e.g. consider that in fact massive galaxies are more likely to host black holes in binaries (Rosado & Sesana 2014). However this would not change the result for shot noise, whose computation is the main goal of our treatment.</footnote>
where we made use of equations (3d) and (D10). Expanding in terms of spherical harmonics

$$ \delta\hat{\Omega}_{\text{GW}}(\mathbf{n}) = \sum_{\ell m} \hat{a}_{\ell m} Y_{\ell m}(\mathbf{n}), \qquad (21) $$

and using

$$ \delta^{(3)}(\mathbf{x} - \mathbf{x}') = \frac{\delta(r - r')}{r^2} \sum_{\ell m} Y_{\ell m}(\mathbf{n}) Y_{\ell m}^*(\mathbf{n}'), \qquad (22) $$

we find that

$$ \begin{aligned} \langle \hat{a}_{\ell m} \hat{a}_{\ell' m'}^* \rangle = & \langle a_{\ell m} a_{\ell' m'}^* \rangle + \delta_{mm'} \delta_{\ell \ell'} \\ & \cdot \int \frac{dr}{(4\pi)^2} \int d\log f \int d\log \mathcal{M} \left| \frac{\partial \hat{\Omega}_{\text{GW}}}{\partial r \partial \log \mathcal{M}} \right|^2 \frac{1}{r^2} \frac{1}{a^3 \bar{n}_s}, \end{aligned} \qquad (23) $$

where the first contribution is the theoretical one, given in equation (5a) and the second the shot noise component. Using the standard definition

$$ \langle \hat{a}_{\ell m} \hat{a}_{\ell' m'}^* \rangle = \hat{C}_{\ell} \delta_{\ell \ell'} \delta_{mm'}, \qquad (24) $$

we immediately find that

$$ \hat{C}_{\ell} = C_{\ell}^{\text{GW}} + N_{\text{shot}}^{\text{GW}}, \qquad (25) $$

where

$$ N_{\text{shot}}^{\text{GW}} = \int \frac{dr}{(4\pi)^2} \int d\log f \int d\log \mathcal{M} \left| \frac{\partial \hat{\Omega}_{\text{GW}}}{\partial r \partial \log \mathcal{M}} \right|^2 \frac{1}{r^2} \frac{1}{a^3 \bar{n}_s}. \qquad (26) $$

One can repeat the same reasoning for the cross-correlation and the galaxy autocorrelation. Using results of the heuristic argument, we can directly see that

$$ N_{\text{shot}}^{\text{GW}\Delta} = \int \frac{dr}{4\pi r^2} \frac{1}{a^3 \bar{n}_G(r)} W(r) \partial_r \bar{\Omega}_{\text{GW}}(r), \qquad (27a) $$

$$ N_{\text{shot}}^{\Delta} = \int \frac{dr}{r^2} \frac{1}{a^3 \bar{n}_G(r)} W^2(r), \qquad (27b) $$

which scale as $1/(a^3\bar{n}_G)$ which is the inverse of the comoving number density of galaxies. In the Appendix E, we explain how the shot noise contribution for the GW autocorrelation can be estimated from a catalogue, and we prove the equivalence between the discrete and continuous descriptions.

To conclude: we found that while shot noise of the AGWB energy density autocorrelation scales as the inverse number of sources (black hole binaries), the shot noise of the cross-correlation with the galaxy distribution scales as the inverse number of galaxies, hence shot noise of a cross-correlation map is much more suppressed than shot noise of the AGWB autocorrelation.

## 2.3 Noise due to finite resolution of the time-residual map

Even for a perfect experiment (i.e. for zero pulsar noise), there is an additional noise component that affects our ability to reconstruct the angular power spectrum of the background intensity (and its cross-correlation with the large-scale structure), due to the limited angular resolution of a given PTA observatory. Indeed the angular power spectrum is built from a map of anisotropies, and the minimum anisotropy that can be detected is set by the map resolution $L_{\text{max}}$, in turn related to the pulsar number $L_{\text{max}} \sim \sqrt{N_p}$, see Hotinli et al. (2019) for a detailed derivation. This in turn propagates on our ability to construct the angular power spectrum. A detailed derivation of this

can be found in section VE of Pitrou & Cusin (2025), and summarized in Appendix B. The resulting noise is given by

$$ N_{\ell}^{\text{GW}}(L_{\text{max}}) = \sqrt{\frac{2}{2\ell + 1} \left(\frac{\bar{\Omega}_{\text{GW}}}{4\pi}\right) B_{\ell m}} \cdot \left[ \left(\frac{\bar{\Omega}_{\text{GW}}}{4\pi}\right)^2 B_{\ell m}^2 + 2C_{\ell}^{\text{GW}} \right]^{1/2}, \qquad (28) $$

where

$$ B_{\ell m}^{-2} = \frac{N_f}{8\pi} \sum_{L_1 L_2} \delta_{\ell L_1 + \ell L_2 + \ell} (2L_1 + 1)(2L_2 + 1) \begin{pmatrix} \ell & L_1 & L_2 \\ 0 & -2 & 2 \end{pmatrix}^2, \qquad (29) $$

where $\delta_{\ell L_1 + \ell L_2 + \ell} \equiv (1 + (-1)^{\ell_1 + \ell_2 + \ell})/2$ and $N_f$ is equal to the number of effective frequencies, which for a perfect experiment is given by $N_f = 2Tf_{\text{max}}$. In a realistic experiment, $N_f$ is suppressed and depends on the pulsar noise level (beside pulsar number and observation time), see Pitrou & Cusin (2025) for a detailed discussion. In all these expressions $0 \le \ell \le 2L_{\text{max}}$ and $L_1, L_2 \le L_{\text{max}}$ with $L_{\text{max}} \sim \sqrt{N_p}$, and $N_p$ is the number of pulsars in the network.

The finite resolution of the GW map also affects the reconstruction of the angular power spectrum of cross-correlation. Assuming the galaxy map to be noiseless (being its angular resolution much better than the GW background one), one finds (see Appendix B)

$$ N_{\ell}^{\text{GW}\Delta}(L_{\text{max}}) = \sqrt{\frac{C_{\ell}^{\Delta}}{2\ell + 1} \left(\frac{\bar{\Omega}_{\text{GW}}}{4\pi}\right) |B_{\ell m}|}. \qquad (30) $$

# 3 RESULTS

Before we embark on assessing the impact of cross-correlations on the detectability of the signal, we note that the weight function, $W(r)$, should be chosen so as to maximize the SNR of cross-correlations. This can be done as long as radial information (i.e. accurate redshifts) are available for all galaxies in the survey we cross-correlate with, which we will assume here. As detailed in Alonso et al. (2020), the optimal weights can be derived in terms of a Wiener filter, leading to the result

$$ W(r) = \frac{\partial_r \bar{\Omega}_{\text{GW}}(r)}{\bar{\Omega}_{\text{GW}}}. \qquad (31) $$

Physically this means that we approximately weight all galaxies by a $1/r^2$ factor, hence mimicking the properties of a background mapped in intensity.

## 3.1 Signal to noise ratio

We now give estimates of the SNR of the autocorrelation and cross-correlation. Following the derivation outlined in Appendix C, the signal to noise of the AGWB autocorrelation is given by

$$ \left(\frac{S}{N}\right)^2_{\text{GW},\ell} = 2(2\ell + 1) \left(\frac{C_{\ell}^{\text{GW}}}{C_{\ell}^{\text{GW}} + N_{\ell}^{\text{GW, tot}}}\right)^2, \qquad (32) $$

where $N_{\ell}^{\text{GW, tot}} = N_{\ell}^{\text{GW}} + N_{\ell}^{\text{GW}}(L_{\text{max}})$ while the one of the cross-correlation is given by

$$ \left(\frac{S}{N}\right)^2_{\text{GW}\Delta,\ell} = \frac{(2\ell + 1)(C_{\ell}^{\text{GW}\Delta})^2}{(C_{\ell}^{\text{GW}\Delta} + N_{\ell}^{\text{GW}\Delta,\text{tot}})^2 + (C_{\ell}^{\text{GW}} + N_{\ell}^{\text{GW, tot}})(C_{\ell}^{\Delta} + N_{\ell}^{\Delta,\text{tot}})}. \qquad (33) $$
where $N_{\ell}^{\text{GW}\Delta, \text{tot}} = N_{\text{shot}}^{\text{GW}\Delta} + N_{\ell}^{\text{GW}\Delta}(L_{\text{max}})$, and we have assumed an ideal galaxy experiment where shot noise is the only noise component. For the AGWB signal, the perfect experiment and our effective treatment of an imperfect experiment are recovered depending on the value of $N_f$ used in the auto and cross-correlation noise. We recall that $L_{\text{max}}$ corresponds to the angular resolution of the time-residual map. One can combine the various multipoles to constrain the global amplitude of the angular power spectrum. Then we define a cumulative SNR as

$$ \left(\frac{S}{N}\right)_{\text{GW}}^2 = \sum_{\ell}^{\ell_{\text{max}}} \left(\frac{S}{N}\right)_{\text{GW},\ell}^2, \qquad (34) $$

and similarly for the cross-correlation. We recall that $\ell_{\text{max}}$ corresponds to the resolution of the energy density map and $\ell_{\text{max}} = 2L_{\text{max}}$.

In a scenario with $N_f \gg 1$, in both equations (32) and (33) the dominant contribution to the denominator is coming from the variance of AGWB autocorrelation $N_{\text{shot}}^{\text{GW}}$, which appears quadratically in equation (32) and linearly in equation (33). Hence, the SNR of cross-correlation will be typically enhanced with respect to the AGWB autocorrelation one. For a finite value of $N_f$, e.g. for $N_f = 10$, with $L_{\text{max}} = 20$, the variance (equation 28) of the autocorrelation estimator is much smaller than the corresponding shot noise when no cut to the catalogue is applied. However, it becomes slightly larger with a $10^{-15}$ cut and even one order of magnitude larger for the $10^{-16}$ cut as can be checked on the left panel of Fig. 2. When considering the cross-correlation with galaxies, the variance (equation 30) is slightly larger than the corresponding shot noise for the full catalogue, and about 1 order of magnitude larger for the cut catalogues, as can be checked on the right panel of Fig. 2. The denominator of equation (33) is none the less still dominated by the product $N_{\ell}^{\text{GW}, \text{tot}} C_{\ell}^{\Delta}$, which is enhanced for the cut catalogues by the aforementioned fact that $N_{\ell}^{\text{GW}}(L_{\text{max}}) > N_{\text{shot}}^{\text{GW}}$.

## 3.2 Astrophysical modelling and numerical results

We make use of the astrophysical BBH populations of Sesana (2013) and Rosado et al. (2015), and we refer the reader to those papers for full details. In short, the BBH cosmic population is derived from observations of the galaxy mass function $\phi(M_g, z)$ and pair fraction $F(z, M_g, q_g)$, as a function of galaxy mass $M_g$, redshift $z$ and mass ratio of galaxy pairs $q_g < 1$. These observables can be combined with a merger time-scale $\tau(z, M_g, q_g)$ inferred by numerical detailed simulations to get a galaxy merger rate as:

$$ \frac{d^3 n_g}{dz dM_g dq_g} = \frac{\phi(M_g, z) F(z, M_g, q_g)}{M_g \ln 10 \tau(z, M_g, q_g)} \frac{dt_r}{dz}, \qquad (35) $$

where $t_r$ is the rest-frame cosmic time and $dt_r/dz$ converts time into redshift for a given cosmology. The BBH merger rate is then computed by populating galaxies with black holes according to observed scaling relations of the form:

$$ \log_{10} M_{\text{BH}} = \alpha + \beta \log_{10} X, \qquad (36) $$

where $X$ can be the galaxy bulge mass, or its mid-infrared luminosity or velocity dispersion $\sigma$ (see Sesana 2013 for a list of those relations).

Here, we take $\phi(M_g, z)$ from Muzzin et al. (2013), $F(z, M_g, q_g)$ from de Ravel et al. (2009), and $\tau(z, M_g, q_g)$ from Kitzbichler & White (2008). Further, we use the $M-\sigma$ relation given by Kormendy & Ho (2013), obtaining a GW background with nominal characteristic strain amplitude at the reference frequency of 1 yr $h_{1\text{yr}} = 2.4 \times 10^{-15}$ where the strain amplitude $h_c(f) = h_{1\text{yr}}(f \times \text{yr})^{-2/3}$ is related to the GW energy density by $\bar{\Omega}_{\text{GW}}(f) =$

$2\pi^2/3(f/H_0)^2 h_c^2(f)$. Using the value $1/H_0 = 14.4 \times 10^9$ yr one finds $\bar{\Omega}_{\text{GW}}(f = \text{yr}^{-1}) \simeq 8 \times 10^{-9}$, consistent with the one inferred from the EPTA DR2 new analysis (Antoniadis et al. 2023).

Given a BBH merger rate, we compute the cosmic population of inspiraling BBHs, $d^3N/(dzdM\,d\ln f)$, by assuming circular-GW driven systems. One can then sample this distribution in the parameter region ($M > 10^7 M_{\odot}) \cap (f > 10^{-9} \text{ Hz}) \cap (z < 1.3)$. One sample draw counts $\approx 100$K binaries, with chirp masses, redshift, and GW emission frequency, that are used to construct the GW background energy density as described in Appendix D.<fnref n="4" />

Finally, in order to compute the statistical properties of the GW background, we also need to choose how to distribute the GW sources across the sky. We assume that the distribution of galaxies follows the distribution of dark matter haloes with bias factor that we assume to be scale-independent and with redshift evolution given by $b(z) = b_0 \sqrt{1+z}$ and $b_0 = 1.5$ (Rassat et al. 2008; Marin et al. 2013). Finally, as mentioned earlier, we assume that the clustering properties of the massive BBH population are identical to those of the underlying galaxy distribution. This is, of course, a simplifying approximation, as we expect more massive galaxies to be more likely to host massive black holes. In other words, the population of massive BBH is expected to be more clustered than the galaxy distribution. Since refining this bias model would likely strengthen the clustering signal, our simplifying assumption is conservative for the purposes of this study. We also use a constant comoving galaxy density $a^3 n_g \sim 0.1 \text{ Mpc}^{-3}$.

As already mentioned, we expect the AGWB strain map to be dominated by a few very bright sources: by resolving them, anisotropies of the AGWB generated by the remaining unresolved population are a representative sample of large-scale structure anisotropies. We filter out from the background budget sources whose strain is higher than the current EPTA sensitivity curve (Agazie et al. 2024), roughly corresponding to $h_{\text{eff}} > 10^{-14}$ around the sweet spot of the sensitivity curve. We have here introduced an effective strain $h_{\text{eff}}$ defined by

$$ h_{\text{eff}}^2(f) = h^2(f) f T_{\text{obs}}, \qquad (37) $$

where $T_{\text{obs}} = 16$ yr and $h(f) = 8/(\sqrt{10}D_L)G M_z [\pi G M_z f]^{2/3}$ with $M_z = (1+z)M$ is the polarization and inclination averaged strain (see e.g. Sesana et al. 2008 with proper factors of redshift accounted for). As the sensitivity of a given PTA observatory will improve, more and more sources will be detectable individually: when deriving forecasts for future experiments, the filtering should be performed following the iterative procedure of Truant et al. (2025), accounting in the computation of the expected noise budget for both a pure noise contribution and contribution of the unresolvable background component.<fnref n="5" /> Since our goal is just to quantify the effect of filtering without targeting a specific future PTA mission, we test the effect of moving the threshold on $h_{\text{eff}}$, filtering out sources with $h_{\text{eff}} > 10^{-15}$ and with $h_{\text{eff}} > 10^{-16}$, the latter corresponding to the plateau of the SKA sensitivity curve (Truant et al. 2025).<fnref n="6" /> We refer to these three values for the filtering threshold as ‘no cut’ (as the

<footnote n="4">⁴A Monte Carlo sampling of the population returns correctly a discrete number of sources, while the semi-analytical result contains spurious contributions from fractional sources that are clearly not present (Sesana, Vecchio & Colacino 2008). In other words, sampling is needed because when using an analytical formula, one would also need to add a ‘fractions of sources,’ while Monte Carlo methods avoid this issue.</footnote>

<footnote n="5">⁵We stress that both these contributions are accounted for when computing the sensitivity curve using real PTA data.</footnote>

<footnote n="6">⁶We stress again that in the context of SKA a proper filtering procedure would need to account not only for the instrumental noise component, but also for</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Background energy density $\hat{\Omega}_{\text{GW}}(f)$ as function of frequency, for the three different cut-off values used in the filtering procedure.</figcaption></figure>

EPTA sensitivity effectively does not allow to resolve any individual source), ‘$10^{-15}$ cut’ and ‘$10^{-16}$ cut’.

In Fig. 1, we show the background energy density as function of frequency, for the three aforementioned cut-off values. In the ‘no cut’ situation, after averaging over modulations, the scaling is roughly $\propto f^{2/3}$ as expected. The bumps in the ‘no cut’ curve are due to the contribution of a few very bright sources that dominate the total signal at high frequencies (plotting source contributions individually, they would appear as spikes at a given frequency). We observe that the filtering procedure has the effect of eliminating these high-frequency bumps and changing the low-frequency scaling. Indeed for a given threshold $h_{\text{eff}}^{\text{cut}}$, we infer from equation (37) that $M_{\text{max}} \propto (h_{\text{eff}}^{\text{cut}})^{3/5} f^{-7/10}$, hence the larger the frequency, the more sources are removed from the catalogue, leaving behind a smaller residual amplitude.

Fig. 2 represents for the three situations the angular power spectrum of auto and cross-correlations, and the corresponding shot noise contributions.⁷ We see that while the autocorrelation map is shot noise dominated ($N_{\text{shot}}^{\text{GW}} \gg C_{\ell}^{\text{GW}}$), the shot noise contribution to the cross-correlation is much suppressed ($N_{\text{shot}}^{\text{GW}\Delta} \ll C_{\ell}^{\text{GW}\Delta}$). The finite resolution noise as well is much suppressed for the cross-correlation than for the autocorrelation. We observe that the filtering procedure reduces both the angular power spectrum of the clustering signal and the noise components, for both auto and cross-angular power spectra.

In Fig. 3 (left panel), we show the SNR for auto and cross-correlations (for the three cut-off values) as a function of the multipole $\ell$. On the right panel, we show the corresponding cumulative SNR, as a function of $\ell_{\text{max}}$. We consider both a scenario with $N_f = 10$ and an ideal scenario where shot noise is the only noise component (bottom plots). As expected, the SNR of the cross-correlation outperforms the autocorrelation one (by more almost two orders of magnitudes in the absence of a cut-off). This is due to the aforementioned hierarchy between signal and noise components along with $N_{\text{shot}}^{\Delta} \ll C_{\ell}^{\Delta}$, i.e. the fact that the autocorrelation SNR scales as $\sqrt{2(2\ell + 1)}C_{\ell}^{\text{GW}}/N_{\ell}^{\text{GW, tot}}$ whereas for cross-correlation it scales as $\sqrt{(2\ell + 1)}C_{\ell}^{\text{GW}\Delta}/\sqrt{C_{\ell}^{\Delta} N_{\ell}^{\text{GW, tot}}}$.

We observe that even in the ideal case where the instrumental noise is set to zero (such that $N_{\ell}^{\text{GW, tot}} = N_{\text{shot}}^{\text{GW}}$), the SNR of single

multipoles is never much larger than 1, unless a very low threshold of detection of single events is chosen. However, the cumulative SNR is above 1 already at $\ell_{\text{max}} = 4$ (corresponding to $\sim 25$ pulsars), for a filtering threshold of $10^{-16}$, and it is above 1 for $\ell_{\text{max}} = 13$ for a threshold of $10^{-15}$. This indicates that, while the reconstruction of the shape of the angular power spectrum will be challenging, it will be potentially possible to constrain the global amplitude of the spectrum, by combining the information from different multipoles. For a realistic experiment (i.e. $N_f$ finite) the situation stays qualitatively the same: when $N_f = 10$, detectability can only be achieved using multipoles up to $\ell_{\text{max}} = 30$, with a threshold of $10^{-15}$.

We observe that to obtain these results, we have fixed the value of the comoving galaxy number density equal to $a^3n_G = 0.1 \text{ Mpc}^{-3}$, however our results for the SNR depend very mildly on this choice. Indeed the SNR of the autocorrelation map does not depend on the galaxy number density assumed, while, for any reasonable value of $a^3n_G$, the dominant contribution to the noise of the cross-correlation map comes from the $\propto N_{\ell}^{\text{GW, tot}}C_{\ell}^{\Delta}$ contribution in the denominator of equation (32), which does not depend on the galaxy number density.

Finally, being the resolution of a PTA map set by the number of pulsars in the network, i.e. $L_{\text{max}} \sim \sqrt{N_{\text{p}}}$, when correlating a PTA map with a galaxy map, it is the former which sets the limiting resolution (current galaxy maps already reach $\ell_{\text{max}} > 10^3$, and the situation will further improve with next generation surveys such as Euclid). The effect of masking is expected to affect our ability to constrain large angular scales. However, as can be seen in Fig. 3, left panel, the SNR has a plateau at small $\ell$, hence even when removing the very few multipoles, our final findings will not be drastically affected.

## 4 DISCUSSION AND CONCLUSIONS

In this article, we explored the potential of cross-correlation between a background map and a dense galaxy map to obtain a first detection of background anisotropies. We considered both the case of a perfect instrument (for which pulsar white noise is negligible), where shot noise is the dominant source of noise, and the case of a realistic survey in which an additional noise component arises due to the finite resolution of a time residual map (finite resolution noise). We derived from first principles the expressions of shot noise and time resolution noise for both the auto and cross-correlation spectrum. The magnitude of the finite resolution noise is set by the effective number of frequencies, which in turn depends on the number of pulsar monitored, pulsar white noise level and observation time, as discussed in detail in Pitrou & Cusin (2025). For a perfect experiment (pulsar white noise is negligible) one has $N_f = 2Tf_{\text{max}} \gg 1$, while for a realistic survey we consider $N_f = 10$.

Using catalogues of massive BBH from Sesana (2013) and Rosado et al. (2015), we showed that, unlike the autocorrelation, the cross-correlation map is not shot noise dominated: for an ideal experiment, the SNR of the cross-correlation outperforms the autocorrelation by almost two orders of magnitude. We also tested the effect of a filtering procedure that removes the brightest sources from the total background: this procedure further reduces shot noise, as the brightest sources dominate the background budget, obscuring the contribution of the unresolved population that traces the large-scale structure. A similar result was recently found in Semenzato et al. (2024), relying on numerical simulations. However, while our findings are consistent with the simulation results of Semenzato et al. (2024), our analytic treatment disagrees with the explanation

the background budget. In this respect, a threshold $h_{\text{eff}} > 10^{-16}$ for filtering out resolvable events with SKA is optimistic.

⁷Our results for the cross-spectrum are compatible in order of magnitude with those of Sah & Mukherjee (2024) (fig. 3), where a different set of astrophysical models is used.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Angular power spectrum of autocorrelation (left panel) and cross-correlation (right panel) and corresponding shot noise and finite resolution map noise contributions. The signal is integrated over the PTA frequency range. We are here using $L_{max} = 20$ and $N_f = 10$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Left: SNR as function of the multipole $\ell$ of cross-correlation (continuous lines) and of autocorrelation (dashed lines) for $N_f = 10$ and for an ideal experiment in which the shot noise is the only noise component. Right: corresponding cumulative SNR as function of $\ell_{max}$. The first line of panels uses $N_f = 10$ whereas the second line is idealized with $N_f = \infty$. Both cases consider an angular resolution set by $L_{max} = 20$.</figcaption></figure>

proposed in that reference.<fnref n="8" /> We show that the finite resolution noise is as well considerably reduced in the cross-correlation.

For both the cases of a perfect experiment and a realistic one, when filtering out bright sources, we had to assume a threshold value to determine whether a source is resolvable individually. To do so, we considered the current EPTA sensitivity curve (effectively resulting in no resolvable sources), a threshold value corresponding to SKA sensitivity, and an intermediate scenario. We found that even in the ideal case where the instrumental noise is set to zero (the noise is just shot noise), the SNR of single multipoles is never much larger

than 1, unless a very low threshold of detection of single events is chosen. However, the cumulative SNR is above 1 already at $\ell_{max} = 4$ (corresponding to ~ 25 pulsars), for a filtering threshold of $10^{-16}$, and it is above 1 for $\ell_{max} = 13$ for a threshold of $10^{-15}$. This indicates that, while the reconstruction of the shape of the angular power spectrum will be challenging, it will be potentially possible to constrain the global amplitude of the spectrum, by combining the information from different multipoles. For a realistic experiment (i.e. $N_f$ finite) the situation stays qualitatively the same: when $N_f = 10$, detectability can be achieved using multipoles up to $\ell_{max} = 30$, with a threshold of $10^{-15}$. We stress again that the case of $N_f = 10$ is rather futuristic: as shown in detail in Pitrou & Cusin (2025) the number of effective frequencies $N_f$ is in fact a decreasing function of the map multipole $L$. Assuming the pulsar noise characteristics of Agazie et al. (2023b), getting $N_f(L = 5) \sim 1$ would require 450 pulsars monitored during

<footnote n="8">⁸Through a detailed analytic derivation, we demonstrate that shot noise is, in fact, never fully eliminated in cross-correlation maps, contrary to what is claimed in equation (7) of Semenzato et al. (2024), which states that the cross-correlation map is shot noise free.</footnote>
<table><caption>**Table 1.** Estimation of $N_f(L)$ with the method of Pitrou & Cusin (2025). The first line corresponds to the set of pulsars of Agazie et al. (2023b), but assuming they are all observed over 25 yr. The second line assumes to have 10 times more pulsars with the same characteristics. The last two lines revisit these scenarios, with all white noises dispersions reduced by a factor 10.</caption><thead><tr><th></th><th>N<sub>f</sub> (L = 2)</th><th>N<sub>f</sub> (L = 3)</th><th>N<sub>f</sub> (L = 4)</th><th>N<sub>f</sub> (L = 5)</th><th>N<sub>f</sub> (L = 10)</th></tr></thead><tbody><tr><td>68 pulsars, 25 yr</td><td>8.6</td><td>1.1</td><td>0.15</td><td>0.03</td><td>0.0002</td></tr><tr><td>680 pulsars, 25 yr</td><td>31</td><td>14</td><td>5.4</td><td>1.8</td><td>0.016</td></tr><tr><td>68 pulsars, 25 yr, improved white noise</td><td>21</td><td>2.3</td><td>0.3</td><td>0.6</td><td>0.0003</td></tr><tr><td>680 pulsars, 25 yr, improved white noise</td><td>83</td><td>35</td><td>12</td><td>3.9</td><td>0.03</td></tr></tbody></table>

25 yr. In Table 1, we collect the estimations of the values of $N_f(L)$ for current and futuristic scenarios. Of course the situation would improve significantly when pulsar white noise is reduced with respect to the model of Agazie et al. (2023b), for example, by increasing the cadence of pulsar observations or by employing new instrumentation. For example, the white noise of pulsars collected in the first data releases by the South African collaboration MeerKAT (Miles et al. 2023) and the Chinese Five-hundred-meter Aperture Spherical Radio Telescope (FAST) (Xu et al. 2023) is already a factor of a few lower than what measured by other telescopes. In particular, FAST displays a median white noise level of about 100 ns across an array of more than 50 pulsars (Xu et al. 2023), and SKA is anticipated to further improve on those numbers. A future study will be dedicated to deriving more realistic forecasts for future PTA experiments and galaxy surveys such as Euclid, properly factoring in the effects of white noise reduction in future PTA surveys and realistic filtering procedures.

## ACKNOWLEDGEMENTS

The work of GC and CP was supported by CNRS. The work of GC and MP has received financial support by the SNSF Ambizione grant PZ00P2_193292 Gravitational wave propagation in the clustered universe. We are very grateful to Marc Kamionkowski for having pointed out that finite resolution noise was missing in our final results. We thank Camille Bonvin and Nastassia Grimm for interesting discussions at different stages of this work. AS acknowledges financial support by the European Research Council Consolidator Grant ‘B Massive’ (Grant Agreement: 818691) and European Research Council Advanced Grant ‘PINGU’ (Grant Agreement: 101142079).

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Agazie G., et al., 2023a, ApJ, 951, L8

* Agazie G., et al., 2023b, ApJ, 951, L9

* Agazie G., et al., 2023c, ApJ, 956, L3

* Agazie G., et al., 2024, ApJ, 966, 105

* Agazie G., et al., 2025, ApJ, 985, 99

* Allen B., 2024, Phys. Rev. D, 110, 043043

* Allen B., Agarwal D., Romano J. D., Valtolina S., 2024, Phys. Rev. D, 110, 123507

* Alonso D., Cusin G., Ferreira P. G., Pitrou C., 2020, Phys. Rev. D, 102, 023002

* Alonso D., Nikjoo M., Renzini A. I., Bellini E., Ferreira P. G., 2024, Phys. Rev. D, 110, 103544

* Anholm M., Ballmer S., Creighton J. D. E., Price L. R., Siemens X., 2009, Phys. Rev. D, 79, 084030

* Antoniadis J. et al., 2023, A&A, 678, A48

* Bernardo R. C., Ng K.-W., 2022, J. Cosmol. Astropart. Phys., 2022, 046

* Bernardo R. C., Ng K.-W., 2023, Phys. Rev. D, 107, 044007

* Bernardo R. C., Ng K.-W., 2025, Int. J. Mod. Phys. D, 34, 2540013

* Bernardo R. C., Liu G.-C., Ng K.-W., 2024, J. Cosmol. Astropart. Phys., 2024, 034

* Cruz N. M. J., Malhotra A., Tasinato G., Zavala I., 2024a, preprint (arXiv:2412.14010)

* Cruz N. M. J., Malhotra A., Tasinato G., Zavala I., 2024b, Phys. Rev. D, 110, 063526

* Cusin G., Tasinato G., 2022, J. Cosmol. Astropart. Phys., 2022, 036

* Cusin G., Pitrou C., Uzan J.-P., 2017, Phys. Rev. D, 96, 103019

* Cusin G., Pitrou C., Uzan J.-P., 2018, Phys. Rev. D, 97, 123527

* Cusin G., Durrer R., Ferreira P. G., 2019a, Phys. Rev. D, 99, 023534

* Cusin G., Dvorkin I., Pitrou C., Uzan J.-P., 2019b, Phys. Rev. D, 100, 063004

* de Ravel L. et al., 2009, A&A, 498, 379

* Depta P. F., Domcke V., Franciolini G., Pieroni M., 2025, Phys. Rev. D, 111, 083039

* Falxa M., Leclere H. Q., Sesana A., 2025, Phys. Rev. D, 111, 023047

* Gair J., Romano J. D., Taylor S., Mingarelli C. M. F., 2014, Phys. Rev. D, 90, 082001

* Grimm N., Pijnenburg M., Cusin G., Bonvin C., 2025a, J. Cosmol. Astropart. Phys., 2025, 011

* Grimm N., Pijnenburg M., Cusin G., Bonvin C., 2025b, J. Cosmol. Astropart. Phys., 2025, 047

* Grunthal K., et al., 2024, MNRAS, 536, 1501

* Hotinli S. C., Kamionkowski M., Jaffe A. H., 2019, Open J. Astrophys., 2, 8

* Jaffe A. H., Backer D. C., 2003, ApJ, 583, 616

* Jenet F. A., Romano J. D., 2015, Am. J. Phys., 83, 635

* Jenkins A. C., Sakellariadou M., 2019, Phys. Rev. D, 100, 063508

* Jenkins A. C., Romano J. D., Sakellariadou M., 2019, Phys. Rev. D, 100, 083501

* Kitzbichler M. G., White S. D. M., 2008, MNRAS, 391, 1489

* Kormendy J., Ho L. C., 2013, ARA&A, 51, 511

* Lamb W. G., Taylor S. R., 2024, ApJ, 971, L10

* Maggiore M., 2007, Gravitational Waves. Vol. 1: Theory and Experiments. Oxford Master Series in Physics. Oxford Univ. Press, Oxford

* Marin F. A., et al., 2013, MNRAS, 432, 2654

* Miles M. T. et al., 2023, MNRAS, 519, 3976

* Mingarelli C. M. F., Sidery T., Mandel I., Vecchio A., 2013, Phys. Rev. D, 88, 062005

* Muzzin A. et al., 2013, ApJ, 777, 18

* Nay J., Boddy K. K., Smith T. L., Mingarelli C. M. F., 2024, Phys. Rev. D, 110, 044062

* Pitrou C., Cusin G., 2025, Phys. Rev. D, 111, 083544

* Pitrou C., Cusin G., Uzan J.-P., 2020, Phys. Rev. D, 101, 081301

* Pol N., Taylor S. R., Romano J. D., 2022, ApJ, 940, 173

* Qin W., Boddy K. K., Kamionkowski M., Dai L., 2019, Phys. Rev. D, 99, 063002

* Rajagopal M., Romani R. W., 1995, ApJ, 446, 543

* Rassat A. et al., 2008, preprint (arXiv:0810.0003)

* Reardon D. J. et al., 2023, ApJ, 951, L6

* Roebber E., Holder G., 2017, ApJ, 835, 21

* Romano J. D., Allen B., 2024, Class. Quantum Gravity, 41, 175008

* Romano J. D., Cornish N. J., 2017, Living Rev. Relativ., 20, 2
Rosado P. A., Sesana A., 2014, MNRAS, 439, 3986

Rosado P. A., Sesana A., Gair J., 2015, MNRAS, 451, 2417

Sah M. R., Mukherjee S., 2024, preprint (arXiv:2407.11669)

Sato-Polito G., Kamionkowski M., 2024, Phys. Rev. D, 109, 123544

Semenzato F., Casey-Clyde J. A., Mingarelli C. M. F., Raccanelli A., Bellomo N., Bartolo N., Bertacca D., 2024, preprint (arXiv:2411.00532)

Sesana A., 2013, MNRAS, 433, L1

Sesana A., Vecchio A., Colacino C. N., 2008, MNRAS, 390, 192

Siemens X., Ellis J., Jenet F., Romano J. D., 2013, Class. Quantum Gravity, 30, 224015

Tian C., Ding R., Kou X.-X., 2024, preprint (arXiv:2412.01219)

Truant R. J., Izquierdo-Villalba D., Sesana A., Mohiuddin Shaifullah G., Bonetti M., 2025, A&A, 694, A282

Wu Y.-M., Bi Y.-C., Huang Q.-G., 2024

Xu H. et al., 2023, Res. Astron. Astrophys., 23, 075024

Yang K. Z., Suresh J., Cusin G., Banagiri S., Feist N., Mandic V., Scarlata C., Michaloliakos I., 2023, Phys. Rev. D, 108, 043025

## APPENDIX A: GENERAL PROPERTIES OF COMPOUND STATISTICS

Let us consider $N$ being a random variable with variance $\langle N^2 \rangle - \langle N \rangle^2$ and mean $\langle N \rangle$. Then for a generic function $f(N)$

$$ \langle f(N) \rangle = \int f(N) p(N) dN. \quad (A1) $$

We are interested in the statistics of the compound random variable $Y$ given by

$$ Y = \sum_{i=1}^{N} X_i, \quad (A2) $$

where $X_i$ is either 1 or 0, i.e. it is a binomial variable $X_i \sim B(1, \beta)$. It has a probability $p(Y)$, such that

$$ p(Y) = \int p(Y|N)p(N)dN. \quad (A3) $$

The conditional probability $p(Y|N)$ can be found from the fact that it is the sum of $N$ random variables (it is a convolution of probabilities), but we do not need to go to this level of sophistication to get convinced that

$$ \langle Y \rangle_N = \int Y p(Y|N) dY = N \langle X \rangle. \quad (A4) $$

The subindex on the left hand side means that an average at $N$ being fixed. For the second moment one has

$$ \langle Y^2 \rangle_N = \int Y^2 p(Y|N) dY = N \langle X^2 \rangle + N(N-1) \langle X \rangle^2, \quad (A5) $$

which is obtained by separation of the $X_i^2$ and the $X_i X_j$ with $i \neq j$.

We can work out the first moments of $Y$ from $p(Y)$. For a generic function of $Y$, we have

$$ \begin{aligned} \langle f(Y) \rangle &= \int f(Y) p(Y) dY \\ &= \iint f(Y) p(Y|N) p(N) dN dY \\ &= \int p(N) dN \left( \int f(Y) p(Y|N) dY \right) \\ &= \int \langle f(Y) \rangle_N p(N) dN, \end{aligned} \quad (A6) $$

where in the last line the remaining randomness of $N$ has to be treated. The average $\langle Y \rangle$ can be easily computed

$$ \langle Y \rangle = \int \langle Y \rangle_N p(N) dN = \int \langle X \rangle N p(N) dN = \langle N \rangle \langle X \rangle. \quad (A7) $$

Similarly for the second moment

$$ \begin{aligned} \langle Y^2 \rangle &= \int \langle Y^2 \rangle_N p(N) dN \\ &= \int (N \langle X^2 \rangle + N(N-1) \langle X \rangle^2) p(N) dN \\ &= \langle N \rangle (\langle X^2 \rangle - \langle X \rangle^2) + \langle N^2 \rangle \langle X \rangle^2. \end{aligned} \quad (A8) $$

Hence we obtain the variance of $Y$

$$ \begin{aligned} \langle Y^2 \rangle - \langle Y \rangle^2 &= \langle N \rangle (\langle X^2 \rangle - \langle X \rangle^2) + (\langle N^2 \rangle - \langle N \rangle^2) \langle X \rangle^2 \\ &= \langle N \rangle \langle X^2 \rangle, \end{aligned} \quad (A9) $$

where the last equality holds when $N$ is a Poisson variable. This derivation is valid whatever the statistics for $X$, and $\langle X^2 \rangle = \langle X \rangle$ if $X \sim B(1, \beta)$. Following the same kind of reasoning we can also get cross-correlations between $Y$ and $N$, and

$$ \begin{aligned} \langle YN \rangle &= \int \langle Y \rangle_N N p(N) dN \\ &= \int \langle X \rangle N^2 p(N) dN \\ &= \langle X \rangle \langle N^2 \rangle \\ &= \langle X \rangle (\langle N \rangle + \langle N \rangle^2), \end{aligned} \quad (A10) $$

where the last equality holds when $N$ is Poisson variable. Then

$$ \langle YN \rangle - \langle Y \rangle \langle N \rangle = \langle X \rangle \mathrm{Var}(N) = \langle X \rangle \langle N \rangle, \quad (A11) $$

where the last equality holds when $N$ is a Poisson variable.

## APPENDIX B: DERIVATION OF THE FINITE RESOLUTION NOISE

This appendix is a summary of the results of Pitrou & Cusin (2025) on anisotropic estimators for a stochastic background. The goal is to derive the variance affecting our ability to construct the angular power spectrum of the background energy density, due to the finite resolution of the time-residual map.

The statistics of the multipoles $\hat{z}_{LM}(f)$ of pulsar redshifts is given by

$$ \langle \hat{z}_{L_1 M_1} (f_1) \hat{z}_{L_2 M_2}^*(f_2) \rangle_A = \delta_{L_1 L_2} \delta_{M_1 M_2} \delta(f_1 - f_2) C_{L_1}^z(f_1), \quad (B1) $$

where the average $\langle \dots \rangle_A$ is over the GW realization and

$$ C_L^z(f) = \hat{P}(f) C_L^{HD}, \quad (B2) $$

where $C_L^{HD}$ is the multipole decomposition of the Hellings-Downs curve while $\hat{P}(f)$ is the (two-sided, meaning that needs to be integrated on both positive and negative frequencies) spectral function, related to the gravitons distribution function as

$$ \hat{P}(f) = 4\pi f^2 \mathcal{A}(2\pi f). \quad (B3) $$

The isotropic part of the spectral density is related to the background energy density, which can be quantified in units of the critical energy density with

$$ \bar{\Omega}_{\mathrm{GW}}(f) = \frac{8\pi^2}{3H_0^2} f^3 \hat{P}(f) = \frac{32\pi^3}{3H_0^2} f^5 \mathcal{A}(2\pi f). \quad (B4) $$

In full generality, the graviton distribution function can depend on directions. Let us parametrize its angular dependence as

$$ \mathcal{A}(\hat{k}) = \mathcal{A}(k) \hat{A}(\hat{k}), \quad \hat{A}(\hat{k}) = \sum_{\ell m} \hat{A}_{\ell m} Y_{\ell}^{m}(\hat{k}), \quad (B5) $$

where $k = |\hat{k}|$ and $\hat{A}_{\ell m} = (-1)^m \hat{A}_{-\ell m}$ since the GW spectrum is real valued. This decomposition assumes the existence of a factorization
between a spectral shape $\mathcal{A}(k)$, and an angular-dependent function $\hat{\mathcal{A}}(\hat{k})$. We get from equations (1) and (B4)

$$ (\delta\Omega_{\text{GW}})_{\ell m} = \frac{8\pi^2}{3H_0^2} f^5 \mathcal{A}(2\pi f) \hat{\mathcal{A}}_{\ell m}. \quad (B6) $$

Let us discuss the effects of an anisotropic spectrum: if $\hat{\mathcal{A}}(\hat{k}) \neq 1$, the correlation in equation (B1) acquires off-diagonal contributions, and one obtains Pitrou & Cusin (2025)

$$ \begin{aligned} \langle \hat{z}_{L_1 M_1}(f_1) \hat{z}_{L_2 M_2}(f_2) \rangle_A = & \delta(f_1 + f_2) \hat{\mathcal{P}}(f_1) F_{L_1} F_{L_2} \\ & \cdot \sum_{\ell m} \delta_{L_1+L_2+\ell}^{\text{even}} \frac{\hat{\mathcal{A}}_{\ell m}^*}{\sqrt{4\pi}} \sqrt{2\ell+1} \begin{pmatrix} \ell & L_1 & L_2 \\ m & M_1 & M_2 \end{pmatrix} \begin{pmatrix} \ell & L_1 & L_2 \\ 0 & -2 & 2 \end{pmatrix}, \end{aligned} \quad (B7) $$

where the quantities $F_L$ are defined as

$$ F_L \equiv \sqrt{4\pi(2L+1)} \sqrt{\frac{(L-2)!}{(L+2)!}}, \quad (B8) $$

and $\delta_{L_1+L_2+\ell}^{\text{even}} = 1$ if $L_1 + L_2 + \ell$ is even and vanishes if it is odd. The isotropic result (equation B1) is easily recovered when using $\hat{\mathcal{A}}_{\ell m} = \delta_\ell^0 \delta_m^0 \sqrt{4\pi}$. The HD correlation is a signature of an isotropic spectral function $\mathcal{A}(k)$. The anisotropic part $\hat{\mathcal{A}}(\hat{k})$ of the GW background is reflected instead in the off-diagonal correlations in the $(L, M)$ space.

Our goal is to isolate the $\hat{\mathcal{A}}_{\ell m}$ modes and this has to be done by exploiting the orthonormality of Wigner 3j. Hence, we introduce the combinations (contractions of the 2-point function with Wigner 3j symbol)

$$ \begin{aligned} \mathcal{E}_{L_1 L_2}^{\ell m}(f_1, f_2) = & \sqrt{2\ell + 1} \sum_{M_1 M_2} \begin{pmatrix} \ell & L_1 & L_2 \\ m & M_1 & M_2 \end{pmatrix} \\ & \cdot \hat{z}_{L_1 M_1}(f_1) \hat{z}_{L_2 M_2}(f_2), \end{aligned} \quad (B9) $$

whose average over the GW realization can be computed using equation (B7) and using properties of Wigner symbols, and is given by

$$ \begin{aligned} \langle \mathcal{E}_{L_1 L_2}^{\ell m}(f_1, f_2) \rangle_A = & \delta_{L_1+L_2+\ell}^{\text{even}} \delta(f_1 + f_2) \hat{\mathcal{P}}(f_1) \\ & \cdot \frac{\hat{\mathcal{A}}_{\ell m}^*}{\sqrt{4\pi}} F_{L_1} F_{L_2} \begin{pmatrix} \ell & L_1 & L_2 \\ 0 & -2 & 2 \end{pmatrix}. \end{aligned} \quad (B10) $$

Hereafter, we consider only parameters such that $L_1 + L_2 + \ell$ is even. In the case where the anisotropy is mild, we can assume the covariance is only due to the isotropic part. In that case, the covariance takes the simple form

$$ \begin{aligned} \text{Var}(\mathcal{E}_{L_1 L_2}^{\ell m}) = & \langle \mathcal{E}_{L_1 L_2}^{\ell m}(f_1, f_2) \mathcal{E}_{L_1^* L_2^*}^{\ell' m^*}(f_1', f_2') \rangle_A \\ & - \langle \mathcal{E}_{L_1 L_2}^{\ell m}(f_1, f_2) \rangle_A \langle \mathcal{E}_{L_1^* L_2^*}^{\ell' m^*}(f_1', f_2') \rangle_A \\ & \approx \delta_\ell^\ell \delta_m^{m'} \left[ \delta_{L_1}^{L_1'} \delta_{L_2}^{L_2'} \delta(f_1 - f_1') C_{L_1}^z(f_1) \delta(f_2 - f_2') C_{L_2}^z(f_2) \right. \\ & \left. + \delta_{L_1}^{L_2'} \delta_{L_2}^{L_1'} \delta(f_1 - f_2') C_{L_1}^z(f_1) \delta(f_2 - f_1') C_{L_2}^z(f_2) \right]. \end{aligned} \quad (B11) $$

If we want to combine the expressions (equation B9) to form an optimal estimator for the $\hat{\mathcal{A}}_{\ell m}$, we must normalize them and weight them by their inverse variance (Hotinli et al. 2019). Explicitly, after having discretized the frequency space with bins of width $1/T$, one finds the following estimator

$$ (\hat{\mathcal{A}}_{\ell m, f}^{L_1 L_2})^e = \frac{\sqrt{4\pi}}{T F_{L_1} F_{L_2}} \frac{1}{\hat{\mathcal{P}}(f)} \begin{pmatrix} \ell & L_1 & L_2 \\ 0 & -2 & 2 \end{pmatrix}^{-1} \mathcal{E}_{L_1 L_2}^{\ell m} (f, -f)^*, \quad (B12) $$

where on the left-hand side we have explicitly indicated that one can build an estimator for each given value of $L_1$ and $L_2$. From equation (B10) the mean of this object is trivially given by $\hat{\mathcal{A}}_{\ell m}^{L_1 L_2}$ as we considered only cases for which $L_1 + L_2 + \ell$ is even. Assuming that anisotropies are mild, and using equation (B11), the variance of this object for $f \neq 0$ is given by

$$ \text{Var} \left( \hat{\mathcal{A}}_{\ell m, f}^{L_1 L_2} \right)^e = \frac{4\pi}{(2L_1 + 1)(2L_2 + 1)} \begin{pmatrix} \ell & L_1 & L_2 \\ 0 & -2 & 2 \end{pmatrix}^{-2}. \quad (B13) $$

We then combine all the estimators $\hat{\mathcal{A}}_{\ell m, f}^{L_1 L_2}$ with inverse variance weighting to obtain the minimum-variance estimator. However since

$$ (\hat{\mathcal{A}}_{\ell m, f}^{L_1 L_2})^e = (\hat{\mathcal{A}}_{\ell m, -f}^{L_1 L_2})^e, \quad (B14) $$

we restrict to $f > 0$ while considering all $\ell_1, \ell_2$ combinations so as to avoid double counting. The estimator is therefore

$$ \hat{\mathcal{A}}_{\ell m}^e = \frac{\sum_{L_1 L_2} \sum_{f>0} \delta_{L_1+L_2+\ell}^{\text{even}} \left(\hat{\mathcal{A}}_{\ell m, f}^{L_1 L_2}\right)^e \text{Var}^{-1}\left(\hat{\mathcal{A}}_{\ell m, f}^{L_1 L_2}\right)}{\sum_{L_1 L_2} \sum_{f>0} \delta_{L_1+L_2+\ell}^{\text{even}} \text{Var}^{-1}\left(\hat{\mathcal{A}}_{\ell m, f}^{L_1 L_2}\right)}. \quad (B15) $$

The variance of this estimator is computed from

$$ \text{Var}^{-1}(\hat{\mathcal{A}}_{\ell m}^e) = \sum_{L_1 L_2} \sum_{f>0} \delta_{L_1+L_2+\ell}^{\text{even}} \text{Var}^{-1}(\hat{\mathcal{A}}_{\ell m}^{L_1 L_2})^e, \quad (B16) $$

and we get

$$ \begin{aligned} \text{Var}^{-1}(\hat{\mathcal{A}}_{\ell m}^e) = & \frac{1}{8\pi} \sum_{L_1 L_2} N_f \delta_{L_1+L_2+\ell}^{\text{even}} (2L_1+1)(2L_2+1) \\ & \times \begin{pmatrix} \ell & L_1 & L_2 \\ 0 & -2 & 2 \end{pmatrix}^2, \end{aligned} \quad (B17) $$

where $N_f \equiv \min(N_f(L_1), N_f(L_2))$ and the sum is over all values of $L_1$ and $L_2$. In this equation, we assumed that in the absence of instrumental noise all frequencies equally contribute, and $N_f$ is equal to the number of frequency modes, i.e. $N_f = 2T f_{\text{max}}$ for a perfect experiment, so that $N_f/2$ is the number of positive frequency modes. In the presence of noise, the effective number of frequency modes has to be evaluated following the procedure described in section VII C2 of Pitrou & Cusin (2025), and is a function of the map multipole $L$, which depends on the pulsar noise level (beside pulsar number and observation time), see Pitrou & Cusin (2025) for a detailed discussion. Hence here $N_f$ represents an effective quantity, which we introduce for practical purposes, even if a realistic evaluation for a given instrument should keep the dependence on $L_1$ and $L_2$. This result (equation B17) agrees with what obtained in Hotinli et al. (2019) (where it is implicitly assumed that $N_f = 1$). Notice that all these sums run over $0 \le \ell \le 2L_{\text{max}}$ and $L_1, L_2 \le L_{\text{max}}$ with $L_{\text{max}} \sim \sqrt{N_p}$, and $N_p$ is the number of pulsars in the network.

We conclude this section by relating this anisotropy estimator to the estimator of the energy density angular power spectrum, which is usually the observables used to quantify anisotropies (Cusin et al. 2017). Using equations (B3), (B4), and (1) we build the estimator for the multipoles defined in equation (6a)

$$ a_{\ell m}^{e} = \frac{8\pi^2}{3H_0^2} f^5 \mathcal{A}(2\pi f) \hat{\mathcal{A}}_{\ell m}^{e}, \quad (B18) $$

whose mean is $a_{\ell m}$ by construction, and whose variance can be trivially computed from the one of $\hat{\mathcal{A}}_{\ell m}^e$. An estimator for the angular
power spectrum of this quantity is given by

$$C_{\ell}^{\text{GWe}} = \frac{1}{2\ell + 1} \sum_{m} a_{\ell m}^{e} a_{\ell m}^{e*}, \qquad (\text{B19})$$

with variance given by

$$\text{Var}(C_{\ell}^{\text{GWe}}) = \frac{2}{(2\ell + 1)^2} \sum_{m} [\text{Var}^4(a_{\ell m}^e) + 2\text{Var}(a_{\ell m}^e a_{\ell m}^*)]. \qquad (\text{B20})$$

We stress again that the average is an average over GW realizations, and there is a residual stochasticity due to cosmological initial conditions (see also Grimm et al. 2025b for a detailed discussion of this double average procedure). In other words, the result above holds for any cosmological realization. To obtain predictions, one has to take an additional average over cosmological initial conditions and we get

$$\text{Var}(C_{\ell}^{\text{GWe}}) = \frac{2}{2\ell + 1} \cdot \left[ \left( \frac{\bar{\Omega}_{\text{GW}}}{4\pi} \right)^4 \text{Var}^4(\hat{A}_{\ell m}^e) + 2\text{Var}(\hat{A}_{\ell m}^e) \cdot \left( \frac{\bar{\Omega}_{\text{GW}}}{4\pi} \right)^2 C_{\ell}^{\text{GW}} \right], \qquad (\text{B21})$$

which can be directly computed from equation (B17). It follows that for a given pulsar number $N_p$ (determining the angular resolution of a time-residual map, as $L_{\text{max}} \sim \sqrt{N_p}$), equation (B21) provides us with the precision with which, in a perfect experiment, the angular power spectrum of the background energy density can be computed from a time-residual map.

We now consider the cross-correlation with the galaxy number counts, under the assumption that the galaxy map has an infinite angular resolution. We compute

$$C_{\ell}^{\text{GW}\Delta\ell} = \frac{1}{2\ell + 1} \sum_{m} a_{\ell m}^{e} b_{\ell m}^{*}, \qquad (\text{B22})$$

where $b_{\ell m}$ are the multipoles of the galaxy number counts, defined in equation (6b). Repeating similar steps as above we get

$$\text{Var}(C_{\ell}^{\text{GW}\Delta\ell}) = \frac{C_{\ell}^{\Delta}}{2\ell + 1} \left(\frac{\bar{\Omega}_{\text{GW}}}{4\pi}\right)^2 \text{Var}(\hat{A}_{\ell m}^{e}), \qquad (\text{B23})$$

which can be directly computed from equation (B17).

## APPENDIX C: SIGNAL TO NOISE DERIVATION

Let us assume that our observables are the angular power spectra of the AGWB and galaxy autocorrelations and their cross-correlation angular spectrum, i.e. $C_{\ell}^{\text{GW}}$, $C_{\ell}^{\Delta}$, and $C_{\ell}^{\text{GW}\Delta}$ are the observables. The estimator (full sky) of the angular power spectrum is

$$\hat{C}_{\ell}^{ab} = \frac{1}{2\ell + 1} \sum_{m} a_{\ell m}^{a} (a_{\ell m}^{b})^*, \qquad (\text{C1})$$

where $a_{\ell m}^{a,b}$ are Gaussian maps and the indices are defined as $a = (\text{GW}, \Delta)$ and $b = (\text{GW}, \Delta)$. We compute the covariance of this estimator. To make the notation compact, we introduce $X = (C_{\ell\min}^{\text{GW}}, C_{\ell\min}^{\Delta}, C_{\ell\min}^{\text{GW}\Delta}, \dots, C_{\ell\max}^{\text{GW}}, C_{\ell\max}^{\Delta}, C_{\ell\max}^{\text{GW}\Delta})$. Then the associated covariance matrix is block diagonal

$$C = \begin{pmatrix} B_{\ell\min} & 0 \\ \vdots & \vdots \\ 0 & B_{\ell\max} \end{pmatrix},$$

with 3×3 blocks

$$B_{\ell} = \frac{1}{2\ell + 1} \cdot \begin{pmatrix} 2(C_{\ell}^{\text{GW}})^2 & 2(C_{\ell}^{\text{GW}\Delta})^2 & 2(C_{\ell}^{\text{GW}})(C_{\ell}^{\text{GW}\Delta}) \\ 2(C_{\ell}^{\text{GW}\Delta})^2 & 2(C_{\ell}^{\Delta})^2 & 2(C_{\ell}^{\Delta})(C_{\ell}^{\text{GW}\Delta}) \\ 2(C_{\ell}^{\text{GW}})(C_{\ell}^{\text{GW}\Delta}) & 2(C_{\ell}^{\Delta})(C_{\ell}^{\text{GW}\Delta}) & C_{\ell}^{\text{GW}}C_{\ell}^{\Delta} + (C_{\ell}^{\text{GW}\Delta})^2 \end{pmatrix},$$

where it is understood that each entry includes a noise part, i.e. one has to replace $C_{\ell}^{ab} \rightarrow C_{\ell}^{ab} + N_{\ell}^{ab}$. We have all the ingredients to compute the (non-Gaussian) likelihood of the angular power spectra $X$ and compute the corresponding Fisher matrix deriving it with respect to the parameters $\theta_a$. The Fisher matrix is defined as

$$F_{ab} = \left\langle -\frac{\partial^2 \ln \mathcal{L}}{\partial \theta_a \partial \theta_b} \right\rangle, \qquad (\text{C2})$$

where $\theta_a$ are a set of parameters of the model. It is possible to show that only the Gaussian part of the likelihood enters the Fisher matrix. We want to estimate the signal to noise of the galaxy map, hence we parametrize (adding a generic noise component to the angular power spectrum of auto and cross-correlations)

$$X_{\ell} = (C^2 C_{\ell}^{\text{GW}} + N_{\ell}^{\text{GW}}, C_{\ell}^{\Delta} + N_{\ell}^{\Delta}, C C_{\ell}^{\text{GW}\Delta} + N_{\ell}^{\text{GW}\Delta}), \qquad (\text{C3})$$

where $C$ is a book keeping parameter for the amplitude. Then

$$(F_C)_\ell = \left(\frac{S}{N}\right)_\ell^2 = (\partial_c X_\ell) \cdot B_\ell^{-1} \cdot (\partial_c X_\ell), \qquad (\text{C4})$$

evaluated in $C = 1$. If we use only the information from the autocorrelation, we get equation (32) while using only cross-correlation we obtain equation (33).

## APPENDIX D: THE GW ENERGY DENSITY

Let us introduce some basic quantities that we extract from the catalogue. The chirp mass of a source $M$ is given by

$$M = \frac{(m_1 m_2)^{3/5}}{(m_1 + m_2)^{1/5}}, \qquad (\text{D1})$$

with $m_{1,2}$ the masses of the two black holes in the binary. We focus on the case of an inspiraling BBH, for which the emitted power is given by (Maggiore 2007)

$$\frac{\mathrm{d}E_s}{\mathrm{d}t_s} = \frac{32}{5Gc^5}(\pi GMf_s)^{10/3} = \frac{32}{5Gc^5}(\pi GM(1+z)f)^{10/3}, \qquad (\text{D2})$$

where $E_s$ is the source-frame energy of one source with mass $M$ and (source-frame) frequency $f_s$, and $t_s$ is the local cosmic time at the source. The source-frame frequency $f_s$ is related to the observed frequency by $f = f_s/(1+z)$. The observed power is found using $E_o = E_s/(1+z)$ and $dt_o = dt_s(1+z)$, hence

$$\frac{\mathrm{d}E_o}{\mathrm{d}t_o} = \frac{1}{(1+z)^2} \frac{\mathrm{d}E_s}{\mathrm{d}t_s}. \qquad (\text{D3})$$

We denote the (comoving) number of sources per comoving volume $dV_c$, per frequency $df$, and per logarithmic chirp mass $d\log M$ by

$$a^3 n_s(n, z, M, f) \equiv f \frac{\mathrm{d}N_s}{\mathrm{d}V_c \mathrm{d}\log(M) \mathrm{d}f} \propto f^{-8/3}, \qquad (\text{D4})$$

where $N_s$ is the number of sources. The previous scaling in frequency is found using that the number of sources in a given frequency band is proportional to the time spent in that band, along with $d\log f/dt \propto$
$f^{8/3}$. We split this quantity into an angle-averaged term $\bar{n}_s$ and an anisotropic contribution

$$n_s(n, z, \mathcal{M}, f) = \bar{n}_s(z, \mathcal{M}, f) + \delta n_s(n, z, \mathcal{M}, f). \quad (D5)$$

To consider the contributions of the different directions $n$ to the total energy density parameter at the observer, we then analogously split $\Omega_{GW}$ in an angle-averaged term $\bar{\Omega}_{GW}$ and an anisotropic fluctuation, with the definition (1). The isotropic part is obtained by summing all contributions as in equation (2a) with

$$\partial_r \bar{\Omega}_{GW}(f, r) = \int d\log \mathcal{M} \frac{\partial \bar{\Omega}_{GW}(f, r, \mathcal{M})}{\partial r \partial \log \mathcal{M}}, \quad (D6)$$

where

$$\frac{\partial \bar{\Omega}_{GW}(f, r, \mathcal{M})}{\partial r \partial \log \mathcal{M}} = \frac{1}{\rho_c (1+z)^2} \frac{1}{dt_s} a^3 \bar{n}_s(z, \mathcal{M}, f), \quad (D7)$$

and the function $z = z(r)$ must be used to relate a given redshift $z$ to the conformal distance $r$.

Denoting the fractional fluctuation in the number density of sources by

$$\delta_s(n, z, \mathcal{M}) = \frac{\delta n_s(n, z, \mathcal{M}, f)}{\bar{n}_s(z, \mathcal{M}, f)}, \quad (D8)$$

and using that the GW sources are biased tracers of the underlying distribution of galaxies,

$$\delta_s(n, z, \mathcal{M}) = b_s(z, \mathcal{M}) \delta_g(n, z), \quad (D9)$$

we obtain

$$\delta \Omega_{GW}(f, \mathbf{n}) = \frac{1}{4\pi} \int dr \int d\log \mathcal{M} \frac{\partial \bar{\Omega}_{GW}(f, r, \mathcal{M})}{\partial r \partial \log \mathcal{M}} b_s(z, \mathcal{M}) \delta_g(\mathbf{n}, z), \quad (D10)$$

where the integration on the conformal distance $r$ can be replaced by an integral on redshifts using the background cosmology and $dr = dz/H(z)$.

Equation (D10) relates the anisotropies in the energy density of an astrophysical SGWB to the anisotropies of the large-scale structure where SGWB sources are found. We assume that GW sources follow the same clustering structure as the galaxies, hence we take $b_s = 1$ in our analysis.

## APPENDIX E: ENERGY DENSITY AND SHOT NOISE FROM A CATALOGUE

For a given source $i$, we introduce the following quantity $P_i \equiv (dE_o/dt_o)_i / \rho_c = (dE_s/dt_s)_i / \rho_c / (1+z_i)^2$, such that for a collection of $N$ sources in directions $\mathbf{n}_i$

$$\Omega_{GW}(\mathbf{n}) = \sum_{i=1}^{N} \left( \frac{P_i}{4\pi r_i^2} \right) \delta^2(\mathbf{n}, \mathbf{n}_i). \quad (E1)$$

Using

$$\delta^2(\mathbf{n}_1, \mathbf{n}_2) = \sum_{\ell m} Y_{\ell m}(\mathbf{n}_1) Y_{\ell m}^*(\mathbf{n}_2), \quad (E2)$$

and an expansion $\Omega_{GW}(\mathbf{n}) = \sum_{\ell m} \Omega_{\ell m}^{GW} Y_{\ell m}(\mathbf{n})$ we get

$$\Omega_{\ell m}^{GW} = \sum_{\ell} \left( \frac{P_{\ell}}{4\pi r_{\ell}^2} \right) Y_{\ell m}^{*}(\mathbf{n}_{\ell}). \quad (E3)$$

An estimator of the $\hat{\mathcal{C}}_{\ell}^{GW}$ defined in equation (5a) is

$$\hat{\mathcal{C}}_{\ell}^{GW} = \frac{1}{2\ell + 1} \sum_{m=-\ell}^{\ell} \Omega_{\ell m}^{GW} \Omega_{\ell m}^{GW*}, \quad (E4)$$

and for random positions of the sources, that is sources which are not located according to clustering properties, it is non-vanishing only due to the Poissonian nature of the source distribution. We want to average this estimator over the direction realizations of the sources, that is over random realizations of the set of directions $\{\mathbf{n}_i\}$, so as to estimate the shot noise. We shall use that [see e.g. Allen et al. (2024) or appendix D of Pitrou & Cusin (2025)]

$$\langle Y_{\ell m}(\mathbf{n}_{i1}) \rangle_{\{\mathbf{n}_i\}} = \frac{\delta_{\ell}^{0} \delta_{m}^{0}}{\sqrt{4\pi}}, \quad (E5)$$

and for $(\ell, m) \neq (0, 0)$

$$\langle Y_{\ell m}(\mathbf{n}_{i1}) Y_{\ell' m'}^{*}(\mathbf{n}_{i1}) \rangle_{\{\mathbf{n}_i\}} = \frac{\delta_{\ell \ell'} \delta_{m m'}}{4\pi}. \quad (E6)$$

Then when averaging over all random positions of these $N$ sources

$$4\pi \langle \Omega_{GW}(\mathbf{n}) \rangle_{\{\mathbf{n}_i\}} = \sqrt{4\pi} \langle \Omega_{00}^{GW} \rangle_{\{\mathbf{n}_i\}} = \sum_{i=1}^{N} \frac{P_i}{4\pi r_i^2}, \quad (E7)$$

$$\langle \hat{\mathcal{C}}_{\ell}^{GW} \rangle_{\{\mathbf{n}_i\}} = \frac{1}{4\pi} \sum_{i=1}^{N} \left( \frac{P_i}{4\pi r_i^2} \right)^2 . \quad (E8)$$

We identify the average GW density as $\bar{\Omega}_{GW} = 4\pi \langle \Omega_{GW}(\mathbf{n}) \rangle_{\{\mathbf{n}_i\}}$ and shot noise with $N_{\text{shot}}^{GW} = \langle \hat{\mathcal{C}}_{\ell}^{GW} \rangle_{\{\mathbf{n}_i\}}$.

In order to check the consistency with the continuous integral formulas we use

$$\sum_{i=1}^{N} \rightarrow \int d\log f \int d\log \mathcal{M} \int dr (4\pi r^2 a^3 \bar{n}_s), \quad (E9)$$

and indeed applying this replacement rule to equation (E7) we get equation (D7). Also for shot noise we get

$$\begin{aligned} \langle \hat{\mathcal{C}}_{\ell}^{GW} \rangle_{\{\mathbf{n}_i\}} &\approx \int d\log f \int d\log \mathcal{M} \int dr a^3 \bar{n}_s r^2 \left( \frac{P}{4\pi r^2} \right)^2 \\ &\approx \int d\log f \int d\log \mathcal{M} \int dr \frac{a^3 \bar{n}_s}{(4\pi)^2 r^2} \left( \frac{dE_s}{dt_s} \right)^2 \end{aligned} \quad (E10)$$

and using equation (D7) we recover equation (26). Hence when using a catalogue, we have to use equation (E7) to estimate the background GW density, and equation (E8) to estimate the shot noise.

This paper has been typeset from a TeX/LA TeX file prepared by the author.