The impact of large-scale galaxy clustering on the variance of the Hellings-Downs correlation: theoretical framework

Recently, pulsar timing arrays (PTAs) were used to deliver evidence for the existence of a stochastic gravitational wave (GW) background NANOGrav:2023gor ; Reardon:2023gzh ; EPTA:2023sfo ; Xu:2023wog . This detection utilizes the angular correlation of timing residuals, i.e. small perturbations in the arrival times of radio pulses caused by the passing of GWs. The origin of the GW background is still unclear, and it might be of both astrophysical or cosmological nature (see e.g. Refs. Shannon:2015ect ; Chen:2019xse ; NANOGrav:2021flc ; NANOGrav:2023hvm ; Vagnozzi:2023lwo ). One possibility to distinguish between these scenarios is to evaluate whether galaxy clustering has an impact on the angular correlation of timing residuals: if sources contributing to the GW background are of astrophysical origin, we expect their distribution to follow the distribution of galaxies. Galaxies are however not perfectly isotropically distributed on the sky: they are grouped in clusters, that are themselves arranged into a large-scale structure. Anisotropies in the galaxy distribution have been observed by various surveys, leading to correlations between galaxy positions up to very large cosmological scales, see e.g. Ref. Blake_2011 ; Howlett:2014opa ; Pezzotta:2016gbo ; Alam2016:1607.03155v1 ; eBOSS:2020yzd . Therefore, any GW background of astrophysical origin is expected to be as well anisotropic: an intrinsic source of anisotropy is due to clustering (the distribution of GW sources follows the one of the large-scale structure), a secondary one is due to the fact that once a GW signal is emitted, it feels the effect of the gravitational potential of structures Cusin:2017fwz ; Cusin:2017mjm ; Cusin:2018rsq . These secondary anisotropies are expected to be subdominant on all scales Pitrou:2019rjz , hence we will focus here on the effect of large-scale galaxy clustering.

Usually, the correlation of pulsar timing residuals due to the passing of GWs, as a function of the pulsar pair separation on the sky, is assumed to be described by the Hellings and Downs (HD) curve Hellings:1983fr . Recent literature has investigated to which extent this curve is precise in the description of correlations of timing residuals Allen:2022dzg ; Allen:2022ksj ; Romano:2023zhb ; Allen:2024rqk ; Bernardo:2022xzl . It has been found that, in practice, PTA observations are not expected to exactly recover the HD curve in any realistic models of the GW background. For instance, the interference of GW sources would cause a departure from this idealized curve Allen:2022bjz ; Allen:2022dzg . In all these studies, the fact that sources are clustered is systematically neglected.¹¹1The possibility of anisotropies in the distribution of GW sources in the PTA band has been pointed out in previous work, however either not in the context of the HD correlation Allen:1996gp ; Cusin:2017fwz ; Cusin:2017mjm ; Cusin:2018rsq ; Cusin:2018avf ; Cusin:2019jpv ; Cusin:2019jhg ; Pitrou:2019rjz ; Jenkins:2019nks ; Alonso:2020mva ; renzini2022 , or without targeting the impact on its variance Mingarelli:2013dsa ; Taylor:2013esa ; Gair:2014rwa ; Ali-Haimoud:2016mbv .

In this work, we present, for the first time, a framework to account for the impact of the cosmological large-scale structure on the HD correlation.²²2Note that in this manuscript, we will use the notion HD correlation to denote the observed correlation of pulsar timing delays, while HD curve denotes the theoretical prediction made by Hellings & Downs in 1983 Hellings:1983fr . This method, novel in the context of PTAs, is based on a two-step averaging approach: first, we describe a Gaussian ensemble of GW sources in a realization of the universe with a given distribution of galaxies. Second, to determine the impact of the galaxy distribution on the HD correlation, we study how it varies not only from ensemble to ensemble of GW sources, but also from realization to realization of the galaxy distribution. Mathematically, the latter is done by taking a second ensemble average over all such realizations. We find that the mean of the HD correlation itself is not affected, which means that anisotropies in the distribution of sources do not generate a systematic bias in the HD correlation. In other words, if we were able to measure the correlations of pulsar timing residuals in many realizations of the universe, we would in average recover the HD curve predicted in Ref. Hellings:1983fr for the case of isotropically distributed sources. However, in practice we have only access to one such realization. To account for this, we calculate the variance of the HD correlation due to the large-scale structure of the Universe, which tells us by how much the measurements in our specific realization of the universe can differ from the mean. In this paper we find that, indeed, large-scale galaxy clustering leads to a new contribution to the HD variance, which needs to be carefully evaluated in light of future observations. Numerical results for specific population models of supermassive black holes are presented in a companion paper Grimm:2024hgi .

In addition to the variance due the clustering of sources, there is as well a variance due to shot noise, i.e. due to the fact that GW sources are not a continuous field, but rather have a discrete Poisson distribution (see e.g. Ref. Jenkins:2019nks ). In this paper we do not take the shot noise contribution into account (see Ref. Allen:2024mtn for a derivation of this contribution for PTAs), and we concentrate only on the modeling of the clustering variance. Consequently, we refer to an isotropic distribution of sources when no clustering is present, and an anisotropic distribution of sources when clustering is accounted for. This use of terminology, although being common in literature, is not fully accurate due to shot noise: since sources have a discrete nature, they will never have a perfectly isotropic distribution, even in the absence of clustering.

The remainder of this work is structured as follows: In section II, we first revise basic notations in the case of an isotropic distribution of GW sources, and then introduce galaxy clustering and discuss the respective ensemble averaging procedure. We present our calculation of the variance of the HD correlation due to large-scale clustering in section III, with our main result given in Eq. (31). We conclude in section IV, and comment on the positivity of this novel contribution to the variance in Appendix A.

Before calculating the impact of large-scale structure on the HD correlation, we first determine how it impacts the two-point function of the GW amplitude. In the following, we first discuss the different layers of stochasticity present in our analyses, then discuss the GW signal for an isotropic distribution of sources, and subsequently introduce the impact of galaxy clustering.

The HD correlation describes the correlation of timing residuals between pulsars, due to the passing of GWs. Following Ref. Allen:2022dzg , we calculate the redshift, $Z_{a}(t)$, induced by the GW on a pulse emitted by pulsar $a$ located in direction ${\bf n}_{a}$ and observed at time $t$ on Earth:

$$Z_{a}(t)\equiv\sum_{A}\int{\rm d}f\int{\rm d}\mathbf{n}\,F_{a}^{A}(\mathbf{n})h_{A}(f,\mathbf{n}){\rm e}^{2\pi ift}\left(1-{\rm e}^{-2\pi if\tau_{a}(1+\mathbf{n}\cdot{\bf n}_{a})}\right)\,.$$

(12)

Here, $\tau_{a}$ is the light travel time between emission and observation, and the antenna pattern functions $F_{a}^{A}(\mathbf{n})$ are defined as

$$F_{a}^{A}(\mathbf{n})\equiv\frac{n_{a}^{i}n_{a}^{j}e_{ij}^{A}(\mathbf{n})}{2(1+\mathbf{n}\cdot\mathbf{n}_{a})}\,.$$

(13)

Note that the property $h_{A}(-f,\mathbf{n})=h_{A}^{\ast}(f,\mathbf{n})$ implies that $Z_{a}(t)\in\mathbb{R}$. Since the redshift of a pulsar is related to its timing residual by time differentiation, the goal of PTA observations can be formulated as the search for correlations between the redshifts of pulsars. We define therefore $\rho_{ab}$ as the product of the redshifts of two pulsar signals averaged over the observation time $T$,⁸⁸8In this paper, $\rho_{ab}$ will always be used as defined in Eq. (14), and should not be confused with a density.

	$\displaystyle\rho_{ab}$	$\displaystyle\equiv\overline{Z_{a}Z_{b}}=\overline{Z_{a}Z_{b}^{\ast}}=\frac{1}{T}\int_{-T/2}^{T/2}{\mathrm{d}}t\ Z_{a}Z_{b}^{\ast}$
		$\displaystyle=\sum_{A,A^{\prime}}\int\mathrm{d}f\int\mathrm{d}f^{\prime}\int\mathrm{d}\mathbf{n}\int\mathrm{d}\mathbf{n}^{\prime}R^{A}_{a}(f,\mathbf{n})^{\ast}R_{b}^{A^{\prime}}(f^{\prime},\mathbf{n}^{\prime})h_{A}^{\ast}(f,\mathbf{n})h_{A^{\prime}}(f^{\prime},\mathbf{n}^{\prime})\mbox{sinc}\left(\pi(f-f^{\prime})T\right)\,,$		(14)

where we followed Ref. Allen:2022dzg in the use of the shorthand notation,

$$R^{A}_{a}(f,\mathbf{n})\equiv\left(1-{\rm e}^{-2\pi if\tau_{a}(1+\mathbf{n}\cdot\mathbf{n}_{a})}\right)F^{A}_{a}(\mathbf{n})\,.$$

(15)

The function $\mbox{sinc}(x)\equiv\sin(x)/x$ in Eq. (14) arises from the time averaging process:

$$\frac{1}{T}\int_{-T/2}^{T/2}{\rm d}t\ {\rm e}^{-2\pi ift}{\rm e}^{2\pi if^{\prime}t}=\frac{\sin\left(\pi(f-f^{\prime})T\right)}{\pi(f-f^{\prime})T}\,.$$

(16)

In the standard case of isotropically distributed sources, the expectation value $\left\langle\rho_{ab}\right\rangle$ over field realizations of $h_{A}$ leads to the standard HD angular correlation, as is done e.g. in Ref. Allen:2022dzg . When sources are anisotropically distributed, however, the expectation value of $\rho_{ab}$ depends on the specific realization through the spectral density $S_{h}(f,\mathbf{n})$, as we see from Eq. (7). We can compute the ensemble average of $\rho_{ab}$ over all realizations of the galaxy density, and using Eq. (11), we see immediately that we recover the standard HD curve predicted in Ref. Hellings:1983fr . This is not surprising, since the ensemble average of the galaxy density is itself homogeneous and isotropic. For future reference, we write this property of $\rho_{ab}$ as

$$\left\langle\rho_{ab}\right\rangle=\left\langle\rho_{ab}\right\rangle^{\text{st.}}\,.$$

(17)

In other words, even when taking anisotropies into account, $\left\langle\rho_{ab}\right\rangle$ takes the “standard” form $\left\langle\rho_{ab}\right\rangle^{\text{st.}}$, namely the one that it has in computations assuming isotropy. The relevant question is to determine to which extent the observed HD correlation in our specific realization of the universe may differ from the expectation value. To assess this, we need to compute the variance of $\rho_{ab}$. Before proceeding with the calculation, we note that $\rho_{ab}$ is defined involving only one pair of pulsars. In practice, PTA surveys track many pulsars forming a large number of pairs, which needs to be taken into account when evaluating the variance for a realistic scenario Allen:2022ksj ; Allen:2022dzg ; Romano:2023zhb . Indeed, the concepts presented in the previous sections are valid independently of how many pulsars are considered. Here, we demonstrate how our formalism is applied in the case of a single pulsar pair. In our companion paper, Ref. Grimm:2024hgi , we present results for the observational limits of a single as well as infinitely many pulsar pairs.

From Eq. (14) we have

	$\displaystyle\rho^{2}_{ab}=\sum_{A,A^{\prime},A^{\prime\prime},A^{\prime\prime\prime}}\iiiint\mathrm{d}f\mathrm{d}f^{\prime}\mathrm{d}f^{\prime\prime}\mathrm{d}f^{\prime\prime\prime}\iiiint\mathrm{d}\mathbf{n}\,\mathrm{d}\mathbf{n}^{\prime}\mathrm{d}\mathbf{n}^{\prime\prime}\mathrm{d}\mathbf{n}^{\prime\prime\prime}\mbox{sinc}\left(\pi(f-f^{\prime})T\right)\mbox{sinc}\left(\pi(f^{\prime\prime}-f^{\prime\prime\prime})T\right)$
	$\displaystyle\times R^{A}_{a}(f,\mathbf{n})^{\ast}R^{A^{\prime}}_{b}(f^{\prime},\mathbf{n}^{\prime})R^{A^{\prime\prime}}_{a}(f^{\prime\prime},\mathbf{n}^{\prime\prime})R^{A^{\prime\prime\prime}}_{b}(f^{\prime\prime\prime},\mathbf{n}^{\prime\prime\prime})^{\ast}h_{A}^{\ast}(f,\mathbf{n})h_{A^{\prime}}(f^{\prime},\mathbf{n}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},\mathbf{n}^{\prime\prime})h_{A^{\prime\prime\prime}}^{\ast}(f^{\prime\prime\prime},\mathbf{n}^{\prime\prime\prime})\,.$		(18)

The variance is given by

$$\left\langle\rho^{2}_{ab}\right\rangle-\left\langle\rho_{ab}\right\rangle^{2}=\left\langle\left\langle\rho^{2}_{ab}\right\rangle_{h|\rho_{g}}\right\rangle_{\mathfrak{p}_{g}}-\left\langle\left\langle\rho_{ab}\right\rangle_{h|\rho_{g}}\right\rangle^{2}_{\mathfrak{p}_{g}}\,.$$

(19)

For the inner ensemble average we use Eq. (7), and the fact that the field $h_{A}$ is Gaussian given any fixed realization of the density field, as discussed in section II.4. We can thus apply Isserlis’ theorem (also known as Wick’s theorem) to the conditioned field. This theorem implies that for a Gaussian field, correlations of four fields can be expressed as sums of products of correlations of two fields. Our calculation is following Appendix C of Ref. Allen:2022dzg with the extension that now, we need to take the anisotropic term of $S_{h}(f,\mathbf{n})$ in Eq. (8) into account. More specifically, Isserlis’ theorem reads

	$\displaystyle\left\langle h_{A}^{\ast}(f,\mathbf{n})h_{A^{\prime}}(f^{\prime},\mathbf{n}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},\mathbf{n}^{\prime\prime})h_{A^{\prime\prime\prime}}^{\ast}(f^{\prime\prime\prime},\mathbf{n}^{\prime\prime\prime})\right\rangle{{}_{h|\rho_{g}}}=$
	$\displaystyle\left\langle h_{A}^{\ast}(f,\mathbf{n})h_{A^{\prime}}(f^{\prime},\mathbf{n}^{\prime})\right\rangle{{}_{h|\rho_{g}}}\left\langle h_{A^{\prime\prime}}^{\ast}(-f^{\prime\prime},\mathbf{n}^{\prime\prime})h_{A^{\prime\prime\prime}}(-f^{\prime\prime\prime},\mathbf{n}^{\prime\prime\prime})\right\rangle{{}_{h|\rho_{g}}}$
	$\displaystyle+\left\langle h_{A}^{\ast}(f,\mathbf{n})h_{A^{\prime\prime}}(f^{\prime\prime},\mathbf{n}^{\prime\prime})\right\rangle{{}_{h|\rho_{g}}}\left\langle h_{A^{\prime}}^{\ast}(-f^{\prime},\mathbf{n}^{\prime})h_{A^{\prime\prime\prime}}(-f^{\prime\prime\prime},\mathbf{n}^{\prime\prime\prime})\right\rangle{{}_{h|\rho_{g}}}$
	$\displaystyle+\left\langle h_{A}^{\ast}(f,\mathbf{n})h_{A^{\prime\prime\prime}}(-f^{\prime\prime\prime},\mathbf{n}^{\prime\prime\prime})\right\rangle{{}_{h|\rho_{g}}}\left\langle h_{A^{\prime}}^{\ast}(-f^{\prime},\mathbf{n}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},\mathbf{n}^{\prime\prime})\right\rangle{{}_{h|\rho_{g}}}\,,$		(20)

where negative signs in the arguments appear making use of the relation $h_{A}(f,\mathbf{n})=h_{A}^{\ast}(-f,\mathbf{n})$. Applying Eq. (7), we find that

	$\displaystyle\left\langle h_{A}^{\ast}(f,\mathbf{n})h_{A^{\prime}}(f^{\prime},\mathbf{n}^{\prime})h_{A^{\prime\prime}}(f^{\prime\prime},\mathbf{n}^{\prime\prime})h_{A^{\prime\prime\prime}}^{\ast}(f^{\prime\prime\prime},\mathbf{n}^{\prime\prime\prime})\right\rangle{{}_{h|\rho_{g}}}=$
	$\displaystyle\frac{1}{4}\delta_{AA^{\prime}}\delta_{A^{\prime\prime}A^{\prime\prime\prime}}\frac{\delta^{2}(\mathbf{n},\mathbf{n}^{\prime})}{4\pi}\frac{\delta^{2}(\mathbf{n}^{\prime\prime},\mathbf{n}^{\prime\prime\prime})}{4\pi}\delta(f-f^{\prime})\delta(f^{\prime\prime}-f^{\prime\prime\prime})S_{h}(f,\mathbf{n})S_{h}(f^{\prime\prime},\mathbf{n}^{\prime\prime})$
	$\displaystyle+\frac{1}{4}\delta_{AA^{\prime\prime}}\delta_{A^{\prime}A^{\prime\prime\prime}}\frac{\delta^{2}(\mathbf{n},\mathbf{n}^{\prime\prime})}{4\pi}\frac{\delta^{2}(\mathbf{n}^{\prime},\mathbf{n}^{\prime\prime\prime})}{4\pi}\delta(f-f^{\prime\prime})\delta(f^{\prime}-f^{\prime\prime\prime})S_{h}(f,\mathbf{n})S_{h}(f^{\prime},\mathbf{n}^{\prime})$
	$\displaystyle+\frac{1}{4}\delta_{AA^{\prime\prime\prime}}\delta_{A^{\prime}A^{\prime\prime}}\frac{\delta^{2}(\mathbf{n},\mathbf{n}^{\prime\prime\prime})}{4\pi}\frac{\delta^{2}(\mathbf{n}^{\prime},\mathbf{n}^{\prime\prime})}{4\pi}\delta(f+f^{\prime\prime\prime})\delta(f^{\prime}+f^{\prime\prime})S_{h}(f,\mathbf{n})S_{h}(f^{\prime},\mathbf{n}^{\prime})\,.$		(21)

Whenever products of the type $R^{A\ast}_{a}R^{A}_{b}$ are considered in integrals, it is common to use the approximation that $\tau_{a,b}f\gg 1$ for situations of practical interest, and thus functions with unequal and highly oscillating phases average out:

$$F^{A}_{a}F^{A}_{b}\ \left[1-{\rm e}^{2\pi if\tau_{a}(1+\mathbf{n}\cdot\mathbf{n}_{a})}\right]\left[1-{\rm e}^{-2\pi if\tau_{b}(1+\mathbf{n}\cdot\mathbf{n}_{b})}\right]\simeq F^{A}_{a}F^{A}_{b}\ (1+\delta_{ab})\,.$$

(22)

This approximation is equivalent to considering only the correlation of the Earth-term for two distinct pulsars. The contribution of the pulsar term is only relevant for the auto-correlation, i.e. with identical pulsars $a=b$, and falls off fast for non-zero angular or radial separations (see e.g. Ref. Mingarelli:2014xfa or Appendix C.2 of Ref. Allen:2022dzg for the contributions of pulsar terms). Thus, this term can be neglected when $a\neq b$. The auto-correlation ($a=b$) term, on the other hand, carries contributions from the Earth and pulsar terms, and therefore the corresponding terms are twice as large as those for $a\neq b$.

With this further approximation, substituting Eq. (21) in Eq. (18) leads to

	$\displaystyle\left\langle\rho_{ab}^{2}\right\rangle{{}_{h|\rho_{g}}}=$	$\displaystyle\frac{1}{4}\iint\mathrm{d}f\,\mathrm{d}f^{\prime}\int\frac{\mathrm{d}\mathbf{n}}{{4\pi}}\,S_{h}(f,\mathbf{n})\chi_{ab}(\mathbf{n})\int\frac{\mathrm{d}\mathbf{n}^{\prime}}{{4\pi}}\,S_{h}(f^{\prime},\mathbf{n}^{\prime})\chi_{ab}(\mathbf{n}^{\prime})$
		$\displaystyle+\iint\mathrm{d}f\,\mathrm{d}f^{\prime}\mathrm{sinc}^{2}\left(\pi(f-f^{\prime})T\right)\int\frac{\mathrm{d}\mathbf{n}}{{4\pi}}\,S_{h}(f,\mathbf{n})\chi_{aa}(\mathbf{n})\int\frac{\mathrm{d}\mathbf{n}^{\prime}}{{4\pi}}\,S_{h}(f^{\prime},\mathbf{n}^{\prime})\chi_{bb}(\mathbf{n}^{\prime})$
		$\displaystyle+\frac{1}{4}\iint\mathrm{d}f\,\mathrm{d}f^{\prime}\mathrm{sinc}^{2}\left(\pi(f-f^{\prime})T\right)\int\frac{\mathrm{d}\mathbf{n}}{{4\pi}}\,S_{h}(f,\mathbf{n})\chi_{ab}(\mathbf{n})\int\frac{\mathrm{d}\mathbf{n}^{\prime}}{{4\pi}}\,S_{h}(f^{\prime},\mathbf{n}^{\prime})\chi_{ab}(\mathbf{n}^{\prime})\,,$		(23)

where

$$\chi_{ab}(\mathbf{n})=\sum_{A=+,\times}F_{a}^{A}(\mathbf{n})F_{b}^{A}(\mathbf{n})\,,$$

(24)

and we focused on non-vanishing pulsar pair separations only, i.e. $\rho_{ab}$ with $a\neq b$ (note that terms $\propto\chi_{aa},\chi_{bb}$ may nevertheless appear, for which we consider contributions from both the pulsar and Earth terms as discussed above).

The GW spectral density, $S_{h}(f,\mathbf{n})$, appearing in Eq. (23) can be expressed as the sum of an isotropic and anisotropic contribution, see Eq. (8). We therefore split the expression for $\left\langle\rho^{2}_{ab}\right\rangle{{}_{h|\rho_{g}}}$ into a “standard” and a “clustering” part,

$\displaystyle\left\langle\rho_{ab}^{2}\right\rangle{{}_{h|\rho_{g}}}=\left\langle\rho_{ab}^{2}\right\rangle^{\mathrm{st.}}_{h|\rho_{g}}+\left\langle\rho_{ab}^{2}\right\rangle^{\mathrm{clust.}}_{h|\rho_{g}}\,.$

(25)

The standard terms,

	$\displaystyle\left\langle\rho_{ab}^{2}\right\rangle^{\mathrm{st.}}_{h|\rho_{g}}=$	$\displaystyle\frac{1}{4}\iint\mathrm{d}f\,\mathrm{d}f^{\prime}\int\frac{\mathrm{d}\mathbf{n}}{{4\pi}}\,\bar{S}_{h}(f)\chi_{ab}(\mathbf{n})\int\frac{\mathrm{d}\mathbf{n}^{\prime}}{{4\pi}}\,\bar{S}_{h}(f^{\prime})\chi_{ab}(\mathbf{n}^{\prime})$
		$\displaystyle+\iint\mathrm{d}f\,\mathrm{d}f^{\prime}\mathrm{sinc}^{2}\left(\pi(f-f^{\prime})T\right)\int\frac{\mathrm{d}\mathbf{n}}{{4\pi}}\,\bar{S}_{h}(f)\chi_{aa}(\mathbf{n})\int\frac{\mathrm{d}\mathbf{n}^{\prime}}{{4\pi}}\,\bar{S}_{h}(f^{\prime})\chi_{bb}(\mathbf{n}^{\prime})$
		$\displaystyle+\frac{1}{4}\iint\mathrm{d}f\,\mathrm{d}f^{\prime}\mathrm{sinc}^{2}\left(\pi(f-f^{\prime})T\right)\int\frac{\mathrm{d}\mathbf{n}}{{4\pi}}\,\bar{S}_{h}(f)\chi_{ab}(\mathbf{n})\int\frac{\mathrm{d}\mathbf{n}^{\prime}}{{4\pi}}\,\bar{S}_{h}(f^{\prime})\chi_{ab}(\mathbf{n}^{\prime})\,,$		(26)

are exactly those specified in Eq. (C25) of Ref. Allen:2022dzg . The clustering terms, on the other hand, arise from the anisotropic part of $S_{h}(f,\mathbf{n})$ and are given by

	$\displaystyle{\left\langle\rho^{2}_{ab}\right\rangle^{\mathrm{clust.}}_{h|\rho_{g}}}=$	$\displaystyle\iint\mathrm{d}z\,\mathrm{d}z^{\prime}\left[G(z)G(z^{\prime})+\Gamma(z,z^{\prime})\right]\iint\frac{\mathrm{d}\mathbf{n}}{4\pi}\frac{\mathrm{d}\mathbf{n}^{\prime}}{4\pi}\delta_{g}(\mathbf{n},z)\delta_{g}(\mathbf{n}^{\prime},z^{\prime})\chi_{ab}(\mathbf{n})\chi_{ab}(\mathbf{n}^{\prime})$
		$\displaystyle+4\iint\mathrm{d}z\,\mathrm{d}z^{\prime}\,\Gamma(z,z^{\prime})\iint\frac{\mathrm{d}\mathbf{n}}{4\pi}\frac{\mathrm{d}\mathbf{n}^{\prime}}{4\pi}\delta_{g}(\mathbf{n},z)\delta_{g}(\mathbf{n}^{\prime},z^{\prime})\chi_{aa}(\mathbf{n})\chi_{bb}(\mathbf{n}^{\prime})$
		$\displaystyle+\mbox{terms linear in $\delta_{g}(\mathbf{n},z)$}\,,$		(27)

where

	$\displaystyle G(z)$	$\displaystyle\equiv\int_{0}^{\infty}\mathrm{d}f\,\bar{S}_{h}(f)b_{\rm GW}(f,z)\,,$
	$\displaystyle\Gamma(z,z^{\prime})$	$\displaystyle\equiv\frac{1}{4}\iint\mathrm{d}f\,\mathrm{d}f^{\prime}\mbox{sinc}^{2}\left(\pi(f-f^{\prime})T\right)\bar{S}_{h}(f)\bar{S}_{h}(f^{\prime})b_{\rm GW}(f,z)b_{\rm GW}({f^{\prime}},z^{\prime})\,,$		(28)

are measures of the strain weighted by the quantity $b_{\rm GW}(f,z)$ relating galaxy overdensity and GW intensity.⁹⁹9We note that, in the case where $b_{\rm GW}(f,z)=1$, $G$ and $\Gamma$ respectively correspond to the quantities $h^{2}$ and $\mathscr{h}^{4}$ defined in Ref. Allen:2022dzg . We note that, while terms linear in $\delta_{g}$ appear in Eq. (31), they vanish in the full ensemble average since $\left\langle\delta_{g}(\mathbf{n},z)\right\rangle_{\mathfrak{p}_{g}}=0$.

Taking the ensemble average of Eqs. (26) and (27) over all realizations of the galaxy density field leads to

$\displaystyle\left\langle\rho_{ab}^{2}\right\rangle=\left\langle\left\langle\rho_{ab}^{2}\right\rangle{{}_{h|\rho_{g}}}\right\rangle_{\mathfrak{p}_{g}}=\left\langle\rho_{ab}^{2}\right\rangle^{\mathrm{st.}}+\left\langle\rho_{ab}^{2}\right\rangle^{\mathrm{clust.}}\,,$

(29)

where, for the standard terms which are independent of the density field,

$$\left\langle\rho_{ab}^{2}\right\rangle^{\mathrm{st.}}=\left\langle\rho_{ab}^{2}\right\rangle^{\mathrm{st.}}_{h|\rho_{g}}$$

(30)

and, for the clustering terms,

	$\displaystyle\left\langle\rho^{2}_{ab}\right\rangle^{\mathrm{clust.}}=$	$\displaystyle\iint\mathrm{d}z\,\mathrm{d}z^{\prime}\left[G(z)G(z^{\prime})+\Gamma(z,z^{\prime})\right]\iint\frac{\mathrm{d}\mathbf{n}}{4\pi}\frac{\mathrm{d}\mathbf{n}^{\prime}}{4\pi}\xi_{g}(\mathbf{n}\cdot\mathbf{n}^{\prime},z,z^{\prime})\chi_{ab}(\mathbf{n})\chi_{ab}(\mathbf{n}^{\prime})$
		$\displaystyle+4\iint\mathrm{d}z\,\mathrm{d}z^{\prime}\Gamma(z,z^{\prime})\iint\frac{\mathrm{d}\mathbf{n}}{4\pi}\frac{\mathrm{d}\mathbf{n}^{\prime}}{4\pi}\xi_{g}(\mathbf{n}\cdot\mathbf{n}^{\prime},z,z^{\prime})\chi_{aa}(\mathbf{n})\chi_{bb}(\mathbf{n}^{\prime})\,.$		(31)

Here, we have used Eq. (5) in the form $\left\langle\delta_{g}(\mathbf{n},z)\delta_{g}(\mathbf{n}^{\prime},z^{\prime})\right\rangle=\xi_{g}(\mathbf{n}\cdot\mathbf{n}^{\prime},z,z^{\prime})$. Given Eq. (17), which states that $\left\langle\rho_{ab}\right\rangle$ is unaffected by clustering, we have

$$\left\langle\rho^{2}_{ab}\right\rangle-\left\langle\rho_{ab}\right\rangle^{2}=\left\langle\rho^{2}_{ab}\right\rangle^{\text{st.}}-\left(\left\langle\rho_{ab}\right\rangle^{\text{st.}}\right)^{2}+\left\langle\rho^{2}_{ab}\right\rangle^{\mathrm{clust.}}\,.$$

(32)

Eq. (31) presents the main result of our work: an analytical expression encoding the effect of galaxy clustering (described by the galaxy correlation function $\xi_{g}$) on the HD variance. In other words, this provides a measure of how much we expect the measured correlation of timing residuals to differ, due to galaxy clustering, from the idealized HD curve in a real measurement. We comment on the positivity of the contribution $\left\langle\rho^{2}_{ab}\right\rangle^{\mathrm{clust.}}$ to the total HD variance in Appendix A.

We close this section by noting that, instead of the galaxy correlation function, an alternative measure of galaxy clustering is provided by the angular power spectrum $C_{l}(z,z^{\prime})$,

$$\xi_{g}(\mathbf{n}\cdot\mathbf{n}^{\prime},z,z^{\prime})=\frac{1}{4\pi}\sum_{l}(2l+1)C_{l}(z,z^{\prime})\mathcal{P}_{l}(\mathbf{n}\cdot\mathbf{n}^{\prime})\,,$$

(33)

with $\mathcal{P}_{l}$ being the Legendre polynomials. This leads to

	$\displaystyle\left\langle\rho^{2}_{ab}\right\rangle^{\mathrm{clust.}}=$	$\displaystyle\iint\mathrm{d}z\,\mathrm{d}z^{\prime}\left[G(z)G(z^{\prime})+\Gamma(z,z^{\prime})\right]\sum_{l}\frac{2l+1}{4\pi}C_{l}(z,z^{\prime})\iint\frac{\mathrm{d}\mathbf{n}}{4\pi}\frac{\mathrm{d}\mathbf{n}^{\prime}}{4\pi}\,\mathcal{P}_{l}(\mathbf{n}\cdot\mathbf{n}^{\prime})\chi_{ab}(\mathbf{n})\chi_{ab}(\mathbf{n}^{\prime})$
		$\displaystyle+4\iint\mathrm{d}z\,\mathrm{d}z^{\prime}\,\Gamma(z,z^{\prime})\sum_{l}\frac{2l+1}{4\pi}C_{l}(z,z^{\prime})\iint\frac{\mathrm{d}\mathbf{n}}{4\pi}\frac{\mathrm{d}\mathbf{n}^{\prime}}{4\pi}\mathcal{P}_{l}(\mathbf{n}\cdot\mathbf{n}^{\prime})\chi_{aa}(\mathbf{n})\chi_{bb}(\mathbf{n}^{\prime})\,.$		(34)

As a further alternative, one can also choose to work with the galaxy power spectrum $P_{g}(k,z,z^{\prime})$,

$$\left\langle\delta_{g}(\mathbf{k},z)\delta_{g}^{\ast}(\mathbf{k}^{\prime},z^{\prime})\right\rangle=(2\pi)^{3}\delta_{D}(\mathbf{k}-\mathbf{k^{\prime}})P_{g}(k,z,z^{\prime})\,.$$

(35)

Here, $\delta_{g}(\mathbf{k},z)$ is the 3D Fourier transform of the galaxy overdensity $\delta_{g}(\mathbf{n},z)$. Using the galaxy power spectrum, Eq. (31) can be simplified and expressed as one instead of two integrals over redshift by using the Limber approximation limber1953analysis , a commonly used approximation in large-scale structure observables. In particular, following Ref. Bartelmann:1999yn , we obtain

	$\displaystyle\left\langle\rho^{2}_{ab}\right\rangle^{\mathrm{clust.}}=$	$\displaystyle\int\mathrm{d}z\left[{G^{2}(z)}+\Gamma(z,z)\right]{H(z)}\int\frac{k\,\mathrm{d}k}{2\pi}P_{g}(k,z)\iint\frac{\mathrm{d}\mathbf{n}}{4\pi}\frac{\mathrm{d}\mathbf{n}^{\prime}}{4\pi}\chi_{ab}(\mathbf{n})\chi_{ab}(\mathbf{n}^{\prime})J_{0}\big{(}kr_{z}\,{|\mathbf{n}-\mathbf{n}^{\prime}|}\big{)}$
		$\displaystyle+4\int\mathrm{d}z\,\Gamma(z,z){H(z)}\int\frac{k\,\mathrm{d}k}{2\pi}P_{g}(k,z)\iint\frac{\mathrm{d}\mathbf{n}}{{4\pi}}\frac{\mathrm{d}\mathbf{n}^{\prime}}{{4\pi}}\chi_{aa}(\mathbf{n})\chi_{bb}(\mathbf{n}^{\prime})J_{0}\big{(}kr_{z}\,{|\mathbf{n}-\mathbf{n}^{\prime}|}\big{)}\,,$		(36)

with $J_{0}$ denoting the 0th order Bessel function and $H(z)$ being the Hubble parameter.¹⁰¹⁰10The factor $H(z)$ appears due to the relation $\mathrm{d}z=H(z)\mathrm{d}r$, i.e. the change of integration variable from redshift $z$ to comoving distance $r$, which is needed to apply the Limber approximation (cf. Eq. (2.83) in Ref. Bartelmann:1999yn ).

This work introduces a rigorous theoretical framework to assess how the HD correlation, describing pulsar pair correlations caused by a stochastic GW background of astrophysical origin, is affected by anisotropies in the distribution of sources due to large-scale galaxy clustering. The presence of such anisotropies is naturally expected, since the GW sources (e.g. supermassive black holes at the center of galaxies) are expected to follow the distribution of matter in the Universe, which is itself anisotropically distributed.

We showed that to properly account for such anisotropies, it is necessary to consider not only ensemble averages over realizations of GW sources, but also ensemble averages over realizations of the galaxy density field. Modeling fluctuations around an isotropic universe by Gaussian random fields and performing ensemble averaging taking two layers of stochasticity into account, we find that the expectation value of the HD correlation is unchanged with respect to that of an isotropic universe. However, we find that anisotropies lead to a new contribution to the variance of the correlation, which has not been considered in previous work Allen:2022dzg ; Allen:2022ksj ; Romano:2023zhb ; Allen:2024rqk . Thereby, large-scale galaxy clustering impacts how much one expects observed timing residual correlations to deviate from the idealized HD curve.

This contribution is important to assess in view of future observations. While the present-day PTA data has error bars sufficiently large to fully overlap with the HD curve NANOGrav:2023gor , the precision of observations is expected to improve babak2024forecasting . Correctly modeling the variance of the HD correlation is therefore important in view of future observations. Indeed, a GW background of cosmological origin is only marginally affected by clustering (i.e. only via relativistic line-of-sight effects, which give very small contributions that we neglect in this work). As such, the clustering induced contribution to the variance evaluated in this work, which is only present in the case of an astrophysical origin, is a new interesting observable which may provide us with precious information on the nature of the stochastic GW background. Moreover, since various work on PTA observations has investigated the possibility to search for new physics (see e.g. Refs. liang2023test ; NANOGrav:2023hvm ; omiya2023hellings ), it is crucial to properly evaluate timing residual correlations accounting for all fundamental contributions to their variance, to enable an accurate comparison between future high-precision data and models.

In a companion paper Grimm:2024hgi , we apply our theoretical formalism to specific population models of supermassive black holes, considering as well the scenario of many pulsar pairs. In particular, in order to precisely quantify the size of the clustering induced variance, we need to determine how the distribution of sources is linked to the distribution of galaxies. This depends closely on the details of the source formation and evolution. We expect that different formation scenarios, leading to different source distributions, may result in different magnitudes of the HD variance, thus influencing how much we expect observations of the HD correlation to differ from the idealized HD curve. This work, in combination with its companion paper Grimm:2024hgi , provides a new framework to study the link between two important observational probes of the properties of the Universe: the large-scale clustering of galaxies and the stochastic GW background observed via PTAs.

We acknowledge discussions with Bruce Allen. 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”). C.B. is also supported by the Swiss National Science Foundation. 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”).