# The impact of large-scale galaxy clustering on the variance of the Hellings-Downs correlation: numerical results

Nastassia Grimm, Martin Pijnenburg, Giulia Cusin, Camille Bonvin

► To cite this version:

Nastassia Grimm, Martin Pijnenburg, Giulia Cusin, Camille Bonvin. The impact of large-scale galaxy clustering on the variance of the Hellings-Downs correlation: numerical results. Journal of Cosmology and Astroparticle Physics, 2025, 04, pp.047. 10.1088/1475-7516/2025/04/047. hal-04805779

HAL Id: hal-04805779

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

Submitted on 2 Nov 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
PAPER • OPEN ACCESS

# The impact of large-scale galaxy clustering on the variance of the Hellings-Downs correlation: numerical results

To cite this article: Nastassia Grimm et al JCAP04(2025)047

View the [article online](#) for updates and enhancements.

You may also like

- [Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I](#)

Heng Xu, Siyuan Chen, Yanjun Guo et al.

- [The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background](#)

Zaven Arzoumanian, Paul T. Baker, Harsha Blumer et al.

- [Cosmological backgrounds of gravitational waves](#)

Chiara Caprini and Daniel G Figueroa
RECEIVED: November 14, 2024

REVISED: February 6, 2025

ACCEPTED: March 15, 2025

PUBLISHED: April 16, 2025

# The impact of large-scale galaxy clustering on the variance of the Hellings-Downs correlation: numerical results

Nastassia Grimm<sup>D,a</sup> Martin Pijnenburg<sup>D,a</sup> Giulia Cusin<sup>D,b,a</sup> and Camille Bonvin<sup>D,a</sup>

E-mail: nastassia.grimm@unige.ch, martin.pijnenburg@unige.ch, giulia.cusin@iap.fr, camille.bonvin@unige.ch

**ABSTRACT:** Pulsar timing array experiments have recently found evidence for a stochastic gravitational wave (GW) background, which induces correlations among pulsar timing residuals described by the Hellings and Downs (HD) curve. Standard calculations of the HD correlation and its variance assume an isotropic background. However, for a background of astrophysical origin, we expect a higher GW spectral density in directions with higher galaxy number densities. In a companion paper, we have developed a theoretical formalism to account for the anisotropies arising from large-scale galaxy clustering, leading to a new contribution to the variance of the HD correlation. In this subsequent work, we provide numerical results for this novel effect. We consider a GW background resulting from mergers of supermassive black hole binaries, and relate the merger number density to the overdensity of galaxies. We find that anisotropies due to large-scale galaxy clustering lead to a standard deviation of the HD correlation at most at percent level, remaining well below the standard contributions to the HD variance. Hence, this kind of anisotropies in the GW source distribution does not represent a substantial contamination to the correlations of timing residuals in present and future PTA surveys. Suitable statistical methods to extract the galaxy clustering signal from PTA data will be investigated in the future.

**KEYWORDS:** galaxy clustering, gravitational waves / sources, gravitational waves / theory

ARXIV EPRINT: 2411.08744

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

<footnote><sup>a</sup>Département de Physique Théorique and Center for Astroparticle Physics, Université de Genève, Quai E. Ansermet 24, CH-1211 Geneve 4, Switzerland</footnote>

<footnote><sup>b</sup>Institut d'Astrophysique de Paris, UMR-7095 du CNRS, 98bis boulevard Arago, F-75014 Paris, France</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Basic concepts related to the Hellings and Downs correlation</td><td>3</td></tr><tr><td>2.1 Description of the signal</td><td>3</td></tr><tr><td>2.2 The Hellings and Downs correlation in the clustered Universe</td><td>4</td></tr><tr><td>2.3 Different types of variance</td><td>5</td></tr><tr><td>2.4 Galaxy clustering variance</td><td>6</td></tr><tr><td>3 Astrophysical modeling</td><td>8</td></tr><tr><td>3.1 Modeling of bGW</td><td>8</td></tr><tr><td>3.2 Reconstruction of bGW from catalogs</td><td>11</td></tr><tr><td>3.3 Discussion on the source bias and galaxy bias</td><td>12</td></tr><tr><td>4 Numerical results for the clustering variance</td><td>13</td></tr><tr><td>4.1 Methods</td><td>13</td></tr><tr><td>4.2 Results</td><td>16</td></tr><tr><td>4.3 Comparison to Allen et al. (2024)</td><td>19</td></tr><tr><td>5 Conclusion</td><td>20</td></tr><tr><td>A Pulsar-averaged galaxy clustering variance</td><td>21</td></tr><tr><td>B Use of the astrophysical catalogs</td><td>24</td></tr></table>

## 1 Introduction

Pulsar timing arrays (PTAs) were used to deliver the first evidence of a stochastic gravitational wave (GW) background in the nHz band [1–4]. The correlation of pulsar timing residuals due to the passing of GWs, as a function of the pulsar pair separation on the sky, is usually assumed to be described by the Hellings and Downs (HD) curve [5].<fnref n="1" /> However, it has been recently shown that the stochasticity of GW sources, and the fact that the number of pulsar pairs at a given separation is finite, cause a departure from this idealized curve [6–11]. In these studies, however, the GW sources are assumed to be isotropically distributed. In ref. [12] we presented a novel framework to account for the impact of the cosmological large-scale structure on the HD correlation.<fnref n="2" /> Indeed, since galaxies are not isotropically distributed on the sky but cluster to form the large-scale structure, we expect an intrinsic isotropy to be present as well for astrophysical GW sources.

<footnote n="1">¹We note that we use the notion *HD curve* to refer to the theoretical prediction by Hellings and Downs [5], while *HD correlation* corresponds to the observed correlation that is in fact measured in our Universe and might differ from the idealized HD curve.</footnote>

<footnote n="2">²Prior to our work in ref. [12], the possibility of anisotropies in the distribution of GW sources in the PTA band has been pointed out in various studies, however either not in the context of the HD correlation [13–23], or without targeting the impact on its variance [24–27]. The detectability of kinematic anisotropies in the SGWB with different PTA experiments has been recently studied in ref. [28].</footnote>
The timing residual measured for a given pulsar is linear in the total GW strain in a given direction, due to a superposition of GWs, emitting incoherently by a large number of sources. Such a signal is stochastic for two distinct reasons: first because we cannot predict the properties of the emitted waves (in particular the phases), and second because we cannot predict the distribution of sources. Usually, the second source of stochasticity is ignored: sources are assumed to be isotropically distributed, and there is therefore only one possible realization of the distribution of sources. The HD curve is then defined as the product of the timing residuals of two pulsars, averaged over all realizations of the GW emission.

As we have shown in ref. [12], to account for anisotropies in the distribution of sources due to large-scale galaxy clustering, one needs to add an additional average over all possible realizations of the distribution of sources. We considered an ensemble of GW sources in one realization of the universe (i.e. a fixed distribution of sources), and then analyzed how the HD correlation varies from realization to realization. The mean of the HD correlation itself is not affected by clustering, which means that anisotropies in the distribution of sources do not generate a systematic bias. However, the large-scale structure of the Universe induces a variance in the measurement of the HD correlation, quantifying to which amount the measurements in our specific realization of the Universe can differ from the mean.<fnref n="3" />

In this work, we numerically evaluate the size of this variance. To do so, we assume that the observed background is given by the GW emission of supermassive black hole (SMBH) binaries in the inspiraling phase. For our case study we choose to work with the catalogs of refs. [29, 30], and we show results for the three population models considered therein: the popIII model based on light growth seeds for SMBHs; and the other two heavy-seed Q3 SMBH models, with (Q3-d) and without (Q3-nd) delays between galaxy merger and SMBH binary merger. We assume that the overdensity of black hole mergers is a biased tracer of the underlying galaxy overdensity, which itself is a biased tracer of the matter distribution. We find that anisotropies due to large-scale galaxy clustering lead to a departure from the HD curve at a sub-percent level, ranging from 0.54% to 2.85% depending on pulsar pair separation, i.e. well below the precision of current PTA surveys [1-4]. Moreover, it is well below the standard variance calculated in ref. [7] under the assumption of an isotropic source distribution. This indicates that the clustering of sources can safely be neglected in current and future analyses.

In addition to these results, which refer to a single pulsar pair, we also compute the variance related to an (infinite) continuous distribution of pulsar pairs, which we refer to as the *pulsar-averaged* case.<fnref n="4" /> Averaging over pairs of pulsars reduces the clustering variance by an order of magnitude, since anisotropies in the distribution of sources affect pulsar pairs

<footnote n="3">³In this work and in the companion paper [12] we consider only the effect of clustering (i.e. fluctuations in the distribution of sources) on the variance of the HD correlation, neglecting relativistic effects (due to gravitational potentials and source velocities) that are expected to be subdominant on all scales, see e.g. the calculation in refs. [14, 20] in the context of ground-based detectors.</footnote>

<footnote n="4">⁴In a previous version of this manuscript, we have referred to the pulsar-averaged as the irreducible case. However, using the term *irreducible* in this context is misleading. As recently shown in ref. [31] (and further supported by ref. [32] for a real-space result), cosmic variance can be mitigated through appropriate frequency weighting. We likewise expect that the clustering-induced variance can also be reduced using a similar weighting approach.</footnote>
at different sky locations differently and are therefore washed out by the average over an infinite number of pulsars. The pulsar-averaged and single pair cases provide, respectively, a lower and upper bound on the clustering variance in a real observation with a finite number of pulsar pairs.

We compare our results with the approximate computations of the clustering variance derived in ref. [33] for the pulsar-averaged case. We find that the pulsar-averaged clustering variance is, with a relative amplitude of 0.10 – 0.14%, an order of magnitude smaller than the upper bound of 1% claimed in ref. [33]. This difference is due to the different modeling of the clustering correlation of GW sources. Ref. [33] assumes a constant value for the angular power spectrum of galaxy correlations (independent on angular separation) and does not take any redshift dependence into account. In contrast, we are using the clustering correlation predicted in a ΛCDM cosmological model, which properly accounts for the non-trivial dependence of clustering correlations on source angular separation and on redshift. Our formalism provides therefore a clear link between our cosmological and astrophysical modeling and the variance of the HD correlation.

The remaining part of this paper is structured as follows. In section 2, we describe basic concepts concerning PTAs and galaxy clustering in the context of this work. We also give an overview over different types of variance contributions to the HD correlation introduced in this and previous works [7, 12], and summarize the analytical results of ref. [12] for the clustering variance. Then, in section 3, we present the modeling and catalog reconstruction of the quantity $b_{GW}$, which relates galaxy density to GW spectral density and modulates the amplitude of the clustering variance. The fact that SMBH mergers are biased tracers of the underlying galaxy and matter fields are discussed as well. Section 4 presents the numerical results for the clustering variance obtained in this work. We also describe our numerical methods and provide a comparison to ref. [33] in this section. Finally, we conclude in section 5. Moreover, in appendix A, we present the analytical results for the pulsar-averaged clustering variance in the case of infinitely many pulsar pairs.

## 2 Basic concepts related to the Hellings and Downs correlation

The theoretical modeling of the HD correlation and its variance, accounting for anisotropies in the distribution of sources, has been developed in ref. [12]. Here, we summarize the key concepts and notations.

### 2.1 Description of the signal

The passing of GWs induces redshifts in pulsar observations, i.e. small shifts in the arrival period between two pulses.<fnref n="5" /> We write the GW propagating in direction $\mathbf{n}$ (due to the sum of waves from various sources propagating in that direction) as

$$h_{ij}(t, \mathbf{n}) \equiv \int df \sum_{A=+, \times} e_{ij}^{A}(\mathbf{n}) h_{A}(f, \mathbf{n}) e^{2\pi i ft}, \qquad (2.1)$$

<footnote n="5">⁵To avoid confusion, we denote with a capital Z the redshift of pulsars induced by the passing of GWs, and with z the cosmological redshift (of the sources) due to the Universe’s expansion. When not explicitly specified, the term “redshift” refers to the latter.</footnote>
where $e_{ij}^A(\mathbf{n})$ is the wave polarization tensor. The pulsar redshift $Z_a$, for a pulsar located in direction $\mathbf{n}_a$, is given by (see e.g. ref. [7]),

$$Z_a(t) = \sum_A \int df \int d\mathbf{n} F_a^A(\mathbf{n}) h_A(f, \mathbf{n}) e^{2\pi i f t} (1 - e^{-2\pi i f \tau_a (1+\mathbf{n}\cdot\mathbf{n}_a)}), \quad (2.2)$$

where $\tau_a$ is the light travel time between emission and observation. Furthermore, the antenna pattern functions $F_a^A(\mathbf{n})$ are defined as

$$F_a^A(\mathbf{n}) = \frac{n_a^i n_a^j e_{ij}^A(\mathbf{n})}{2(1 + \mathbf{n} \cdot \mathbf{n}_a)}. \quad (2.3)$$

The pulsar redshift is directly proportional to the amplitude (or strain) of the GWs, $h_A(f, \mathbf{n})$, where $A$ denotes the polarization of the wave, + and $\times$.

The aim of PTA observations is to detect correlated shifts in the pulses of two distinct pulsars, due to the passing of GWs. We therefore introduce the quantity $\rho_{ab}$, defined as the product of $Z_a$ and $Z_b$ for two pulsars located at $\mathbf{n}_a$ and $\mathbf{n}_b$, averaged over the observation time $T$,

$$\rho_{ab} \equiv \overline{Z_a Z_b}. \quad (2.4)$$

The HD curve describes the theoretical expectation of this quantity, and it is obtained by taking the ensemble average of eq. (2.4) accounting for both the stochasticity of GWs and the fact that sources are distributed in a stochastic way. More details on these different layers of stochasticity can be found in ref. [12]. Therein, we have used distinct notations to highlight the difference between full ensemble averages and ensemble averages assuming a given density distribution. In this work, as we focus on the numerical predictions, we omit this distinction in our notation and refer to our previous work for the formal details.

## 2.2 The Hellings and Downs correlation in the clustered Universe

PTA surveys measure the angular correlations of pulsar redshifts, with theoretical expectation value given by $\langle \rho_{ab} \rangle$. Here, the average is in fact a double average: a first average over GW realizations and a second average over stochastic initial conditions of the matter overdensity (of which the distribution of GW sources is a biased tracer). From eqs. (2.2) and (2.4), we see that the quantity $\langle \rho_{ab} \rangle$ contains ensemble averages over the product of two $h_A(f, \mathbf{n})$. For a given realization of the density field, this ensemble average over GW realizations is given by

$$\langle h_A^*(f, \mathbf{n}) h_B(f', \mathbf{n}') \rangle = \frac{1}{2} \delta_{AB} \delta(f - f') S_h(f, \mathbf{n}) \frac{\delta^2(\mathbf{n}, \mathbf{n}')}{4\pi}, \quad (2.5)$$

where the delta functions ensure that only GWs with the same polarization state, at the same frequency and from the same direction are correlated.<fnref n="6" /> The spectral density $S_h(f, \mathbf{n})$ is usually assumed to be isotropic, $S_h(f, \mathbf{n}) = \bar{S}_h(f)$. Then, the HD correlation is determined by

$$\langle \rho_{ab} \rangle = h^2 \mu_{\text{HD}}(\gamma), \quad h^2 = \int df \bar{S}_h(f), \quad (2.6)$$

<footnote n="6">⁶This is another way of stating that, when we sum over a large number of products of incoherent signals, non-vanishing contributions arise only from products of signals with equal phases, see e.g. ref. [14].</footnote>
where $h^2$ is the squared strain, $\gamma$ denotes the angle between the pulsar directions $\mathbf{n}_a$ and $\mathbf{n}_b$, and the HD curve $\mu_{\text{HD}}(\gamma)$ is given by

$$ \mu_{\text{HD}}(\gamma) = \int \frac{d\mathbf{n}}{4\pi} \chi_{ab}(\mathbf{n}) = \frac{1}{4} + \frac{1}{12}\cos(\gamma) + \frac{1}{2}(1 - \cos(\gamma))\log\left(\frac{1 - \cos(\gamma)}{2}\right). \quad (2.7) $$

Indeed, due to the statistical isotropy of $h_A$, the HD curve depends on $\mathbf{n}_a$ and $\mathbf{n}_b$ only via their angular separation $\gamma$. More specifically, it is determined via the sky-average of the geometrical objects $\chi_{ab}(\mathbf{n})$,

$$ \chi_{ab}(\mathbf{n}) = \sum_{A=+,\times} F_a^A(\mathbf{n}) F_b^A(\mathbf{n}). \quad (2.8) $$

In ref. [12] and in this work, we account for the fact that galaxies follow the cosmological large-scale structure. This leads to higher galaxy densities and therefore, in case of an astrophysical origin of the stochastic GW background, to higher spectral densities in certain directions. Splitting $S_h(f, \mathbf{n})$ into an isotropic and anisotropic part, we can relate it to the fractional galaxy overdensity $\delta_g(\mathbf{n}, z)$ as

$$ S_h(f, \mathbf{n}) = \bar{S}_h(f) + \delta S_h(f, \mathbf{n}) = \bar{S}_h(f) \left[ 1 + \int dz \, b_{\text{GW}}(z) \delta_g(\mathbf{n}, z) \right], \quad (2.9) $$

where $b_{\text{GW}}(z)$ is some redshift dependent factor$^7$ that relates GW anisotropies to galaxy anisotropies. It depends on the details of the astrophysical model in consideration and we will compute it below for the case where the SGWB is due to SMBH mergers. We note that $b_{\text{GW}}(z)$ does not depend on the direction $\mathbf{n}$: since we assume that the distribution of GW sources follows the distribution of galaxies, the anisotropic behavior of $S_h(f, \mathbf{n})$ is determined by $\delta_g(\mathbf{n}, z)$ alone.

To obtain the expectation value of the HD correlation while taking galaxy clustering into account, we need to take the ensemble average of eq. (2.5) over all realizations of the galaxy density field $\delta_g$. In ref. [12], we showed that since $\rho_{ab}$ is linear is $\delta_g$ and since this quantity has zero mean, the anisotropic part $\delta S_h(f, \mathbf{n})$ vanishes when averaged over all realizations. The mean of the HD correlation is therefore unaffected by anisotropies in the distribution of sources. This is, however, not true for the variance of the HD correlation, which obtains an additional contribution due to galaxy clustering. In the following, we will summarize different types of variances, before discussing analytical and numerical results in the subsequent sections.

## 2.3 Different types of variance

The aim of our work is to provide a complete computation of the variance of the HD correlation, i.e. to assess by how much a measurement of angular correlations of pulsar redshifts in a single realization of the universe can differ from the expectation value $\langle \rho_{ab} \rangle$, taking all sources of variance in the presence of anisotropies into account. We distinguish different types of variances, due to the different sources of stochasticity in the signal.

<footnote>$^7$We note that in ref. [12], to be as general as possible, we had introduced $b_{\text{GW}}(f, z)$ as a redshift and frequency dependent factor. In this work, however, we will see that $b_{\text{GW}}(z)$ appears as a frequency independent factor in our catalog-based astrophysical modeling.</footnote>
The standard variance of the HD correlation computed in refs. [6–10] is evaluated assuming a fixed, isotropic realization of the distribution of sources, i.e. neglecting the stochastic nature of the source distribution. It is related to the 4-point function of the strain averaged over all GW realizations. We refer to this variance as *standard variance*. It describes how the signal can differ from the HD curve due to the fact that we observe only one realization of the GW strain. In contrast, in ref. [12], we also accounted for the stochasticity in the source distribution. We refer to the resulting variance as *clustering variance*. It describes how the signal can differ from the HD curve due to the fact that we observe only one realization of the density field.

In eq. (2.4), the signal is computed for a single pair of pulsars *a* and *b*. In practice however, one can average over all pairs of pulsars in the array that share the same angular separation $\gamma$. This reduces the variance, since anisotropies in the GW signal affect different pulsar pairs (in different directions) differently. Averaging over many pairs therefore erases part of the anisotropies. The reduction of the variance depends on the number of pulsar pairs per angular separation $\gamma$.<fnref n="8" /> The limiting case, corresponding to a continuous (infinite) distribution of pulsars, provides a lower bound to the variance. We call these variances *pulsar-averaged standard variance* and *pulsar-averaged clustering variance*. A summary of the different contributions to the variance of the HD correlation is presented in table 1. In the main part of this work, we show the computation of the clustering variance in the case of a single pulsar pair per angular separation $\gamma$ (which provides an upper bound). Results for the pulsar-averaged variance, obtained by averaging over a continuous, infinite number of pulsars, will be derived and presented in appendix A.

Finally, we note that in principle, an additional shot noise contribution arises from the fact that the signal is due to a discrete number of GW sources (see, e.g., ref. [21]). This contribution accounts for the fact that even if the distribution of matter would be homogeneously distributed, the number of discrete sources tracing this homogeneous background is not perfectly homogeneous, but it follows a Poisson distribution around the mean. This also impacts the PTA signal, see ref. [33], but is not subject of this work. Here, we focus on the galaxy clustering contribution, which constitutes an intrinsic source of anisotropy independent of shot noise.

## 2.4 Galaxy clustering variance

In ref. [12], we have derived an analytical expression for the contribution of large-scale galaxy clustering on the variance of the HD correlation, denoted by $\langle \rho_{ab}^2 \rangle^{\text{clust}}$.<fnref n="9" /> This calculation is based on a two-step averaging procedure where we first assume a fixed galaxy density distribution and then take ensemble averages over all such distributions. The resulting

<footnote n="8">⁸Investigations of the standard variance for various finite sets of pulsar pairs have been performed in ref. [8]. In particular, figure 9 therein shows that, for 88 pulsars corresponding to the specifications of the International Pulsar Timing Array [1–4], the variance is already considerably reduced compared to the case of only one pair per angular separation. Future data sets, containing an even larger number of pulsars, will approach the pulsar-averaged limit (with an infinite set of pulsars) even more. We stress again that this limit can be further reduced with a proper frequency weighting, as recently shown in ref. [31].</footnote>

<footnote n="9">⁹Note again that there is no contribution of galaxy clustering to the mean of the HD correlation, which can be expressed as $\langle \rho_{ab} \rangle^{\text{clust.}} = 0$.</footnote>
<table>
  <tr>
    <td>Standard variance</td>
    <td>Variance due to the stochasticity of the GW emission, assuming an isotropic distribution of sources and one pulsar pair</td>
    <td>Eq. (C28) in [7]</td>
  </tr>
  <tr>
    <td>Clustering variance</td>
    <td>Variance due to the stochasticity in the anisotropic distribution of GW sources induced by galaxy clustering, for one pulsar pair</td>
    <td>Eq. (2.10),<br>Eq. (31) in [12]</td>
  </tr>
  <tr>
    <td>Pulsar-averaged standard variance (cosmic variance)</td>
    <td>Standard variance in the continuous pulsar pair limit</td>
    <td>Eq. (A.4),<br>Eq. (C45) in [7]</td>
  </tr>
  <tr>
    <td>Pulsar-averaged clustering variance</td>
    <td>Clustering variance in the continuous pulsar pair limit</td>
    <td>Eq. (A.9)</td>
  </tr>
</table>

**Table 1.** We list the different contributions to the variance of the HD correlation, including a description and their defining equation. The first two lines correspond to the case where the signal arises from only one pulsar pair. The remaining lines correspond to the observational limit of a continuous sample of pulsars on the sky.

variance is given by

$$ \langle \rho_{ab}^2 \rangle^{\text{clust.}} = \iint dz \, dz' [G(z)G(z') + \Gamma(z, z')] \iint \frac{dn}{4\pi} \frac{dn'}{4\pi} \xi_g(\mathbf{n} \cdot \mathbf{n}', z, z') \chi_{ab}(\mathbf{n})\chi_{ab}(\mathbf{n}') \\ + 4 \iint dz \, dz' \Gamma(z, z') \iint \frac{dn}{4\pi} \frac{dn'}{4\pi} \xi_g(\mathbf{n} \cdot \mathbf{n}', z, z') \chi_{aa}(\mathbf{n})\chi_{bb}(\mathbf{n}'). \quad (2.10) $$

Here, $\xi_g(\mathbf{n} \cdot \mathbf{n}', z, z')$ denotes the galaxy correlation function

$$ \xi_g(\mathbf{n} \cdot \mathbf{n}', z, z') \equiv \langle \delta_g(\mathbf{n}, z) \delta_g(\mathbf{n}', z') \rangle. \quad (2.11) $$

It encodes the correlation between overdensities of galaxies in two different directions **n** and **n'** of the sky (depending only on the angle between them due to statistical isotropy) and at redshifts *z* and *z'*. The form and amplitude of this correlation function directly governs the clustering variance. The functions *G(z)* and *Γ(z, z')* in eq. (2.10) are measures of the GW strain weighted by $b_{GW}(z)$,

$$ G(z) = \int_{0}^{\infty} df \bar{S}_{h}(f)b_{GW}(z), \\ \Gamma(z, z') = \frac{1}{4} \int_{0}^{\infty} \int_{0}^{\infty} df \, df' \, \text{sinc}^2(\pi(f-f')T) \bar{S}_{h}(f)\bar{S}_{h}(f')b_{GW}(z)b_{GW}(z'), \quad (2.12) $$

where sinc(*x*) = sin(*x*)/*x*. We note that the frequency range of integration is theoretically infinite for a detector sensitive to arbitrary frequencies. In practice, however, the frequency bounds of integration should be chosen to correspond to the frequency band of the experiment considered. For our numerical evaluations, we choose a lower bound $f_{min} = 2.1$ nHz (corresponding to $1/T$ for an observational time $T = 15$ yr, as in the latest NanoGrav data release [1]) and an upper bound $f_{max} = 100$ nHz. We note that the dependence of our results on the precise upper bound is small. Moreover, we note that in a realistic scenario, a sensitivity function should as well be included. Here, we aim to provide a theoretical prediction independently of the observational method, and thus omit any detector-specific sensitivity function.<fnref n="10" />

<footnote n="10">¹⁰We see from eqs. (2.6) and (2.12) that $G(z)$ and the squared strain $h^2$ have in fact the same frequency dependence. Thus, for the terms proportional to $G(z)$, the relative size of the variance compared to the signal will not be affected by a sensitivity function or a different choice of frequency bounds. The terms proportional to $\Gamma(z, z')$ are in principle affected by these details, although we do not expect them to significantly alter our results given their overall small amplitude (see section 4.2).</footnote>
Eq. (2.10) constitutes the main result of our previous work, ref. [12]. It corresponds to the case of observing a single pulsar pair. In appendix A, we show the analogous derivations for the pulsar-averaged variances, assuming a continuous distribution of infinitely many pulsars on the sky. Note that similar derivations have recently been presented in refs. [34, 35]. To move from this theoretical framework to numerical evaluations, we need to specify an astrophysically realistic model for the function $b_{GW}(z)$. This is done in the following section, applying three different models of SMBH mergers.

# 3 Astrophysical modeling

The clustering variance depends on the mass and redshift distribution of sources and on the astrophysical kernel $b_{GW}(z)$, defined through eq. (2.9). Moreover, the mean spectral density $\bar{S}_h$ affects the normalization of the HD correlation and its variance. We now model these quantities considering different astrophysical models for SMBH formation and evolution.

## 3.1 Modeling of $b_{GW}$

In order to derive an expression for the astrophysical kernel $b_{GW}(z)$ we follow different steps: first, we derive an expression for the background energy density, in terms of source distribution and emission rate. Then we relate energy density to spectral density $S_h$, and we determine the kernel by comparing with the parametrization in eq. (2.9).

Any GW carries a certain energy flux $\Phi$ (flow of energy per unit time per unit area), and the local flux can be related to the gauge invariant gravitational energy density $\rho_{GW}$ by [36]

$$ \rho_{GW} = c^{-1}\Phi. \qquad (3.1) $$

The flux associated with an astrophysical SGWB is the cumulative flux of all astrophysical sources. Thus, to each comoving volume element $dV_c$ located at comoving distance $r_z$, we may assign a certain GW luminosity $L_{GW}$, corresponding to the GW energy released per unit source time within this volume, see refs. [14, 37, 38] for more details. The energy density (or flux) is then the integral over comoving distance $r_z$ of the luminosity $L_{GW}$, with an inverse square luminosity distance factor $d_L^{-2}$ accounting for the fact that source frequencies and the energies are redshifted and spread out onto a sphere,<fnref n="11" /> that is:

$$ \rho_{GW}(\mathbf{n}) = c^{-1} \int dr_z r_z^2 \frac{1}{4\pi d_L^2} \frac{dL_{GW}}{dV_c} = \frac{1}{4\pi} \int dz \frac{1}{(1+z)^2 H(z)} \frac{dL_{GW}}{dV_c}(\mathbf{n}, z), \quad (3.2) $$

where we introduced the Hubble parameter $H(z)$ by changing variable from $r_z$ to $z$. Since GW sources are not isotropically distributed on the sky, the luminosity $L_{GW}$ and thus the energy density $\rho_{GW}$ are functions of the sky direction $\mathbf{n}$.

<footnote n="11">¹¹Note that the redshift $z$ and the luminosity distance $d_L$ in eq. (3.2) also depend on $\mathbf{n}$, since in an anisotropic universe, these quantities are not only due to the expansion of the Universe, but also to Doppler effects, gravitational redshift and integrated Sachs-Wolfe, see e.g. ref. [39]. The anisotropy due to the source distribution dominates however over the redshift and luminosity distance anisotropies and we neglect the latter here.</footnote>
The luminosity per comoving volume element in eq. (3.2) can be expressed as the integral over all masses and frequencies of the individual power $dE_s/dt_s$ radiated by each source within this volume element

$$ \frac{dL_{\text{GW}}}{dV_c} = \iint d\log(\mathcal{M}) df_s \frac{dN_s}{dV_c d\mathcal{M} d\mathcal{f}_s}(\mathbf{n}, z, \mathcal{M}, f_s) \frac{dE_s}{dt_s}(\log(\mathcal{M}), f_s) \quad (3.3) $$

$$ = \iint d\log(\mathcal{M}) df \left. \left( \frac{dN_s}{dV_c d\log(\mathcal{M}) df} \frac{dE_s}{dt_s} \right) \right|_{(\mathbf{n}, z, \mathcal{M}, f_s = f(1+z))}, \quad (3.4) $$

where $\mathcal{M}$ is the chirp mass of the source given by

$$ \mathcal{M} = \frac{(m_1 m_2)^{3/5}}{(m_1 + m_2)^{1/5}}, \quad (3.5) $$

with $m_{1,2}$ the masses of the two BHs in the binary. $N_s$ denotes the number of sources, $E_s$ is the source-frame energy of one source with mass $\mathcal{M}$ and (source-frame) frequency $f_s$, and $t_s$ is the local cosmic time at the source. In eq. (3.4), we changed variable from the source-frame frequency $f_s$ to the observed frequency $f = f_s/(1+z)$.

Denoting the number of sources per comoving volume $dV_c$, per frequency $df$ and per logarithmic chirp mass $d\log \mathcal{M}$ by

$$ n_s(\mathbf{n}, z, \mathcal{M}, f) \equiv \frac{dN_s}{dV_c d\log(\mathcal{M}) df}, \quad (3.6) $$

we then split it into an angle-averaged term $\bar{n}_s$ and an anisotropic contribution

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

We introduce the energy density parameter $\Omega_{\text{GW}}$,

$$ \Omega_{\text{GW}} \equiv \frac{1}{\rho_c} \frac{d\rho_{\text{GW}}}{d\log f}, \quad (3.8) $$

as the ratio of the energy density of the astrophysical SGWB per logarithmic observed frequency, $d\rho_{\text{GW}}/d\log f$, to the critical energy density $\rho_c$ of the Universe. To consider the contributions of the different directions $\mathbf{n}$ to the total energy density parameter at the observer, we then analogously split $\Omega_{\text{GW}}$ in an angle-averaged term $\bar{\Omega}_{\text{GW}}$ and an anisotropic fluctuation, as

$$ \Omega_{\text{GW}}(f, \mathbf{n}) = \bar{\Omega}_{\text{GW}}(f) + \delta\Omega_{\text{GW}}(f, \mathbf{n}). \quad (3.9) $$

We now focus on the clustering contribution $\delta n_s$ to the anisotropies in the distribution of sources in order to extract the quantity $b_{\text{GW}}$. Denoting the fractional fluctuation in the number density of sources by

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

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

$$ \delta_s(\mathbf{n}, z, \mathcal{M}) = b_s(z, \mathcal{M}) \delta_g(\mathbf{n}, z), \quad (3.11) $$
we obtain

$$ \delta\Omega_{\text{GW}}(f, \mathbf{n}) = \frac{f}{4\pi\rho_c} \iint dz d\log(\mathcal{M}) \frac{1}{(1+z)^2 H(z)} \frac{dE_s}{dt_s} \bar{n}_s(z, \mathcal{M}, f) b_s(z, \mathcal{M}) \delta_g(\mathbf{n}, z). \quad (3.12) $$

Eq. (3.12) relates the anisotropies in the energy density of an astrophysical SGWB to the anisotropies of the large-scale structure where SGWB sources are found. The bias $b_s(z, \mathcal{M})$, defined through eq. (3.11), relates the fractional galaxy anisotropies to the fractional source anisotropies. It does not depend on direction, since we assume that GW sources follow the same clustering structure as the galaxies, but it can depend on mass and redshift. Note that in eq. (3.10), we assumed that given a mass and a redshift, the clustering contribution shares the same frequency distribution as the background, so that $\delta_s$ does not depend on $f$. In other words, a source with parameters $(z, \mathcal{M})$ traces the galaxy identically no matter at what frequency it is emitting.

The spectral density of a SGWB given in eq. (2.9) is directly related to its energy content via the relation [40]

$$ \Omega_{\text{GW}} = c^2 f^3 S_h / (8G\rho_c). \quad (3.13) $$

In particular, for the anisotropic part, we arrive at

$$ \delta S_h = \delta\Omega_{\text{GW}} \frac{8G\rho_c}{f^3 c^2} = \bar{S}_h(f) \int dz b_{\text{GW}}(z) \delta_g(\mathbf{n}, z), \quad (3.14) $$

where we have used eq. (2.9) for the definition of $\delta S_h$. The form of the gravitational energy released by astrophysical sources was so far generic. We now focus on the case of an inspiraling BBH, for which we have [36]

$$ \frac{dE_s}{dt_s} = \frac{96}{15} \pi^{10/3} \frac{1}{Gc^5} (GM)^{10/3} f_s^{10/3} = \frac{96}{15} \pi^{10/3} \frac{1}{Gc^5} (GM)^{10/3} (1+z)^{10/3} f^{10/3}. \quad (3.15) $$

Combining eqs. (3.14), (3.15) and (3.12) we obtain

$$ b_{\text{GW}}(z) = \frac{1}{\bar{S}_h(f)} \int d\log(\mathcal{M}) b_s(z, \mathcal{M}) \alpha(z, \mathcal{M}, f), \quad (3.16) $$

where

$$ \alpha(z, \mathcal{M}, f) = \frac{64}{5} \frac{\pi^{7/3}}{c^7 H(z)} (GM)^{10/3} f^{4/3} (1+z)^{4/3} \bar{n}_s(z, \mathcal{M}, f). \quad (3.17) $$

Eq. (3.16) provides an expression for the gravitational kernel $b_{\text{GW}}(z)$ in terms of the distribution of sources $\bar{n}_s(z, \mathcal{M}, f)$, the mean spectral density $\bar{S}_h(f)$ and the source bias $b_s(z, \mathcal{M})$. We have so far allowed the $f$-dependency in the distribution of sources $\bar{n}_s$ to be generic. However, note that in the scenario of a distribution of supermassive binaries evolving via quasi-circular inspirals, the frequency distribution is, in the limit of many very massive sources [41],

$$ \bar{n}_s(z, \mathcal{M}, f) = \bar{n}_s(z, \mathcal{M}, f = 1\text{Hz}) \left(\frac{f}{1\text{Hz}}\right)^{-11/3}, \quad (3.18) $$

corresponding to the idea that the probability for a binary to be found in a frequency interval $\Delta f$ is proportional to the time it spends in the $\Delta f$ interval. In this case, we see that $b_{\text{GW}}(z)$
does not depend on frequency $f$, since in eq. (3.16) both $\alpha(z, \mathcal{M}, f)$ and $\bar{S}_h(f)$ scale as $f^{-7/3}$. As a consequence, $\bar{S}_h(f)$ and $\delta S_h(f)$ have the same frequency dependence, see eq. (2.9).<fnref n="12" />

We can extract the sky-averaged quantities by repeating the same calculation and substituting the quantity $\bar{n}_s(z, \mathcal{M}, f)b_s(z, \mathcal{M})\delta_g(\mathbf{n}, z)$ with the sky-averaged quantity $\bar{n}_s(z, \mathcal{M}, f)$ in eq. (3.12). This leads to

$$ \bar{S}_h(f) = \int dz \int d\log(\mathcal{M}) \alpha(z, \mathcal{M}, f). \qquad (3.19) $$

## 3.2 Reconstruction of $b_{\text{GW}}$ from catalogs

From eqs. (3.16) and (3.17), the contribution of source anisotropies to the variance of the HD correlation must rely on a model for the (isotropic) distribution in chirp mass and volume of inspiraling SMBHs, $\bar{n}_s(z, \mathcal{M}, f)$. This distribution follows from the distribution of individual SMBHs, combined with some modeling of the binary formation process. The former has observational constraints from electromagnetic signals of high redshift SMBH-powered quasars, other active galactic nuclei emission processes, as well as baryon dynamics, see refs. [42–46] for a selection of historical approaches. The theoretical modeling of this SMBH population is however a long-term astrophysical challenge closely tied to the mechanism via which such SMBH may arise, see refs. [47, 48] for reviews. In particular, from naive models it is hard to predict the observations of heavy SMBH at redshifts as early as $z = 6$ [49].

The merging process is a further challenge due to the wide range of scales involved, ranging from Mpc scale for galaxy halo mergers to $10^{-6}$pc for the subsequent SMBH merger. Owing to this difficulty, we rely here on the numerical catalogs of SMBH binaries described in refs. [29, 30]. These catalogs explore two different scenarios for the production of seeds of individual SMBH:

1. heavy ($\sim 10^5 M_\odot$) seeds following the collapse of protogalactic disks at $10 \lesssim z \lesssim 15$,

2. light ($\sim 150 M_\odot$) seeds following the collapse of population III stars at $z \sim 20$.

These seeds are then evolved within a model for galaxies that includes prescriptions for the evolution and mergers of dark matter halos, baryonic gas dynamics, and environmental effects, thereby providing a description of the co-evolution of SMBHs and their galactic hosts. We refer to refs. [29, 30] for more details on the catalog implementation.

For the case of seeds from the collapse of protogalactic disks, the catalog exists in two variants. In one case, a non-zero time delay $\Delta z_{\text{del}}$ between galaxy and SMBH merger is implemented, while in the other it is explicitly set to $\Delta z_{\text{del}} = 0$. These variants are named Q3-d and Q3-nd, with Q3 referring to the value of the Toomre stability parameter $Q = 3$ of the collapsing protogalactic disk, while d and nd respectively refer to the assumption of a delay or no delay. For the case of population III star seeds, the assumptions on the delay $\Delta z_{\text{del}}$ have been found to not appreciably alter the population of sources, hence the population III seeds model is only provided in one variant, labeled as popIII.

<footnote n="12">¹²Even in scenarios where $\bar{S}_h(f)$ does not scale as $f^{-7/3}$, i.e. where the GW background is not perfectly described by quasi-circular inspiraling SMBH mergers, we generally expect $\bar{S}_h(f)$ and $\delta S_h(f)$ to scale equally with frequency and thus $b_{\text{GW}}(z)$ to be frequency-independent. A frequency dependence of $b_{\text{GW}}(z)$ could only arise if regions with higher source density would exhibit a different frequency scaling as regions with lower source density.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The function $b_{\text{GW}}(z)$ is shown for the three different population models in consideration.</figcaption></figure>

In appendix B, we outline how to use these numerical catalogs of SMBH binaries to reconstruct the isotropic GW luminosity over $z$, assuming that binaries follow quasi-circular orbits. We then extract from it a numerical expression of the astrophysical kernel $b_{\text{GW}}$ of eq. (3.16). We plot our results in figure 1 for the three catalog models, which all show a similar behavior for the astrophysical kernel. From figure 1, we see that the astrophysical signal to which PTA are sensitive arises from sources at $z \lesssim 2$, with a peak at $z \sim 0.1$. The popIII model produces a catalog of lighter but more numerous binary sources compared to the Q3 models, such that, after performing the integral in eq. (3.16), their astrophysical kernels become comparable. This is true whether or not one assumes $\Delta z_{\text{del}} = 0$.

## 3.3 Discussion on the source bias and galaxy bias

In our numerical computation of $b_{\text{GW}}$ from eq. (3.16) we have set the source bias to $b_s(z, \mathcal{M}) = 1$. This source bias encodes a) the way the distribution of GW sources (inspiraling SMBH in our case) traces the distribution of galaxy mergers, and b) the way the distribution of galaxy mergers traces the distribution of galaxies. The merger rate of two galaxies, resulting from the merger of the two host dark matter halos, is actually mass dependent as can be predicted e.g. with the extended Press-Schechter formalism [50]. More specifically, it is higher for higher mass galaxies and halos, and in principle depends on redshift as well. From the Millennium numerical simulations, this dependency has been found to be of the type $\sim (1+z)^{0.1} M_h^{0.13}$, with $M_h$ being the mass of the merged halo [51]. Given this weak dependency, one may simplify the description by assuming that galaxy mergers trace the galaxy population identically at all masses and redshifts.

Using this simplification, the galaxy fractional number overdensity $\delta_g$ then equals the fractional overdensity of galaxy merger number counts, as no significant bias arises from b), i.e. the way the distribution of galaxy mergers traces the distribution of galaxies. We further note that the overdensity of galaxy merger counts corresponds to the fractional number overdensity of inspiraling SMBH $\delta_s^{\text{clust.}}$, with a possible delay $\Delta z_{\text{del}}$ owing to the typical time delay between galaxy merger and SMBH merger, such that $\delta_s^{\text{clust.}}(\mathbf{n}, z) \simeq \delta_g(\mathbf{n}, z + \Delta z_{\text{del}})$. Setting $\Delta z_{\text{del}} = 0$, assuming that the time span of such mergers is small compared to the cosmological time span over which the galaxy overdensity evolves, we arrive at $\delta_s^{\text{clust.}}(\mathbf{n}, z) \simeq \delta_g(\mathbf{n}, z)$. In other words, also no significant bias arises from a), i.e. the way the distribution of galaxy
mergers traces the distribution of galaxies, and we arrive at $b_s(z, \mathcal{M}) \simeq 1$. The fact that, as seen above, catalogs with and without delay $\Delta z_{\text{del}}$ lead to similar values of the mass integral over $\alpha$ indicates that this approximation is reasonable.

Another bias factor that also affects our numerical results is the galaxy bias $b_g(z)$. Indeed, the quantity that we use in eq. (2.10) to describe the impact of galaxy clustering on the HD correlation is the galaxy correlation function $\xi_g(\mathbf{n} \cdot \mathbf{n}', z, z')$. However, cosmological models instead make prediction on the matter correlation function, $\xi_m(\mathbf{n} \cdot \mathbf{n}', z, z')$. Since galaxy over- or underdensities are not perfect tracers of the underlying matter field, these quantities are related via the galaxy bias:

$$ \xi_g(\mathbf{n} \cdot \mathbf{n}', z, z') = b_g^2(z) \xi_m(\mathbf{n} \cdot \mathbf{n}', z, z'). \qquad (3.20) $$

The galaxy bias itself can be expressed as a number-weighted integral over the halo bias $b_h(M_h)$, depending on halo mass $M_h$ [52]. Thus, our astrophysical modeling would in principle need to be extended by choosing a model for the halo bias (see e.g. figure 15 in ref. [52]) and relating halo masses to central BH masses as outlined in ref. [53]. This bias, depending on BH mass, could then be included in eq. (3.16) for $b_{\text{GW}}$, and the integral over chirp mass $\mathcal{M}$ would be replaced by integrals over the individual BH masses available from the catalogs [29, 30]. However, as we can infer from figure 15 in ref. [52], halo bias reaches a maximum value of 4 – 5 for the most massive halos only, and thus we expect the overall boost of our signal due to this bias to be below this value. In the following, we set the galaxy bias to 1 for simplicity. As we will see, since the clustering variance is roughly two orders of magnitude smaller than the standard variance, even a boost of a factor 5 would not make this contribution relevant for observations.

# 4 Numerical results for the clustering variance

Having established our astrophysical modeling, we are now equipped to perform numerical evaluations for the galaxy clustering variance of the HD correlation. In the following, we will first describe our numerical methods and then present our results. A comparison to ref. [33], a recent work with similar considerations, is provided as well.

## 4.1 Methods

The galaxy clustering variance in eq. (2.10) poses considerable challenges for numerical evaluation, since it contains a 6 dimensional integral (four over angles and two over redshifts) of the galaxy correlation function $\xi_g(\mathbf{n} \cdot \mathbf{n}', z, z')$, weighted by the highly non-trivial functions $\chi_{ab}(\mathbf{n})$. In ref. [12], we presented alternative expressions using either the galaxy power spectrum $P_g(k, z, z')$ or the angular power spectrum $C_l(z, z')$ defined via

$$ \begin{gathered} \langle \delta_g(\mathbf{k}, z) \delta_g^*(\mathbf{k}', z') \rangle = (2\pi)^3 \delta_D(\mathbf{k} - \mathbf{k}') P_g(k, z, z'), \\ \xi_g(\mathbf{n} \cdot \mathbf{n}', z, z') = \frac{1}{4\pi} \sum_l (2l + 1) C_l(z, z') P_l(\mathbf{n} \cdot \mathbf{n}'), \end{gathered} \qquad (4.1) $$

where $\delta_g(\mathbf{k}, z)$ is the Fourier transform of the galaxy overdensity and $P_l$ represents the $l$-th Legendre polynomial. Indeed, the integral over $\mathbf{n}$ and $\mathbf{n}'$ in the second line of eq. (2.10) has a
convenient analytical solution when using the angular power spectrum $C_l(z, z')$. We find that

$$\iint \frac{dn}{4\pi} \frac{dn'}{4\pi} P_l(\mathbf{n} \cdot \mathbf{n}') \chi_{aa}(\mathbf{n}) \chi_{bb}(\mathbf{n}') = \frac{1}{9} \left( \delta_{l0} + \frac{1}{4} \cos \gamma \delta_{l1} + \frac{1}{400} (1 + 3 \cos(2\gamma)) \delta_{l2} \right), \quad (4.2)$$

where $\gamma$ denotes the angle between the two pulsars at a and b.<fnref n="13" /> As a consequence, the second line in eq. (2.10) only contributes to the variance if a monopole ($l=0$), dipole ($l=1$) or quadrupole ($l=2$) in the galaxy overdensities would be present. However, the monopole is set to zero by construction (any non-zero monopole can be reabsorbed in the strain and thus in the normalization of the HD curve), while contributions at ultra-large angular scales $l=1, 2$ are small unless some mechanism breaking the cosmological principle is present in the universe.<fnref n="14" /> Here, we stay within the scope of the standard $\Lambda$CDM cosmological model obeying the cosmological principle. This means that the second line of eq. (2.10) cannot lead to any sizable contribution, and we thus focus on the first line.

In the following, we numerically evaluate the remaining term in eq. (2.10), having contributions from all galaxy clustering multipoles. To compute this term numerically, it is more convenient to rewrite the correlation function in terms of the power spectrum $P_g(k, z)$. We can then apply the Limber approximation [62, 63] to reduce from two to one integral over redshift,

$$\langle \rho_{ab}^2 \rangle^{\text{clust.}} = \int dz [G^2(z) + \Gamma(z, z)] \frac{H(z)}{c} \iint \frac{dn}{4\pi} \frac{dn'}{4\pi} \chi_{ab}(\mathbf{n}) \chi_{ab}(\mathbf{n}') \\ \times \int \frac{k \, dk}{2\pi} P_g(k, z) J_0(kr_z \cos^{-1}(\mathbf{n} \cdot \mathbf{n}')) . \quad (4.3)$$

Here, $J_0$ is the 0th order Bessel function, $H(z)$ is the Hubble parameter at redshift $z$, and $r_z$ denotes the comoving distance given by $r_z = \int_0^z dz' c H(z')^{-1}$. The Limber approximation has been shown to work well for galaxy clustering at small separation and in the case of wide redshift bins [64]. Moreover, it is more precise at low than at high redshifts [64, 65]. Since the gravitational kernel $b_{\text{GW}}$ has support up to redshift $z \sim 2$ and peaks at a low redshift $z \sim 0.1$ (see figure 1), we expect the Limber approximation to have a high precision. Only for very large angular separations, this approximation will not be accurate. However, at those large scales the amplitude of the galaxy correlations is very small, and we expect no significant contribution to the variance from such less accurate modes.

Our numerical evaluation is performed using the following steps:

1. First, we determine the galaxy power spectrum $P_g(k, z)$ at redshifts $\{0, 0.001, 0.01, 0.05, 0.1, 0.5, 1, 5, 10\}$ by using the publicly available<fnref n="15" /> CLASS code [66] to obtain the matter

<footnote n="13">¹³This expression can be shown by first choosing sky coordinates such that the pulsar a aligns with the z-direction, $\mathbf{n}_a = \mathbf{e}_z$, and the pulsar b lies on the $x-z$ axis, $\mathbf{n}_b = \sin \vartheta \mathbf{e}_x + \cos \vartheta \mathbf{e}_z$. Then, the quantities $\rho_{aa}(\mathbf{n})$ and $\rho_{bb}(\mathbf{n}')$ as well the Legendre polynomial $P_l(\mathbf{n} \cdot \mathbf{n}')$ can be decomposed into spherical harmonics (similar in spirit to ref. [34]), and applying the orthogonality of spherical harmonics finally leads to eq. (4.2).</footnote>

<footnote n="14">¹⁴Some observations indicate the presence of structure at ultra large angular scales, such as the cosmic dipole that has been measured from the overabundance of radio surveys [54–58]. However, this anomaly might be due to systematic effects and independent measurements of the cosmic dipole, e.g. with ground based GW observatories [59–61], could bring more clarity.</footnote>

<footnote n="15">¹⁵https://github.com/lesgourg/class_public.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left panel: Galaxy power spectrum $P_g(k, z)$ (set equal to the matter power spectrum, ignoring galaxy bias) as a function of $k$ for 3 redshift values $z = 0.01, 0.5, 1$. Right panel: The function $F(d, z)$ at the same 3 redshifts. Dashed lines denote negative values.</figcaption></figure>

power spectrum,<fnref n="16" /> and ignoring the bias $b_g$ between matter and galaxy correlations as discussed above. We assume cosmological parameters measured from the Planck survey [69]. While we consider a limited amount of redshifts, we will see below that this choice is sufficient to describe the redshift-behavior of the integrals over angles $\mathbf{n}$ and $\mathbf{n}'$, which appear roughly linear on a double-logarithmic scale. The galaxy power spectrum $P_g(k, z)$, as a function of $k$ and for three values of $z$, is plotted on the left panel of figure 2.<fnref n="17" />

2. The integral over $k$ in eq. (4.3) to obtain the function

$$F(d, z) = \int \frac{k \, dk}{2\pi} P_g(k, z) J_0(kd), \qquad (4.4)$$

is then performed using an FFTLog code.<fnref n="18" /> We show the result, as a function of $d$ for three values of $z$, in the right panel of figure 2.

3. Then, we numerically evaluate the integral over angles $\mathbf{n}$ and $\mathbf{n}'$,

$$\mathrm{Int}(\gamma, z) = \iint \frac{\mathrm{d}\mathbf{n}}{4\pi} \frac{\mathrm{d}\mathbf{n}'}{4\pi} \chi_{ab}(\mathbf{n}) \chi_{ab}(\mathbf{n}') F(r_z \cos^{-1}(\mathbf{n} \cdot \mathbf{n}'), z), \qquad (4.5)$$

for values {$0.01, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1$}π of separation $\gamma$ between the pulsars a and b. We perform a Monte Carlo integration 1000 times using 10⁷ sampling points each, giving us both a mean and an error estimate (well within the marker size in

<footnote n="16">¹⁶The CLASS code provides the options to use either the linear or non-linear matter power spectrum. Here, we choose the non-linear option, as it is known to work well for high-mass halos [67], corresponding to high-mass central BHs which contribute more to the GW signal. Another possible approach would be to use the linear power spectrum and then apply a non-linear halo bias [68], although we choose to avoid such complications given the overall small amplitude of the clustering variance.</footnote>

<footnote n="17">¹⁷The power spectrum is typically given in units of $(\mathrm{Mpc}/h)^3$, where $h$ (unrelated to the GW strain) sets the value of the Hubble parameter today, $H_0 = 100h \, \mathrm{km/s/Mpc}$. Our final results do not depend on $h$, as this factor cancels out in eq. (4.3).</footnote>

<footnote n="18">¹⁸https://github.com/JCGoran/fttlog-python.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Results for Int($\gamma$, $z$) defined in eq. (4.5) obtained from Monte Carlo integration at various pulsar pair separations $\gamma$, plotted as a function of redshift. For better comprehensibility of this figure, we omit the pulsar pair separations $\{0.1, 0.3, 0.5, 0.7, 0.9\}\pi$ and the redshift $z = 0.0001$, where no unexpected behavior has been found.</figcaption></figure>

figure 3).<fnref n="19" /> The results Int($\gamma$, $z$) are plotted in figure 3 as a function of $z$ for the various values of $\gamma$. We see that the functions exhibit a smooth redshift behavior. Therefore, no increase in the amount of redshifts is needed to obtain a reliable estimate of the galaxy clustering variance. We also note that these functions drop off fast with increasing $z$. Thus, we will find that ~95% of the clustering contribution to the variance arises from redshifts $z \le 0.1$, where the astrophysical kernel $b_{\text{GW}}(z)$ roughly peaks.

4. Finally, to perform the integral over redshift in eq. (4.3), we need to evaluate the kernels $G(z)^2$ and $\Gamma(z, z)$. These depend on the function $b_{\text{GW}}(z)$ via eqs. (2.12). After obtaining $b_{\text{GW}}(z)$ from the catalogs as described in section 3.2, we can combine with the redshift-interpolated functions Int($\gamma$, $z$) (eq. (4.5)), see the solid curves in figure 3) and perform the integral over redshift using standard quadrature integration methods.

Finally, we note that the expectation value for the HD correlation itself is given by the HD curve times the squared strain $h^2$ as a normalization factor, see eq. (2.6). Therefore, in the numerical results presented, all variance contributions are divided through this normalization factor for comparison with the HD curve itself.

## 4.2 Results

In the left panel of figure 4, we show the standard deviation caused by galaxy clustering in percent of the HD curve, for the three different astrophysical models in consideration. In the right panel, we compare the different variance contributions to the HD curve itself: for all pulsar pair separations, the standard variance is the dominant contribution, of the order of the HD curve itself, while the clustering variance reaches at most percent-level. These results are shown for the Q3-d model, but the same holds for all the three astrophysical models under

<footnote n="19">¹⁹For each value of $\gamma$ and $z$, this computation needs roughly 3 hours runtime with the usage of 10 CPUs on the Yggdrasil cluster of the University of Geneva.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left panel: We show the standard deviation caused by galaxy clustering, divided by the squared strain $h^2$, in percent of the HD curve for the three different astrophysical models in consideration. Right panel: We show the absolute value of the HD curve, the size of the standard deviation (square root of the variance) due to the standard variance calculated in ref. [7], along with the galaxy clustering contribution (for the Q3-d model) calculated in this work.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as figure 4, but for the case of infinitely many pulsar pairs instead of on pair; Left panel: We show the *pulsar-averaged* standard deviation caused by galaxy clustering, divided by the squared strain $h^2$, in percent of the HD curve for the three different astrophysical models in consideration. Right panel: We show the absolute value of the HD curve, the size of the standard deviation (square root of the variance) due to the *pulsar-averaged* standard variance calculated in ref. [7], along with the *pulsar-averaged* galaxy clustering contribution (for the Q3-d model) calculated in this work.</figcaption></figure>

study. Indeed, we see from the left panel that the differences between the three models are small, at most of the order of 10%, with the Q3-d model leading to a slightly larger clustering variance than the other two models. We have also evaluated that the term proportional to $G(z)^2$ in eq. (4.3) is about 4 times larger than the term proportional to $\Gamma(z, z)$. Since the term proportional to $G(z)^2$ is independent of observational time $T$, stochasticity remains even in the limit of infinite observational time due to galaxy clustering, i.e. due to the fact that the GW intensity is not uniform on the sky.

In appendix A, we show in detail the calculation of the pulsar-averaged galaxy clustering variance (obtained by averaging over an infinite number of pulsar pairs, at fixed angular separation). Here, we present the numerical results for this case in figure 5. As can be seen in
<table><caption>**Table 2.** Standard deviation in percent of the HD curve caused by the various contributions summarized in table 1. The 11 different pulsar pair separations are given in units of $\pi$. For the clustering variance and pulsar-averaged clustering variance, we assume the Q3-d model, leading to slightly higher results than the popIII and Q3-nd models (see left panels of figures 4 and 5). We only report 11 distinct pulsar pair separations, as the integrals over angles for the clustering variance and pulsar-averaged clustering variance (eqs. (4.3) and (A.11), respectively) are performed numerically.</caption><thead><tr><th>Pair separation</th><th>0.01</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.8</th><th>0.9</th><th>1.0</th></tr></thead><tbody><tr><td>Standard variance</td><td>56</td><td>74</td><td>181</td><td>629</td><td>184</td><td>174</td><td>314</td><td>938</td><td>183</td><td>116</td><td>103</td></tr><tr><td>Clustering variance</td><td>0.63</td><td>0.62</td><td>0.91</td><td>2.84</td><td>0.97</td><td>0.88</td><td>1.23</td><td>2.85</td><td>0.69</td><td>0.56</td><td>0.54</td></tr><tr><td>Irr. standard variance</td><td>7</td><td>8</td><td>12</td><td>9</td><td>9</td><td>12</td><td>16</td><td>13</td><td>12</td><td>14</td><td>15</td></tr><tr><td>Irr. clustering variance</td><td>0.10</td><td>0.11</td><td>0.13</td><td>0.11</td><td>0.12</td><td>0.13</td><td>0.14</td><td>0.12</td><td>0.13</td><td>0.13</td><td>0.13</td></tr></tbody></table>

table 2, we find a suppression by factors 4 – 26 for the clustering variance, depending on the pulsar separation (and assuming the Q3-d model, although results are similar for all models). Hence, performing the pulsar pair averaging largely counterbalances the stochasticity arising from the anisotropy in GW intensity. More specifically, the HD curve would be exactly recovered in the limit of infinitely many pulsar pairs and an infinitely large observational time. The remaining term $\propto \Gamma(z, z)$, which depends on observational time, is as well altered comparing the single pulsar pair and pulsar-averaged case. In particular, the dependence on the pulsar pair separation is not the same, leading to the variation in suppression factors. In practice, the variance of the HD correlation in a realistic PTA experiment (with a finite number of pulsar pairs for each angular separation) will lie in between the results of figure 4 and the pulsar-averaged variance results in figure 5.

We stress that while the equations that we derived are generic, some astrophysical simplifications have been made in our analysis. More specifically, we have assumed no galaxy bias, and that galaxy mergers directly trace mergers of central black holes. Hence, SMBH mergers are treated as perfect tracers of the underlying matter field in our computation. However, as discussed in section 3.3, we do not expect an increase of the clustering variance to a level comparable to the standard variance (which is roughly 2 orders of magnitude larger than the clustering variance, see table 2) by lifting these assumptions. Likewise, while we did not include a frequency sensitivity function, it is unlikely that such detector specifications would substantially change the order of magnitude of the clustering variance compared to the standard variance. Moreover, even though we applied a standard procedure for large-scale structure observables, small errors could e.g. arise due to the fact that the Limber approximation does not precisely describe clustering on very large angular scales. In the case of galaxy clustering, the Limber approximation has been well-studied and is known to work particularly well for broad bins at low redshift, reaching percent-level precision at angular scales $l \gtrsim 10$ (see e.g. figure 2 of ref. [64]). Even at the largest scales, $l < 10$, the difference between the Limber approximated and precise results is at most at the level of 10%. Thus, as the Limber approximation has a high precision for large-scale structure observables, we do not expect it to have a large impact on our results. Another imprecision arises as we apply a sparse sampling in $z$ of the functions $\text{Int}(\gamma, z)$. Since they show a smooth behavior on a
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left panel: Angular power spectrum $C_\ell(z = z')$ of matter density correlations at various values of $z$. The constant value $C_\ell = 0.002$ and the cut-off of $\ell_{\max} = 50$ suggested in ref. [33] is shown as well. Right panel: Same as left panel, but with the angular power spectra multiplied by the factor $(2\ell + 1)$.</figcaption></figure>

double-logarithmic scale, the resulting numerical errors are expected to be small. To conclude, given that the galaxy clustering variance remains significantly below the standard variance, all simplifications made in our analysis are well justified.

## 4.3 Comparison to Allen et al. (2024)

A recent work, ref. [33], has as well studied the impact of galaxy clustering on the HD correlation, applying a spherical harmonic analysis. In particular, they choose to work with the angular power spectrum $C_\ell$ of the matter correlations (see eq. (4.1)). An expression for an upper limit of the clustering contribution in the pulsar-averaged case has been derived under the crude assumption that $C_\ell$ is constant for all $\ell$. Then, setting this constant value to $C_\ell = 0.002$, an upper bound of one percent has been found. Employing a simple model where this constant value is cut off and set to zero above $\ell_{\max} = 50$, the authors argue that the galaxy clustering contribution needs to stay below this percent level. In the left panel of figure 6, we plot the $C_\ell$ obtained from CLASS at redshifts $z = 0.1, 0.2, 0.4$ (setting $z' = z$ in eq. (4.1)). The horizontal line corresponding to a constant $C_\ell = 0.002$ is as well included, along with the dashed vertical line corresponding to a cut-off at $\ell = 50$. We see that, for the lower two redshifts, the $C_\ell$ are above the horizontal line for low angular multipoles. Let us recall that, due to the steep decline of the integrals $\text{Int}(\gamma, z)$ with $z$ and the fact that the astrophysical kernel $b_{\text{GW}}(z)$ peaks at $z \sim 0.1$, about 95% of the contribution to the clustering variance arises from redshifts $z \le 0.1$. The dependence on redshift is not evident in ref. [33], while in principle, the clustering contribution should be obtained as the integral over $z$ and $z'$ (see eq. (36) in ref. [12]), taking the redshift dependence of both $C_\ell(z, z')$ and $b_{\text{GW}}(z)b_{\text{GW}}(z')$ into account. Therefore, a constant $C_\ell = 0.002$ does not properly capture the redshift-dependence and the angular dependence of the angular power spectrum.

We also note that the angular power spectrum $C_\ell$ cannot be cut at any low multipole. Indeed, while the $C_\ell$ appear to be strictly decreasing with higher $l$, the strength of galaxy correlation is affected by the product $(2\ell + 1)C_\ell$, as can be seen from eq. (4.1). This product,
shown in the right panel of figure 6, increases up to $l \approx 10^3$, reflecting the fact that correlations
are stronger at small angular separation (i.e. large $l$). For this reason, a priori we do not
expect low multipoles $l \lesssim 50$ to be the major contribution to the clustering contribution.
All multipoles should be consistently included, weighted by the angular dependence of the
functions $\chi_{ab}(\mathbf{n})$ that will determine the impact of different multipoles (see eq. (36) in ref. [12]
or eq. (6.2) in ref. [33]). In our case, we work with the galaxy power spectrum in $k$-space
instead of the angular power spectrum, and we take the full shape and redshift-dependence
of the matter power spectrum $P_m(z, k)$ into account. With this we can precisely compute the
galaxy clustering variance, which we find to be at the level of 0.1% in the pulsar-averaged
case, i.e. one order of magnitude below the upper bound found in [33]. Our methods are
readily extendable to study more practical statistical tools to extract the galaxy clustering
information from PTA surveys. For this future task, the scale-dependence of galaxy clustering
will be crucial.

# 5 Conclusion

In this work, based on the theoretical formalism established in ref. [12], we performed
numerical calculations for the galaxy clustering contribution to the variance of the HD
correlation. For the case of one pulsar pair per angular separation, we have found this
contribution to be at the level of 0.5 – 3%, depending on the pulsar pair separation. This is
substantially smaller than the standard variance, calculated in the absence of anisotropies,
which is of the order of the HD curve itself.

In our work we considered three different astrophysical models to describe the distribution
of SMBH GW sources, and we found qualitatively similar results for all three models. Our
computation in the main part of this work is done for the case of one pulsar pair for each
separation, hence the estimate we provide for the galaxy clustering contribution to the
variance are an upper bound to what one would obtain in realistic PTA surveys, which
typically contain more than one pulsar pair at a given separation. We have also calculated,
in appendix A, the pulsar-averaged galaxy clustering variance, obtained by averaging over an
infinite number of pulsar pairs at each angular separation. In this case, providing a lower
limit for the variance in a realistic observation, both variance contributions (standard and
clustering) are suppressed by roughly an order of magnitude compared to the case of one
pulsar pair. As a consequence, the galaxy clustering contribution is small compared to the
standard contribution as well in this case. Realistic estimates of the variance for finite pulsar
pair numbers will lie in between the single pair and infinite pair case (see figure 9 in ref. [8]),
but the exact number of pulsar pairs observed will not change the fact that the clustering
contribution to the variance remains subdominant with respect to the standard variance.

Potential sources of astrophysical and numerical errors are discussed at the end of
section 4, although none of them are expected to lift this contribution above the percent
level. Therefore, we can safely conclude that the clustering variance is not a major source
of uncertainty for current and future PTA surveys. In fact, even if future surveys make a
significant improvement in precision with an increase of pulsar pairs and observational time
(see e.g. ref. [70]), the variance on the HD correlation will be dominated by the standard
variance (calculated in ref. [7] for an isotropic distribution of GW sources).
In this work we have treated GW source clustering as a source of contamination. However, it is interesting to rather consider it as a signal: the fact that sources are clustered leads to anisotropies in the observed energy density. Different source populations are expected to lead to different levels of anisotropy, hence anisotropies are an interesting observable to look at to distinguish different scenarios for the origin of the SGWB. On the observational side, mapping techniques are already in place, and forecasts have shown that the ability of PTAs to map a background will significantly increase in the future, as the number of pulsar monitored and the observation time will increase [71]. Predictions for the angular power spectrum for different astrophysical scenarios are under study. Moreover, another interesting avenue lies in the cross-correlation between PTAs and large-scale structure tracers, see ref. [72] for a simulation-based approach. Thus, while galaxy clustering does not appreciably affect the variance of the HD correlation, further investigations on synergies between PTAs and galaxy clustering may lead to promising results.

## Acknowledgments

We thank Ruth Durrer, Gabriele Franciolini, Cyril Pitrou, Alberto Sesana, Gianmassimo Tasinato and Marta Volonteri for valuable discussions. We are grateful to Maria Berti, Enea Di Dio and Francesco Iacovelli for useful advice on numerical methods. We acknowledge the use of the HPC cluster Yggdrasil at the University of Geneva for conducting our numerical evaluations. N.G. and C.B. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 863929; project title "Testing the law of gravity with novel large-scale structure observables"). The work of G.C. is supported by CNRS and G.C. and M.P. acknowledge support from the Swiss National Science Foundation (Ambizione grant, "Gravitational wave propagation in the clustered universe").

## A Pulsar-averaged galaxy clustering variance

In section 2.4, we have presented results on the variance of the HD correlation when considering a single pulsar pair at fixed positions $n_a$ and $n_b$. However, one can as well first average over all possible locations $n_a$ and $n_b$ separated by a fixed angle $\gamma$, see ref. [7], corresponding to the limit of a survey with infinitely many observed pulsars. Thereby, one obtains a different observable, with expectation value as well given by the standard HD curve but having a reduced variance. More precisely, the HD correlation $\rho_{ab}$ depending on the pulsars $a$ and $b$ is replaced by the pulsar averaged correlation $\rho(\gamma)$,<fnref n="20" />

$$ \rho(\gamma) = \langle \rho_{ab} \rangle_p = \langle \overline{Z_a Z_b} \rangle_p, \qquad (A.1) $$

<footnote n="20">²⁰A possible source of confusion could arise from the fact that while $\rho(\gamma)$ depends only on the pulsar pair separation $\gamma$, this is true for the quantity $\rho_{ab}(\gamma)$ as well. Indeed, if we consider only one pulsar pair $a$ and $b$, we can choose a coordinate system such that the pulsar $a$ is aligned with the $z$-direction, while $b$ lies in the $x-z$ plane. However, the quantity $\rho(\gamma)$ represents an average over all possible pulsar pairs. As soon as more than one pair is involved, a unique choice of coordinate system does not allow to align one pulsar in each pair with the $z$-direction. Therefore, even though both $\rho(\gamma)$ and $\rho_{ab}(\gamma)$ only depend on $\gamma$, these two are distinct quantities.</footnote>
where $\langle \dots \rangle_p$ denotes the integral over all sky directions $\mathbf{n}_a$ and $\mathbf{n}_b$ at a fixed angular separation. In ref. [7], the variance of $\rho(\gamma)$ is calculated in the absence of inhomogeneities.<fnref n="21" /> Here, we calculate the additional galaxy clustering contribution (see also refs. [34, 35] for similar recent work).

Following ref. [7], we omit the exponential term in the round brackets of eq. (2.2), since this term will be rapidly oscillating and average to zero for pulsar distances much larger than $1/f$. Then, $\rho(\gamma)$ reads:

$$ \begin{aligned} \rho(\gamma) &= \langle \overline{Z_a Z_b} \rangle_p \\ &= \sum_{A,A'} \int df \int df' \int d\mathbf{n} \int d\mathbf{n}' \langle F_a^A(\mathbf{n}) F_b^{A'}(\mathbf{n}') \rangle_p h_A^*(f, \mathbf{n}) h_{A'}(f', \mathbf{n}') \mathrm{sinc}(\pi(f-f')T). \end{aligned} \quad (A.2) $$

For brevity of notation, one can define

$$ \mu_{AA'}(\gamma, \beta(\mathbf{n}, \mathbf{n}')) \equiv \langle F_a^A(\mathbf{n}) F_b^{A'}(\mathbf{n}') \rangle_p, \quad (A.3) $$

where $\beta(\mathbf{n}, \mathbf{n}') = \cos^{-1}(\mathbf{n} \cdot \mathbf{n}')$ denotes the separation angle between the two GW sources at $\mathbf{n}$ and $\mathbf{n}'$. This quantity reflects the fact that the antenna pattern functions depend on the exact geometrical constellation, i.e. the positions of both pulsars and the two GW source directions. After varying the pulsar pair, at fixed angular separation $\gamma$, over the sky, this dependence is reduced to a dependence on the two separation angles $\gamma$ and $\beta$ only.

Starting from eq. (A.2), we can repeat the procedure in section III of ref. [12], replacing any factors $R_a^A(f, \mathbf{n}) R_b^{A'}(f', \mathbf{n}')$ by $\mu_{AA'}(\gamma, \beta(\mathbf{n}, \mathbf{n}'))$. For the standard terms, i.e. without the impact of inhomogeneities, this leads to

$$ \begin{aligned} \langle \rho^2 \rangle^{\text{st.}} &= \frac{1}{4} \iint df \, df' \bar{S}_h(f) \bar{S}_h(f') \iint \frac{d\mathbf{n}}{4\pi} \frac{d\mathbf{n}'}{4\pi} \sum_{A,A'} \mu_{AA}(\gamma, 0) \mu_{A'A'}(\gamma, 0) \\ &\quad + \frac{1}{2} \iint df \, df' \bar{S}_h(f) \bar{S}_h(f') \mathrm{sinc}^2(\pi(f-f')T) \\ &\quad \times \iint \frac{d\mathbf{n}}{4\pi} \frac{d\mathbf{n}'}{4\pi} \sum_{A,A'} \mu_{AA'}(\gamma, \beta(\mathbf{n}, \mathbf{n}')) \mu_{AA'}(\gamma, \beta(\mathbf{n}, \mathbf{n}')), \end{aligned} \quad (A.4) $$

corresponding to eq. (C44) in ref. [7]. For the integral over angles in the first line, we note that

$$ \sum_A \mu_{AA}(\gamma, 0) = \sum_A \langle F_a^A(\mathbf{n}) F_b^A(\mathbf{n}) \rangle_p = \mu_{\mathrm{HD}}(\gamma) \quad (A.5) $$

is the usual HD curve (corresponding to having one GW point source and averaging over pulsar pair positions, see section II in ref. [7]). Therefore, the first line in eq. (A.4) vanishes when subtracting the square of the mean, and the second line constitutes the *pulsar-averaged standard variance* (an analytical solution for the integral over angles is given in eq. (G11) of ref. [7]).

<footnote n="21">²¹In ref. [7] the standard variance of $\rho(\gamma)$ (denoted by $\Gamma(\gamma)$ therein), i.e. without taking the presence of inhomogeneities into account, is referred to as “cosmic variance”, and the difference between the variance of $\rho_{ab}(\gamma)$ and $\rho(\gamma)$ is referred to as “pulsar variance”. Here, to avoid ambiguity with large-scale structure terminology, we refer to the variance of $\rho(\gamma)$ as the pulsar-averaged variance, occurring after averaging over infinitely many pulsar pairs.</footnote>
Let us now apply the same considerations to the clustering terms. We obtain:

$$ \langle \rho^2 \rangle^{\text{clust.}}_{h|\rho_g} = \iint dz \, dz' G(z) G(z') \iint \frac{dn}{4\pi} \frac{dn'}{4\pi} \delta_g(\mathbf{n}, z) \delta_g(\mathbf{n}', z') \sum_{A,A'} \mu_{AA}(\gamma, 0) \mu_{A'A'}(\gamma, 0) \\ + 2 \iint dz \, dz' \Gamma(z, z') \iint \frac{dn}{4\pi} \frac{dn'}{4\pi} \delta_g(\mathbf{n}, z) \delta_g(\mathbf{n}', z') \\ \times \sum_{A,A'} \mu_{AA'}(\gamma, \beta(\mathbf{n}, \mathbf{n}')) \mu_{A'A'}(\gamma, \beta(\mathbf{n}, \mathbf{n}')). \quad (A.6) $$

which we can further write as

$$ \langle \rho^2 \rangle^{\text{clust.}}_{h|\rho_g} = \mu_{\text{HD}}^2(\gamma) \iint dz \, dz' G(z) G(z') \iint \frac{dn}{4\pi} \frac{dn'}{4\pi} \delta_g(\mathbf{n}, z) \delta_g(\mathbf{n}', z') + \iint dz \, dz' \Gamma(z, z') \\ \times \iint \frac{dn}{4\pi} \frac{dn'}{4\pi} \delta_g(\mathbf{n}, z) \delta_g(\mathbf{n}', z') [\mu^2(\gamma, \beta(\mathbf{n}, \mathbf{n}')) + \mu^2(\pi - \gamma, \pi - \beta(\mathbf{n}, \mathbf{n}'))]. \quad (A.7) $$

Here $\mu(\gamma, \beta)$ is defined as

$$ \mu(\gamma, \beta) \equiv \mu_{++}(\gamma, \beta) + \mu_{\times\times}(\gamma, \beta), \quad (A.8) $$

and we have applied the parity transformation properties $\mu_{++}(\gamma, \beta) = \mu_{++}(\pi - \gamma, \pi - \beta)$ and $\mu_{\times\times}(\gamma, \beta) = -\mu_{\times\times}(\pi - \gamma, \pi - \beta)$. Taking the full ensemble average and performing a change of variables, this leads to

$$ \langle \rho^2 \rangle^{\text{clust.}} = \frac{1}{2} \mu_{\text{HD}}^2(\gamma) \iint dz \, dz' G(z) G(z') \int_0^\pi d\beta \sin\beta \xi_g(\cos\beta, z, z') \\ + \frac{1}{2} \iint dz \, dz' \Gamma(z, z') \int_0^\pi d\beta \sin\beta \xi_g(\cos\beta, z, z') [\mu^2(\gamma, \beta) + \mu^2(\pi - \gamma, \pi - \beta)]. \quad (A.9) $$

We note that the first term, proportional to the squared HD curve itself, is only vanishing if

$$ \int_0^\beta d\beta \sin\beta \xi_g(\cos\beta, z, z') \neq 0, \quad (A.10) $$

or equivalently (considering the first line in eq. (A.7)), if $\int dn \delta_g(\mathbf{n}, z) \neq 0$. This term would correspond to a monopole fluctuation in the galaxy overdensity, and we omit it in our calculation. Cosmological fluctuations at very large angular scales are predicted to be small in a $\Lambda$CDM universe. Moreover, a monopole fluctuation in the number of galaxies and thus the number of GW sources, increases the sky-averaged spectral density $\bar{S}_h$ and can be reabsorbed in the GW strain.<fnref n="22" /> In other word, it increases or decreases the overall amplitude of the HD correlation at any pulsar pair separation, but does not alter its shape, which is reflected in the proportionality to $\mu_{\text{HD}}(\gamma)^2$. We therefore focus on the second line, which, using the Limber approximation as before, can be simplified as

$$ \langle \rho^2 \rangle^{\text{clust.}} = \frac{1}{2} \iint dz \, \Gamma(z, z) \frac{H(z)}{c} \int_0^\pi d\beta \sin\beta [\mu^2(\gamma, \beta) + \mu^2(\pi - \gamma, \pi - \beta)] \\ \times \int \frac{k \, dk}{2\pi} P_g(k, z) J_0(kr_z\beta). \quad (A.11) $$

<footnote n="22">²²This is analogous to the usual procedure in galaxy surveys, where galaxy overdensity is defined with respect to the average number density observed on the sky. This procedure is absorbing any monopole fluctuation. Indeed, such monopole would be completely degenerate with the number density of galaxies that would form in an isotropic universe.</footnote>
This integral can be solved numerically, using the analytical expression for $\mu^2(\gamma, \beta)$ given in eq. (G5) of ref. [7].

## B Use of the astrophysical catalogs

We extract the GW luminosity in a representative comoving volume element of the universe following the merger tree philosophy of refs. [29, 30], for the three different scenarios outlined in section 3.2.

The merger tree algorithms of refs. [29, 30] provide us lists of synthetic mergers events, with each event $i$ being associated to a merger redshift $z_m^i$ and a chirp mass $M^i$. Each event is also attributed a weight $W^i$ (with dimension of inverse volume) which gives the catalog the interpretation of a box of comoving volume $V_c^{\text{box}} = 1 \text{ Mpc}^3$ whose content evolves over $z$. The weighting of events by $W^i$ follows from a Press-Schechter prescription, we refer to refs. [29, 30] for more details.

Assuming that the catalog binaries are emitting GW in the PTA frequency band following a quasi-circular relative orbital motion and a standard cosmology, the pair $(z_m^i, M^i)$ fully determines $f_s^i(z)$ and $dE_s^i(z)/dt_s$, that is, the whole history (over redshift $z$) of the frequency and power emission of source $i$, via the chirp equation [36]:

$$ \frac{d f_s^i}{dt_s} = \frac{96}{5} \pi^{8/3} \left( \frac{GM^i}{c^3} \right)^{5/3} (f_s^i)^{11/3}. \qquad (B.1) $$

We associate the local time of the source at redshift $z$ to be the cosmic time $t_s(z)$ of the Universe.

The Press-Schechter weights ensure that the box comoving volume at redshift $z$ is a fair representation of the Universe at that redshift. The sky-averaged GW luminosity density of the Universe is therefore assumed to be equal to the one of the box, which is given by the cumulative power of all catalog sources $i$ at redshift $z$,

$$ \frac{d\bar{L}_{\text{GW}}}{dV_c}(z) = \frac{L_{\text{GW}}^{\text{box}}(z)}{V_c^{\text{box}}} = \sum_i W_i \frac{dE_s^i}{dt_s}(z), \qquad (B.2) $$

containing the contribution of all sources regardless of their frequencies. This function may in practice not be very smooth, since, at each redshift $z$, it consists of a sum over the contributions of a discrete and possibly small amount of sources. We therefore smooth it by taking a (sliding) average at each redshift $z$. In the following, $\Delta t_s$ is the source time interval over which the smoothing average is taken:

$$ \sum_i W_i \frac{dE_s^i}{dt_s} \approx \frac{1}{\Delta t_s} \int_{\Delta t_s} dt_s \sum_i W_i \frac{dE_s^i}{dt_s} = \frac{1}{\Delta t_s} \sum_i W_i \int_{\Delta t_s} dt_s \frac{d f_s^i}{d f_s^i} \qquad (B.3) $$

$$ \begin{aligned} &= \frac{\pi^{2/3}}{3G} \frac{1}{\Delta t_s} \sum_i W_i (GM^i)^{5/3} \int_{\Delta f_s^i} \frac{d f_s^i (f_s^i)^{-1/3}}{d f_s^i} \\ &= \frac{\pi^{2/3}}{2G} \frac{1}{\Delta t_s} \sum_i W_i (GM^i)^{5/3} \left[ (f_s^i)^{2/3} \right]_{\min(\Delta f_s^i)}^{\max(\Delta f_s^i)}, \end{aligned} \qquad (B.4) $$
where, for each $i$, we changed the integration variable from the cosmic time at the source $t_s$ to $f_s^i$. Moreover, $\Delta f_s^i$ is the source frequency interval of integration corresponding to the so far generic smoothing time interval $\Delta t_s$. In practice, we use a smoothing cosmic time interval at the source such that it corresponds to a redshift interval $\Delta z = 0.1$, that is, $\Delta t_s(z) = |t_s(z + 0.05) - t_s(z - 0.05)|$. We then refine the analysis with a logarithmic $\Delta z$ smoothing at $z < 0.05$ to capture, with a finer resolution, the most recent catalog events, since the luminosity must vanish at $z = 0$ (no source is present directly at the observer).

The above luminosity contains the contributions of all frequencies. The analysis is then divided into the 14 different observed frequency bands of the NANOGrav collaboration [1], that is, $f_j = j/T$, $j = 1, \dots, 14$, with $\Delta f = 1/T$ and $T$ the observation time. For each catalog event $i$, we numerically solve the chirp equation (B.1) to define

$$ \Delta z^{i,j} = \left\{ z : f_s^i(z)/(1+z) \in [f_j - \Delta f/2, f_j + \Delta f/2] \right\}, \qquad (B.5) $$

as the redshift interval over which source $i$ emits at a frequency that is observed in band $j$.

We then associate a source frequency interval $\Delta f_s^{i,j}$ to $\Delta z^{i,j}$ in the same way as we have previously associated a frequency interval $\Delta f_s^i$ to a generic source time interval $\Delta t_s$ or redshift interval $\Delta z$. If we then substitute

$$ \left[ (f_s^i)^{2/3} \right]_{\min(\Delta f_s^i)}^{\max(\Delta f_s^i)} \qquad (B.6) $$

by

$$ \sum_{j=1}^{14} \left[ (f_s^i)^{2/3} \right]_{\min(\Delta f_s^i \cap \Delta f_s^{i,j})}^{\max(\Delta f_s^i \cap \Delta f_s^{i,j})}, \qquad (B.7) $$

in eq. (B.4), we restrict the luminosity to the frequency contributions in the 14 observed PTA bands. The restriction to a single $j$-term in the sum (B.7) defines the luminosity in the observed band $j$.

## References

[1] NANOGRAV collaboration, *The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background*, Astrophys. J. Lett. 951 (2023) L8 [arXiv:2306.16213] [INSPIRE].

[2] D.J. Reardon et al., *Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array*, Astrophys. J. Lett. 951 (2023) L6 [arXiv:2306.16215] [INSPIRE].

[3] EPTA collaboration, *The second data release from the European Pulsar Timing Array — I. The dataset and timing analysis*, Astron. Astrophys. 678 (2023) A48 [arXiv:2306.16224] [INSPIRE].

[4] H. Xu et al., *Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I*, Res. Astron. Astrophys. 23 (2023) 075024 [arXiv:2306.16216] [INSPIRE].

[5] R. Hellings and G. Downs, *Upper limits on the isotropic gravitational radiation background from pulsar timing analysis*, Astrophys. J. Lett. 265 (1983) L39 [INSPIRE].

[6] B. Allen, *Will pulsar timing arrays observe the hellings and downs correlation curve?*, Frascati Phys. Ser. 74 (2022) 65 [INSPIRE].
[7] B. Allen, *Variance of the Hellings-Downs correlation*, Phys. Rev. D **107** (2023) 043018
[arXiv:2205.05637] [INSPIRE].

[8] B. Allen and J.D. Romano, *Hellings and Downs correlation of an arbitrary set of pulsars*, Phys. Rev. D **108** (2023) 043026 [arXiv:2208.07230] [INSPIRE].

[9] J.D. Romano and B. Allen, *Answers to frequently asked questions about the pulsar timing array Hellings and Downs curve*, Class. Quant. Grav. **41** (2024) 175008 [arXiv:2308.05847] [INSPIRE].

[10] B. Allen and S. Valtolina, *Pulsar timing array source ensembles*, Phys. Rev. D **109** (2024) 083038 [arXiv:2401.14329] [INSPIRE].

[11] R.C. Bernardo and K.-W. Ng, *Pulsar and cosmic variances of pulsar timing-array correlation measurements of the stochastic gravitational wave background*, JCAP **11** (2022) 046 [arXiv:2209.14834] [INSPIRE].

[12] N. Grimm, M. Pijnenburg, G. Cusin and C. Bonvin, *The impact of large-scale galaxy clustering on the variance of the Hellings-Downs correlation: theoretical framework*, JCAP **03** (2025) 011 [arXiv:2404.05670] [INSPIRE].

[13] B. Allen and A.C. Ottewill, *Detection of anisotropies in the gravitational wave stochastic background*, Phys. Rev. D **56** (1997) 545 [gr-qc/9607068] [INSPIRE].

[14] G. Cusin, C. Pitrou and J.-P. Uzan, *Anisotropy of the astrophysical gravitational wave background: Analytic expression of the angular power spectrum and correlation with cosmological observations*, Phys. Rev. D **96** (2017) 103019 [arXiv:1704.06184] [INSPIRE].

[15] G. Cusin, C. Pitrou and J.-P. Uzan, *The signal of the gravitational wave background and the angular correlation of its energy density*, Phys. Rev. D **97** (2018) 123527 [arXiv:1711.11345] [INSPIRE].

[16] G. Cusin, I. Dvorkin, C. Pitrou and J.-P. Uzan, *First predictions of the angular power spectrum of the astrophysical gravitational wave background*, Phys. Rev. Lett. **120** (2018) 231101 [arXiv:1803.03236] [INSPIRE].

[17] G. Cusin, R. Durrer and P.G. Ferreira, *Polarization of a stochastic gravitational wave background through diffusion by massive structures*, Phys. Rev. D **99** (2019) 023534 [arXiv:1807.10620] [INSPIRE].

[18] G. Cusin, I. Dvorkin, C. Pitrou and J.-P. Uzan, *Properties of the stochastic astrophysical gravitational wave background: astrophysical sources dependencies*, Phys. Rev. D **100** (2019) 063004 [arXiv:1904.07797] [INSPIRE].

[19] G. Cusin, I. Dvorkin, C. Pitrou and J.-P. Uzan, *Stochastic gravitational wave background anisotropies in the mHz band: astrophysical dependencies*, Mon. Not. Roy. Astron. Soc. **493** (2020) L1 [arXiv:1904.07757] [INSPIRE].

[20] C. Pitrou, G. Cusin and J.-P. Uzan, *Unified view of anisotropies in the astrophysical gravitational-wave background*, Phys. Rev. D **101** (2020) 081301 [arXiv:1910.04645] [INSPIRE].

[21] A.C. Jenkins, J.D. Romano and M. Sakellariadou, *Estimating the angular power spectrum of the gravitational-wave background in the presence of shot noise*, Phys. Rev. D **100** (2019) 083501 [arXiv:1907.06642] [INSPIRE].

[22] D. Alonso, G. Cusin, P.G. Ferreira and C. Pitrou, *Detecting the anisotropic astrophysical gravitational wave background in the presence of shot noise through cross-correlations*, Phys. Rev. D **102** (2020) 023002 [arXiv:2002.02888] [INSPIRE].
[23] A.I. Renzini, J.D. Romano, C.R. Contaldi and N.J. Cornish, Comparison of maximum-likelihood mapping methods for gravitational-wave backgrounds, Phys. Rev. D 105 (2022) 023519 [arXiv:2107.02292] [INSPIRE].

[24] C.M.F. Mingarelli, T. Sidery, I. Mandel and A. Vecchio, Characterizing gravitational wave stochastic background anisotropy with pulsar timing arrays, Phys. Rev. D 88 (2013) 062005 [arXiv:1306.5394] [INSPIRE].

[25] S.R. Taylor and J.R. Gair, Searching For Anisotropic Gravitational-wave Backgrounds Using Pulsar Timing Arrays, Phys. Rev. D 88 (2013) 084001 [arXiv:1306.5395] [INSPIRE].

[26] J. Gair, J.D. Romano, S. Taylor and C.M.F. Mingarelli, Mapping gravitational-wave backgrounds using methods from CMB analysis: Application to pulsar timing arrays, Phys. Rev. D 90 (2014) 082001 [arXiv:1406.4664] [INSPIRE].

[27] Y. Ali-Haïmoud and M. Kamionkowski, Cosmic microwave background limits on accreting primordial black holes, Phys. Rev. D 95 (2017) 043534 [arXiv:1612.05644] [INSPIRE].

[28] N.M.J. Cruz, A. Malhotra, G. Tasinato and I. Zavala, Measuring kinematic anisotropies with pulsar timing arrays, Phys. Rev. D 110 (2024) 063526 [arXiv:2402.17312] [INSPIRE].

[29] E. Barausse, The evolution of massive black holes and their spins in their galactic hosts, Mon. Not. Roy. Astron. Soc. 423 (2012) 2533 [arXiv:1201.5888] [INSPIRE].

[30] A. Klein et al., Science with the space-based interferometer eLISA: Supermassive black hole binaries, Phys. Rev. D 93 (2016) 024003 [arXiv:1511.05581] [INSPIRE].

[31] C. Pitrou and G. Cusin, Mitigating cosmic variance in the Hellings-Downs curve: a Cosmic Microwave Background analogy, arXiv:2412.12073 [INSPIRE].

[32] B. Allen and J.D. Romano, Optimal Reconstruction of the Hellings and Downs Correlation, Phys. Rev. Lett. 134 (2025) 031401 [arXiv:2407.10968] [INSPIRE].

[33] B. Allen, D. Agarwal, J.D. Romano and S. Valtolina, Source anisotropies and pulsar timing arrays, Phys. Rev. D 110 (2024) 123507 [arXiv:2406.16031] [INSPIRE].

[34] B. Allen, Pulsar timing array harmonic analysis and source angular correlations, Phys. Rev. D 110 (2024) 043043 [arXiv:2404.05677] [INSPIRE].

[35] D. Agarwal and J.D. Romano, Cosmic variance of the Hellings and Downs correlation for ensembles of universes having nonzero angular power spectra, Phys. Rev. D 110 (2024) 043044 [arXiv:2404.08574] [INSPIRE].

[36] M. Maggiore, Gravitational Waves. Vol. 1: Theory and Experiments, Oxford Master Series in Physics, Oxford University Press (2007).

[37] I. Dvorkin, J.-P. Uzan, E. Vangioni and J. Silk, Synthetic model of the gravitational wave background from evolving binary compact objects, Phys. Rev. D 94 (2016) 103011 [arXiv:1607.06818] [INSPIRE].

[38] E.S. Phinney, A practical theorem on gravitational wave backgrounds, astro-ph/0108028 [INSPIRE].

[39] C. Bonvin and R. Durrer, What galaxy surveys really measure, Phys. Rev. D 84 (2011) 063505 [arXiv:1105.5280] [INSPIRE].

[40] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology, Oxford University Press (2018).
[41] A. Sesana, A. Vecchio and C.N. Colacino, The stochastic gravitational-wave background from massive black hole binary systems: implications for observations with Pulsar Timing Arrays, Mon. Not. Roy. Astron. Soc. 390 (2008) 192 [arXiv:0804.4476] [INSPIRE].

[42] W.L.W. Sargent et al., Dynamical evidence for a central mass concentration in the galaxy M87, Astrophys. J. 221 (1978) 731 [INSPIRE].

[43] R.D. Blandford and C.F. McKee, Reverberation mapping of the emission line regions of Seyfert galaxies and quasars, Astrophys. J. 255 (1982) 419.

[44] B.M. Peterson, Reverberation mapping of active galactic nuclei, Publ. Astron. Soc. Pac. 685 (1993) 247 [INSPIRE].

[45] SDSS collaboration, A survey of $z > 5.8$ quasars in the Sloan Digital Sky Survey I: Discovery of three new quasars and the spatial density of luminous quasars at $z \sim 6$, Astron. J. 122 (2001) 2833 [astro-ph/0108063] [INSPIRE].

[46] L. Ferrarese and D. Merritt, A fundamental relation between supermassive black holes and their host galaxies, Astrophys. J. Lett. 539 (2000) L9 [astro-ph/0006053] [INSPIRE].

[47] K. Inayoshi, E. Visbal and Z. Haiman, The Assembly of the First Massive Black Holes, Ann. Rev. Astron. Astrophys. 58 (2020) 27 [arXiv:1911.05791] [INSPIRE].

[48] M. Volonteri, M. Habouzit and M. Colpi, The origins of massive black holes, Nature Rev. Phys. 3 (2021) 732 [arXiv:2110.10175] [INSPIRE].

[49] E.L. Turner, Quasars and galaxy formation. i-the z greater than 4 objects, Astron. J. 101 (1991) 5.

[50] C. Lacey and S. Cole, Merger rates in hierarchical models of galaxy formation, Mon. Not. Roy. Astron. Soc. 262 (1993) 627 [INSPIRE].

[51] O. Fakhouri, C.-P. Ma and M. Boylan-Kolchin, The Merger Rates and Mass Assembly Histories of Dark Matter Haloes in the Two Millennium Simulations, Mon. Not. Roy. Astron. Soc. 406 (2010) 2267 [arXiv:1001.2304] [INSPIRE].

[52] V. Desjacques, D. Jeong and F. Schmidt, Large-Scale Galaxy Bias, Phys. Rept. 733 (2018) 1 [arXiv:1611.09787] [INSPIRE].

[53] L. Ferrarese, Beyond the bulge: a fundamental relation between supermassive black holes and dark matter halos, Astrophys. J. 578 (2002) 90 [astro-ph/0203469] [INSPIRE].

[54] J. Colin, R. Mohayaee, M. Rameez and S. Sarkar, High redshift radio galaxies and divergence from the CMB dipole, Mon. Not. Roy. Astron. Soc. 471 (2017) 1045 [arXiv:1703.09376] [INSPIRE].

[55] C.A.P. Bengaly, R. Maartens and M.G. Santos, Probing the Cosmological Principle in the counts of radio galaxies at different frequencies, JCAP 04 (2018) 031 [arXiv:1710.08804] [INSPIRE].

[56] N.J. Secrest et al., A Test of the Cosmological Principle with Quasars, Astrophys. J. Lett. 908 (2021) L51 [arXiv:2009.14826] [INSPIRE].

[57] T.M. Siewert, M. Schmidt-Rubart and D.J. Schwarz, Cosmic radio dipole: Estimators and frequency dependence, Astron. Astrophys. 653 (2021) A9 [arXiv:2010.08366] [INSPIRE].

[58] N.J. Secrest et al., A Challenge to the Standard Cosmological Model, Astrophys. J. Lett. 937 (2022) L31 [arXiv:2206.05624] [INSPIRE].

[59] N. Grimm et al., Combining chirp mass, luminosity distance, and sky localization from gravitational wave events to detect the cosmic dipole, Mon. Not. Roy. Astron. Soc. 526 (2023) 4673 [arXiv:2309.00336] [INSPIRE].
[60] S. Mastrogiovanni, C. Bonvin, G. Cusin and S. Foffa, Detection and estimation of the cosmic dipole with the einstein telescope and cosmic explorer, Mon. Not. Roy. Astron. Soc. 521 (2023) 984 [arXiv:2209.11658] [INSPIRE].

[61] B. Cousins, A. Dhani, B.S. Sathyaprakash and N. Yunes, Finding cosmic anisotropy with networks of next-generation gravitational-wave detectors, arXiv:2406.15550 [INSPIRE].

[62] D.N. Limber, The analysis of counts of the extragalactic nebulae in terms of a fluctuating density field, Astrophys. J. 117 (1953) 134.

[63] M. Bartelmann and P. Schneider, Weak gravitational lensing, Phys. Rept. 340 (2001) 291 [astro-ph/9912508] [INSPIRE].

[64] M. LoVerde and N. Afshordi, Extended Limber Approximation, Phys. Rev. D 78 (2008) 123506 [arXiv:0809.5112] [INSPIRE].

[65] M. Martinelli et al., Ultralarge-scale approximations and galaxy clustering: Debiasing constraints on cosmological parameters, Mon. Not. Roy. Astron. Soc. 510 (2022) 1964 [arXiv:2106.15604] [INSPIRE].

[66] J. Lesgourgues, The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview, arXiv:1104.2932 [INSPIRE].

[67] Y. Zu et al., Cosmological Constraints from the Large Scale Weak Lensing of SDSS MaxBCG Clusters, Mon. Not. Roy. Astron. Soc. 439 (2014) 1628 [arXiv:1207.3794] [INSPIRE].

[68] A.J. Mead and L. Verde, Including beyond-linear halo bias in halo models, Mon. Not. Roy. Astron. Soc. 503 (2021) 3095 [arXiv:2011.08858] [INSPIRE].

[69] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[70] S. Babak, M. Falxa, G. Franciolini and M. Pieroni, Forecasting the sensitivity of pulsar timing arrays to gravitational wave backgrounds, Phys. Rev. D 110 (2024) 063022 [arXiv:2404.02864] [INSPIRE].

[71] P.F. Depta, V. Domcke, G. Franciolini and M. Pieroni, Pulsar timing array sensitivity to anisotropies in the gravitational wave background, arXiv:2407.14460 [INSPIRE].

[72] F. Semenzato et al., Cross-Correlating the Universe: The Gravitational Wave Background and Large-Scale Structure, arXiv:2411.00532 [INSPIRE].