Forecasting pulsar timing array sensitivity to anisotropy in the stochastic gravitational wave background

Introduction

All of the long-baseline pulsar timing arrays (PTAs)—namely the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, McLaughlin [42]), the European Pulsar Timing Array (EPTA, Kramer & Champion [39]) and the Parkes Pulsar Timing Array (PPTA, Hobbs [31])—have now reported evidence for the presence of a spectrally-common process in their latest datasets ([9, 16, 25]). However, the evidence for Hellings & Downs (HD) cross correlations (Hellings & Downs [30]), which is considered the definitive signature of a stochastic gravitational wave background (GWB), was not significant in any of these datasets. A spectrally-common process was also reported by the International Pulsar Timing Array (IPTA, Hobbs et al. [32]) consortium in their second data release (DR2, Perera et al. [50]), which used older versions of datasets from the aforementioned PTAs. However this dataset also lacked definitive evidence for the HD signature (Antoniadis et al. [7]). Despite this, the detection of such a spectrally-common process is considered as the first step towards the eventual detection of a GWB (Romano et al. [55]). Based on simulations of the NANOGrav 12.5 yr dataset (Alam et al. [2]), and if the signal observed in the 12.5 yr dataset is an astrophysical GWB signal, NANOGrav is expected to have sufficient evidence to report the detection of HD-consistent correlations within the next few years (Pol et al. [52]).

Once this signal is confirmed to be a GWB signal in future PTA datasets, the onus will be on characterizing the source of the GWB. An important part of this analysis will be measuring the spatial distribution of power in the GWB. For example, if the source of the GWB is a population of inspiraling supermassive black hole binaries (SMBHBs, Rajagopal & Romani [53]; Jaffe & Backer [36]; Sesana et al. [57]; Burke-Spolaor et al. [14]), then their spatial distribution might track that of the local matter distribution. In particular, nearby galaxy clusters that host an overabundance of SMBHBs may show up as a hotspot of GW emission on the sky. For example, the Virgo cluster (Sandage [56]) has an angular diameter of $\sim10^\circ$, and could show up as a hotspot on the GW sky at a multipole of $l_{\mathrm{Virgo}} = 180^\circ/ \theta\approx18$. On the other hand, a single SMBHB that is louder than the GWB will show up as a point source anisotropy, with multiple such single sources producing a “popcorn”- or pixel-style spatial distribution of GWB power. However, if the GWB is produced by a cosmological source, such as cosmic strings (e.g., Olmez et al. [47]) or primordial GWs (e.g., Grishchuk [27]), then the GWB power distribution on the sky may not display the same anisotropies as that from a SMBHB-produced GWB. Thus, the anisotropy of the GWB, or the lack of it, allows us to make inferences about the source of the GWB.

Multiple techniques have been developed to probe the anisotropy of a GWB with PTAs (e.g., [18, 44, 61, 62]). While these methods differ in their choice of basis for modeling anisotropy, they all take the pulsar times-of-arrival (TOAs) as their initial data, and employ Bayesian techniques to constrain parameters that describe anisotropy-induced deviations away from the HD signature. Constraints on anisotropy were first presented by the EPTA as part of the analysis of their first data release ([63]), where they showed that the strain amplitude in $l > 0$ spherical-harmonic multipoles is less than 40% of the monopole value ($l = 0$, i.e., isotropy).

As PTA datasets grow longer in timespan, add more pulsars, and become denser with higher cadence observations, the analysis time for Bayesian methods based on TOAs are going to increase dramatically. Additionally, Bayesian model selection of anisotropy is predicated on having an appropriate hypothesis of the anisotropy, e.g., a power map built from galaxy catalogs, statistically populated with inspiraling SMBHBs. In other words, Bayesian model selection always requires two models to compare; a null and signal model. By contrast, frequentist techniques allow one to reject a null hypothesis (in this case, isotropy) if the data are in sufficient tension with it. Importantly, rejecting a null hypothesis at a probability given by some $p$-value is not the same as accepting a signal hypothesis with a certainty of $1-p$. However significant tension with the assumption of isotropy is an important indicator of beyond-HD signatures in the cross-correlation data.

To overcome these challenges, we develop a frequentist framework that employs the cross-correlations between pulsar timing residuals across a PTA as sufficient data with which to search for anisotropy. These cross correlations are measured as part of the standard GWB detection pipelines (e.g., [9, 16, 25]), using established optimal two-point correlation techniques ([5, 6, 19, 67]). Thus our framework can be easily incorporated into ongoing analysis campaigns. The

data volume in our framework is also significantly lower than analyses starting at the TOA-data level, enabling rapid estimation of anisotropy. Finally, as mentioned above, the frequentist framework allows us to infer (although not necessarily claim detection of) anisotropy via rejection of the null hypothesis of isotropy, thereby simplifying the process of constraining anisotropy.

Other frequentist techniques for searching for anisotropy with PTAs have been proposed. These include, for example, the Fisher matrix formalism developed in Ali-Haïmoud et al. [3] and Ali-Haïmoud et al. (2021) where they employ “principal maps”, which are the eigenmaps of the Fisher matrix, to search for anisotropy. Hotinli et al. [33], on the other hand, decompose the timing residual power spectrum onto bipolar spherical harmonics to search for anisotropy using the correlations between pulsars in a PTA. These frequentist techniques also focus on the detection of anisotropy through rejection of the null hypothesis of isotropy. The use of cross correlations between detector baselines is also commonly used when searching for anisotropic GWs in the LIGO (e.g., [64, 1]) and LISA bands (e.g., Banagiri et al. [12]).

We use our framework on simulations of idealized-PTA and near-future NANOGrav data, presenting, for the first time, projections on the sensitivity of NANOGrav and other PTAs to GWB anisotropy in the mid-to-late 2020s. As shown in Pol et al. [52], if the common-spectrum signal observed in the NANOGrav 12.5 yr dataset is a GWB, then NANOGrav can be expected to have sufficient statistical evidence to claim a detection of HD correlations within the next few years. By employing the same assumptions and simulation pipeline as Pol et al. [52] to generate our simulations, we forecast the level of anisotropy that would result in a statistically significant rejection of the null hypothesis of isotropy by the NANOGrav data.

The paper is organised as follows: in § 2, we describe the measurement of cross-correlations and their uncertainties, along with the maximum-likelihood framework and detection statistics that are used to constrain anisotropy, while § 3 connects the cross correlation uncertainty to the noise properties, cadence, and timing baseline of the PTA. § 4 presents scaling relations for anisotropy decision thresholds of an “ideal PTA” as a function of cross-correlation uncertainty and number of pulsars in a PTA, while § 5 presents them for a realistic PTA that is generated using the NANOGrav 12.5 yr dataset. Finally, we present a discussion of results and prospects for the future in § 6. In Appendix § 8 and Appendix B we present derivations for computing the cross correlations in real PTA datasets and scaling relations for the cross correlation uncertainty with respect to PTA specifications like timing baseline, white noise and observation cadence.

Methods

The optimal cross correlation statistic and overlap reduction function

A GWB can be uniquely identified through the cross correlations it induces between the times-of-arrivals (TOAs) of pulsars in a PTA [30, 65]. The timing cross correlations between pulsar pairs, $\rho_{ab}$, and their uncertainties, $\sigma_{ab}$, can be written as [19, 15, 58, 67]

where $\delta\mathbf{t}_a$ is a vector of timing residuals for pulsar-$a$, $\mathbf{P}_a = \langle\delta\mathbf{t}_a\delta\mathbf{t}_a^T\rangle$ is the measured autocovariance matrix of pulsar-$a$, and $\hat{\mathbf{S}}_{ab}$ is the template scaled-covariance matrix between pulsar-$a$ and pulsar-$b$. This scaled-covariance matrix is a template for the GWB spectral shape only, and is independent of the GWB’s amplitude and the cross-correlation signature that it induces. It is related to the full covariance matrix by $\mathbf{S}_{ab} = \langle\delta\mathbf{t}_a\delta\mathbf{t}_b^T\rangle= A_{\mathrm{gwb}}^2\chi_{ab}\hat{\mathbf{S}}_{ab}$, where $A_{\mathrm{gwb}}$ is the GWB amplitude corresponding to a given strain-spectrum template, and $\chi_{ab}$ is the GWB-induced cross-correlation value for this pair of pulsars, e.g., the Hellings & Downs factor in the case of an isotropic GWB.

The cross-correlation statistic accounts for fitting and marginalization over the timing ephemeris of each pulsar, as well as its intrinsic white and red noise characteristics, where a power-law spectrum is usually assumed for the intrinsic red noise. Implementations of the statistic also usually assume a power-law strain spectrum template for the GWB following $f^{-2/3}$ [51] filtering it across the pairwise-correlated timing residuals in order to extract an optimal measurement of the GWB amplitude [15, 67]. However, we note that this statistic is flexible enough to allow for different parametrized spectral templates. See Appendix § 8 for details of how this is implemented in real PTA searches for the GWB, including determining the noise weighting for each pulsar, setting the GWB’s spectral template, and how the parameters of the pulsar timing ephemeris are marginalized over. In this paper we simulate and analyze data at the level of $\{\rho_{ab},\sigma_{ab}\}$, then connect our results later to the underlying geometry and noise characteristics of the PTA (see Appendix B).

A PTA with $N_{\mathrm{psr}}$ pulsars has $N_{\mathrm{cc}} = N_{\mathrm{psr}}(N_{\mathrm{psr}} - 1)/2$ distinct cross-correlation values. The angular dependence of these empirically measured cross-correlations can be modelled by the detector overlap reduction function (ORF)¹ [22, 44, 61, 24, 62] such that, for an unpolarized GWB,

where $P(\hat{\Omega})$ is the angular power of the GWB in direction $\hat{\Omega}$, normalized such that $\int_{S^2}d^2\hat{\Omega}\,P(\hat{\Omega}) = 1$, and $\mathcal{F}^{A}(\hat{p},\hat{\Omega})$ is the antenna response pattern of a pulsar in unit-vector direction $\hat{p}_a$ to each GW polarization $A \in\{+,\times\}$:

where $e^{A}_{ij}(\hat{\Omega})$ are polarization basis tensors, and $(i,j)$ are spatial indices.

We can recast the sky integral in Equation (2) as a sum over equal-area pixels [24, 62]. Assuming an unpolarized GWB, and ignoring random pulsar term contributions to the cross correlations, this can be written as

where $k$ denote pixel indices. Or, in a general matrix form

where $\mathbf{\Gamma}$ is an $N_{\mathrm{cc}}$ vector of ORF values for all distinct pulsar pairs, $\mathbf{P}$ is an $N_{\mathrm{pix}}$ vector of GWB power values at different pixel locations, and $\mathbf{R}$ is a $(N_{\mathrm{cc}} \times N_{\mathrm{pix}})$ overlap response matrix given by

where the normalization is chosen so that the ORF matches the Hellings & Downs values in the case of an isotropic GWB with $P_k = 1\ \forall k$.

Thus, in our notation the expected value of $\rho_{ab}$ is such that $\langle\rho_{ab}\rangle= A^2\Gamma_{ab}$. However, in the remainder of this paper we will deal with amplitude-scaled cross-correlation values, $\rho_{ab}/A^2$, where we assume that in a real search an initial fit for $A^2$ will be performed on $\{\rho_{ab}\}$ under the assumption of isotropy. These amplitude-scaled cross-correlation values can then be directly compared with the ORF model to probe anisotropy. We suppress the explicit amplitude scaling in the remainder of our notation, such that $\rho_{ab}$ henceforth implies amplitude-scaled cross-correlation values.

Assuming a stationary Gaussian distribution for the cross correlation uncertainty, the likelihood function for the cross correlations can be written as

where $\boldsymbol{\Sigma}$ is the diagonal covariance matrix of cross-correlation uncertainties. We now discuss several different bases on which to decompose the angular power vector, $\mathbf{P}$.

Pixel basis

The GWB power can be parametrized using HEALPix sky pixelization Górski et al. [26], where each equal-area pixel is independent of the pixels surrounding it

The number of pixels is set by $N_{\mathrm{pix}} = 12N_{\mathrm{side}}^{2}$ and $N_{\mathrm{side}}$ defines the tessellation of the healpix sky Górski et al. [26]. The rule of thumb for PTAs is to have $N_{\mathrm{pix}} \le N_{\mathrm{cc}}$ [18, 54]. This basis is well suited for detection of pixel-scale anisotropy, which can arise from individual sources of GWs, and where an isotropic GWB would be represented by equal power (within uncertainties) in each pixel on the sky.

Spherical harmonic basis

Alternatively, the GWB power can be decomposed onto the spherical harmonic basis (e.g. Thrane et al. [64]), where the lowest order multipole ($l = 0$) defines an isotropic background, while higher multipoles add anisotropy. The GWB power in this basis can be written as

where $Y_{lm}$ are the real valued spherical harmonics and $c_{lm}$ are the spherical harmonic coefficients of $P(\hat{\Omega})$. In this basis, the ORF anisotropy coefficient for the $l,m$ components between pulsars $a,b$ can be written as Taylor et al. [62],

where $k$ represents the pixel index corresponding to $\hat{\Omega}$ and the constant $\kappa$ accounts for the pixel area in the healpix sky tessellation.

This basis representation, contrary to the pixel basis, is better suited for modeling large-scale anisotropies in the GWB. Based on diffraction-limit arguments, the highest order mode, $l_{\max}$, that can be used for modeling the anisotropy depends on the number of pulsars in the PTA, $l_{\max} \sim\sqrt{N_{\mathrm{psr}}}$ [13, 54]. However, Floden et al. [23] have shown that while the diffraction limit is attuned to maximizing the significance of the detection of anisotropy, values of $l > l_{\max}$ can be included in spherical harmonic decompositions to improve the localization of any anisotropy after its detection.

The results can be expressed in terms of $C_l$, which is the squared angular power in each mode $l$

Physically, $C_l$ represents the amplitude of statistical fluctuations in the angular power of the GWB at scales corresponding to $\theta= 180^\circ/l$. An isotropic background in this basis will contain power only in the $l = 0$ multipole, thus filling the entire sky, while an anisotropic background will have power in the higher $l$ multipoles. On the other hand, the variance of the angular power distribution can be written as Jenkins [37]

The quantity $l(l+1)C_l/2\pi$ thus represents the variance per logarithmic multipole bin, and is what we use to present our results in this work. As pointed out in Jenkins [37], reporting $C_l$ is analogous to reporting the GWB strain power spectral density, with $C_l = \text{constant}$ representing a white angular power spectrum, while reporting $l(l+1)C_l/2\pi$ is analogous to reporting the GWB energy density spectrum $\Omega_{\mathrm{gw}}(f)$, with $l(l+1)C_l/2\pi= \text{constant}$ representing a scale-invariant angular power spectrum.

For these bases, since the problem is linear in the regression coefficients, the maximum likelihood solution can be derived analytically [64, 54, 35]

where $\mathbf{M} \equiv\mathbf{R}^{\mathsf{T}}\boldsymbol{\Sigma}^{-1}\mathbf{R}$ is the Fisher information matrix, with the uncertainties on the $c_{lm}$ coefficients given by the diagonal elements of $\mathbf{M}^{-1}$, and $\mathbf{X} \equiv\mathbf{R}^{\mathsf{T}}\boldsymbol{\Sigma}^{-1}\boldsymbol{\rho}$ is the “dirty map”, an inverse-noise weighted representation of the total power on the sky as “seen” through the response of the pulsars in the array.

Square-root spherical harmonic basis

A drawback of both the pixel and spherical harmonic bases is that they allow the GWB power to assume negative values, which is an unphysical realization of the GWB. While these tendencies can be curbed through the use of regularization techniques or rejection priors [61], this results in the addition of a hyperparameter to the analysis that requires further optimization or non-analytic priors [61, 62]. A more elegant solution is to use a basis that intrinsically conditions the GWB power to be positive over the whole sky.

Such a basis can be generated by modeling the square-root of the GWB power, $P(\hat{\Omega})^{1/2}$, rather than modeling the power itself. This technique was introduced in a Bayesian context in Payne et al. [49] for LIGO, Banagiri et al. [12] for LISA, and in Taylor et al. [62] for PTAs. Decomposing the square-root power onto spherical harmonics, the GWB power can be written as

where $Y_{LM}$ are real valued spherical harmonics and $b_{LM}$ are the search coefficients. Banagiri et al. [12] showed that the search coefficients in this basis can be related to the spherical harmonic coefficients via

where $\beta_{lm}^{LM,L'M'}$ is defined as

with $C_{LM,L'M'}^{lm}$ being the Clebsch-Gordon coefficients. Taylor et al. [62] showed that even though a full reconstruction of $c_{lm}$ requires an infinite sum over the $b_{LM}$ coefficients, restricting the maximum mode to $L \leq L_{\max}=l_{\max}$ is sufficient to produce an accurate reconstruction of the GWB power.

Since the problem in this basis is non-linear in the regression coefficients, the likelihood in Equation (7) cannot be maximized analytically. The maximum likelihood solution thus has to be calculated through numerical optimization techniques. In this work, we use the lmfit (Newville et al. [45]) Python package, where we use the Levenberg-Marquardt (LM) optimization algorithm [40, 41] to determine the maximum likelihood solution. The goodness-of-fit is assessed through the $\chi^2_{\mathrm{dof}}$ that is reported by LMFIT.

Figure 1 shows an example of recovering an anisotropic background using this basis. To produce the simulated cross-correlation data, we inject a GW power map corresponding to the synthesized population of SMBHBs from Taylor et al. [62] into an “ideal PTA” consisting of 100 pulsars, with a constant cross-correlation measurement uncertainty of 0.01. Note that this is an uncertainty on cross correlations that can assume any value between $-0.2$ and 0.5, and thus represents extremely accurate measurement of the cross correlations between pulsars in the PTA (see § 4.1 for our definition of an “ideal PTA”). We see that this basis is capable of reproducing the injected angular power spectrum, which translates into an accurate identification of the anisotropy in the GWB.

Detection statistics

In addition to estimating the values of the anisotropy search coefficients (i.e., pixels or spherical harmonic coefficients), we also need to quantify the evidence for the presence of anisotropy in the cross correlation data. As described earlier, when searching for anisotropy, we seek to reject the null hypothesis of isotropy. In this section, we present two frequentist detection statistics that can quantify the evidence for anisotropy by measuring the (in)compatibility of the data with the null hypothesis.

Signal-to-noise ratio

Confidence in the detection of anisotropy can be quantified through the signal-to-noise (S/N) ratio, which is defined as the ratio of the maxima of the likelihood functions between any two models

where $\Lambda_{\mathrm{ML}} = p(\boldsymbol{\rho}|\mathbf{P}_{\mathrm{ML},1})/p(\boldsymbol{\rho}|\mathbf{P}_{\mathrm{ML},2})$ is the ratio of the maxima of the likelihood functions for the two models that are under consideration.

When searching for anisotropy, we can define three S/N statistics which together provide a complete description of the evidence for the GWB signal present in the cross correlation data:

1. “Total S/N”: This is defined as the ratio of the maximum likelihood values of an anisotropic model (i.e. $p(\boldsymbol{\rho}|\mathbf{P}_{\mathrm{ML}}(l_{\max}>0))$) to a model with only spatially-uncorrelated noise (i.e., $p(\boldsymbol{\rho}|\mathbf{P}_{\mathrm{ML}}=0)$). The total S/N quantifies the evidence for the presence of any signal in the cross correlations.

2. “Isotropic S/N”: This is defined as the ratio of the maximum likelihood values of an isotropic model (i.e., $p(\rho|\mathbf{P}_{\mathrm{ML}}(l_{\max}=0))$) to a model with noise only (i.e., $p(\rho|\mathbf{P}_{\mathrm{ML}}=0)$). Note that the isotropic S/N is equivalent to the optimal S/N statistic defined in Chamberlin et al. [15]. The isotropic S/N ratio quantifies how well the cross correlations are described by a purely isotropic model.

3. “Anisotropic S/N”: This is defined as the ratio of the maximum likelihood values of a model with anisotropy (i.e., $p(\rho|\mathbf{P}_{\mathrm{ML}}(l_{\max}>0))$) to an isotropic model (i.e., $p(\rho|\mathbf{P}_{\mathrm{ML}}(l_{\max}=0))$). The anisotropic S/N ratio quantifies the evidence in favor of inclusion of modes $l>0$.

Decision threshold

Another method for determining the significance of possible anisotropy is to assess the certainty with which the null hypothesis of isotropy can be rejected. If we can quantify the distribution of the angular power, $C_l$, under the null hypothesis, then we can also quantify how (in)consistent the measured angular power is with isotropy through a test statistic like the $p$-value.

To calculate the null distribution of $C_l$ under the null hypothesis, we generate many realizations of cross-correlation data, where we assume that the measurements are Gaussian distributed around the Hellings & Downs curve, with the spread of the distribution given by the uncertainty on the cross correlation values. We define the decision threshold, $C_l^{\mathrm{th}}$, as the value of $C_l$ corresponding to a $p$-value of $3 \times10^{-3}$, where a measurement of angular power greater than this threshold would indicate a tension with the null hypothesis at the $\sim3\sigma$ level.

Cross-correlation uncertainties from PTA design

The uncertainty on the cross-correlation measurements introduced in Equation (1) depends on the pulsars that are used to construct the PTA ([6, 15, 58]). As shown in Siemens et al. [58], the trace in Equation (1) can be written as

where $T$ is the duration of time that a pulsar is observed, i.e., the timing baseline, $A$ is the amplitude of the GWB, $P_g(f)$ is the power spectrum of the timing residuals induced by the GWB, $P_a(f)$ and $P_b(f)$ are the intrinsic power spectra of pulsars $a$ and $b$, $f_l=1/T$ and $f_h$ are the low and high frequency cutoffs used in the GWB analysis.

Using the same analysis as presented in Siemens et al. [58] and Chamberlin et al. [15], we can show that (see Appendix B) in the weak-signal regime, the cross-correlation uncertainty between a given pair of pulsars scales as

while in the strong-signal regime

where $w$ is the white noise RMS, $c = 1/\Delta t$ is the observing cadence, $\gamma= 3 - 2\alpha$ Arzoumanian et al. [9] is the slope of the timing-residual power spectrum induced by the GWB, and $A$ is the amplitude of the GWB whose characteristic strain spectrum is given by

with $\alpha= -2/3$ for a SMBHB GWB.

These scaling relations imply that the uncertainty on the cross correlations can be reduced by observing pulsars for longer duration (timing baseline), or with a higher cadence. The effect of the timing baseline and cadence is strongest in the weak-signal regime while it becomes weaker as we move into the strong-signal regime. Similarly, in the weak-signal regime, the cross-correlation uncertainty can be reduced by decreasing the intrinsic white noise by, for example, increasing the receiver bandwidth or increasing the integration time for the pulsars. The white noise does not affect the cross-correlation uncertainty in the strong-signal regime, where the GWB signal dominates over the intrinsic noise in the pulsars at all frequencies.

Simulations with an ideal PTA

Defining an ideal PTA

In the framework described above, the rejection of isotropy in the GWB depends on three variables: the uncertainty on the cross correlation values, the number of pulsars in the PTA (which defines the number of cross correlation values that are measured), and the distribution of the pulsars on the sky. The first two variables primarily dictate the strength of the rejection of isotropy, while the third variable is important for the characterization of the anisotropy. In this section, we examine how the cross correlation uncertainty and the number of pulsars in the PTA affect an ideal PTA, which we define as a PTA that has pulsars distributed uniformly on the sky, and all pulsars in the PTA having identical noise properties. The latter constraint implies that all measured cross-correlations have the same uncertainty. This is different from current, real PTAs, where each pulsar is unique and thus the uncertainties on the cross correlations between each pulsar pair are different. We examine realistic PTAs in §5.

Scaling relations

For a given level of anisotropy, its detection significance will depend on the number of pulsars in the array, as well as on the accuracy with which the cross correlations between different pulsar pairs can be measured. Romano & Cornish [54] showed that the total S/N for linear anisotropy models can be written as

where $\hat{\mathbf{P}}$ is the maximum-likelihood estimate of the GWB power. Since $\mathbf{\Sigma}$ is a diagonal matrix of the squared cross correlation uncertainties, the total S/N $\propto\sigma^{-1}$. Similarly, the total S/N is proportional to the square root of the number of data points available for inference (assuming all cross correlation uncertainties are the same for all pairs), i.e.

where for sufficiently large $N_{\mathrm{psr}}$, the total S/N will scale linearly with the number of pulsars in the PTA.

We show that both of these scaling relations for the total S/N are also satisfied when using the non-linear maximum likelihood approach described in § 2.1.3. Figure 2 and the left panel in Figure 3 show that the total S/N scales inversely with the uncertainty on the cross correlations, while the right hand panel of Figure 3 show that the total S/N scales proportionally to the number of pulsars in the PTA. Since the injected signal here is an isotropic GWB, Figure 2 shows that the total and isotropic S/N values are identically large, and decrease as the uncertainty on the cross correlations increases. The anisotropic S/N shows little evolution across uncertainties, though the better fit provided by the additional degrees of freedom in an anisotropic model prevent it from being consistent with zero for small cross-correlation uncertainties. For large uncertainties, Figure 2 shows that we lose the ability to detect and distinguish between isotropic or anisotropic signals in the data.

We can similarly compute scaling relations for the decision threshold as a function of the cross correlation uncertainty and number of pulsars in the PTA. Since this is an empirically constructed detection statistic, we do not have analytical expressions for its scaling relations, though we can derive the scaling expressions computationally, as shown in Figure 4. As expected, as we increase the number of pulsars in the PTA, the decision threshold decreases across all multipoles. Similarly, as the uncertainty on the cross correlation measurements decreases, so do the multipole-dependent decision thresholds, corresponding to an improved sensitivity to deviations away from isotropy.

Sensitivity figure of merit

Rather than treating the number of pulsars and the cross-correlation uncertainty as separate variables, we can consider a combination that is inspired by the weighted arithmetic mean of cross-correlation measurements involved in, e.g., S/N calculations. In such calculations, operations like $\sum_{ab}(\cdots)/\sigma_{ab}^{2}$ are proportional to $N_{\mathrm{cc}}/\sigma^{2} \propto(N_{\mathrm{psr}}/\sigma)^{2}$ in the limit of equal cross-correlation measurement uncertainties, or where we can characterise the distribution of uncertainties by its mean or median over pairs. Therefore we define a sensitivity figure of merit (FOM), $N_{\mathrm{psr}}/\sigma$, and quantify the dependence of our detection statistics with respect to this. This also allows us to quantify the trade-off between the number of pulsars in the PTA and the cross-correlation uncertainties, which, in turn, are related to the noise characteristics of the pulsars.

The relation of the total S/N and decision threshold to the sensitivity FOM is shown in Figure 5. We confirm that the total S/N is proportional to the sensitivity FOM, $N_{\mathrm{psr}}/\sigma$ in logarithmic space. This implies that, as expected, fewer pulsars or larger uncertainties on the cross correlation measurements result in a reduced total S/N, while larger numbers of pulsars or smaller uncertainties result in an increase in the total S/N value. Similarly, Figure 5 shows the dependence of the decision threshold for each anisotropy multipole on $N_{\mathrm{psr}}/\sigma$. Similar to the total S/N, a PTA with fewer pulsars or larger uncertainties on the cross correlation measurements will be able to reject the null hypothesis with lower significance than a PTA with more pulsars and/or smaller uncertainties on the cross correlations.

Simulations with a realistic PTA

The ideal PTA described in § 4 was useful in discerning scaling relationships that map between the PTA design and detection statistics. However, unlike the ideal PTA in § 4, real PTAs do not (yet) consist of pulsars distributed uniformly across the sky, nor are all the pulsars in the array identical. The latter fact implies that the cross-correlation uncertainties in a real PTA will be described by a distribution, rather than a constant value as assumed in § 4.

To simulate a realistic PTA, we use the methods developed in Pol et al. [52]. We base our simulations on the NANOGrav 12.5 year dataset (Alam et al. [2]), and extend the dataset to a 20-year timing baseline to forecast the sensitivity of NANOGrav to anisotropies in the GWB. The TOA timestamps of the initial 12.5 year portion are the same as those in the real NANOGrav dataset, while the radiometer uncertainties and pulse-phase jitter noise that are injected are obtained from the maximum-likelihood pulsar noise analysis performed as part of the NANOGrav 12.5 year analysis (Arzoumanian et al. [9]). The injection values for the intrinsic per-pulsar red noise were taken from a global PTA analysis that also modeled a common-spectrum process. This is done to isolate the intrinsic red noise in each pulsar’s dataset so that it is not contaminated by the common process reported in Arzoumanian et al. [9].

Once the 45 simulated pulsars from the NANOGrav 12.5 year dataset are generated using the above recipe, the dataset is then extended into the future by generating distributions for the cadence and measurement uncertainties using the last year’s worth of data for each pulsar. We then draw TOAs using these distributions until the dataset has a maximum baseline of 20 years. Finally, we inject 100 statistically random realizations of an isotropic gravitational wave background with an amplitude of $A_{\mathrm{GWB}} = 2 \times10^{-15}$ and spectral index $\alpha= -2/3$, consistent with the common process observed in Arzoumanian et al. [9].

We then pass all 100 realizations of the dataset through the standard NANOGrav detection pipeline (Arzoumanian et al. [8, 9]; ); Pol et al. [52]) to calculate the cross correlations and their uncertainties between all pairs in the 45 pulsar dataset (see also Appendix § 8). The evolution of these cross-correlation uncertainties across the 100 realizations as a function of the

iliar to the total S/N, a PTA with fewer pulsars or larger uncertainties on the cross correlation measurements will be able to reject the null hypothesis with lower significance than a PTA with more pulsars and/or smaller uncertainties on the cross correlations. timing baseline is shown in Figure 6. As we can see, the median cross-correlation uncertainty reduces from $\sim5$ at 13 years (similar to the total baseline of the 12.5 year dataset) to $\sim1$ at 20 years, implying a scaling relation $\sigma\propto T^{-7/2}$. This is shallower than the spectral index predicted for the weak-signal regime in § 4.2, which implies that the NANOGrav PTA is in the intermediate-signal regime, which is corroborated by the fact that the lowest frequencies of the PTA are dominated by a common-spectrum process (interpreted as a GWB) as shown in Arzoumanian et al. [9]. Combining the median uncertainty with the 45 pulsars in the PTA, we obtain values for our sensitivity FOM $N_{\mathrm{psr}}/\sigma\approx9$ at 13 years, and $N_{\mathrm{psr}}/\sigma\approx45$ at 20 years.

We pass the cross correlations measured from these 100 realizations through the statistical framework described in § 2 to search for the presence of anisotropy under realistic PTA and data-quality conditions. The evolution of the three S/N statistics as a function of time are shown in Figure 7. Since the injected GWB is isotropic, the total and isotropic S/N increase with the timing baseline. This is consistent with the reduction in the uncertainties on the cross correlations allowing for a stronger detection of the isotropic background. By contrast, the anisotropic S/N does not increase with time, and has support at S/N = 0 for all baselines. Note that the total S/N seen in these realistic simulations is consistent with the prediction made in Figure 5 and Pol et al. [52]. Similarly, Figure 8 shows the evolution of the decision threshold as a function of spherical harmonic multipole $l$ and the timing baseline. We find that the anisotropy decision threshold is such that, in terms of the $C_l$ values, GWB anisotropies at levels $C_{l=1}/C_{l=0} \gtrsim0.3$ (i.e. greater than 30% of the power in the monopole) would be inconsistent with the null hypothesis of isotropy at the $p = 3 \times10^{-3}$ level for the 20 yr baseline.

For comparison, in Figure 8 we also plot the Bayesian 95% upper limits on GWB anisotropy using six pulsars from the EPTA’s first data release (Taylor et al. [63]) for a model extending to $l_{\max} = 4$. This dataset had a maximum baseline of 17.7 years, which is toward the upper end of the baselines that we simulate for the NANOGrav data. However, the number of pulsars (6) in the EPTA analysis is significantly lower than the number of pulsars in our simulations (45). The longer EPTA timing baseline allows the Bayesian EPTA upper limits and NANOGrav anisotropy decision thresholds to be comparable at low multipoles until NANOGrav’s timing baseline exceeds that of the EPTA. However, the larger number of pulsars in NANOGrav not only gives it higher spatial resolution (and thus access to higher multipoles), but also improves the sensitivity of NANOGrav at higher multipoles relative to the EPTA 2015 limit. As shown in Floden et al. [23], access to these higher multipoles can aid in the localization of GWB anisotropies caused by individual GW sources, or due to finiteness in the source population constituting the GWB. This highlights the importance of including more pulsars in a

PTA for spatial resolution, even with the shorter timing baseline of such new additions.

Discussion and Conclusion

We have explored the detection and characterization of stochastic GWB anisotropy through pulsar cross correlations in a PTA. Using a frequentist maximum likelihood approach, we can search for anisotropy by modeling GWB power in individual sky pixels or through a weighted sum of spherical harmonics. Anisotropy would then manifest in measured cross-correlations between pulsar timing residuals through a power-weighted

overlap of pulsar GW antenna response functions. As a refinement on previous approaches, we prevent the GWB power from assuming unphysical negative values by adopting a model that naturally restricts this; we have referred to this as the square-root spherical harmonic basis throughout our analysis. We have also defined two detection metrics: (a) the signal-to-noise ratio, S/N, defined as the ratio between the maximum likelihood values between a signal and noise model, and (b) the anisotropy decision threshold, $C_l^{\mathrm{th}}$, defined as the level at which the measured angular power is inconsistent with isotropy at the $p = 3 \times10^{-3}$ ($\sim3\sigma$) level. The S/N comes in three flavors: (i) the total S/N, which measures the strength of an anisotropic GWB model against noise alone; (ii) the isotropic S/N which measures the strength of an isotropic GWB model against noise alone (this directly corresponds to the usual optimal statistic S/N used in PTA data analysis); and (iii) the anisotropic S/N, which measures the statistical preference for anisotropy against isotropy.

We examined the evolution of these detection statistics as a function of the uncertainty on the measured cross correlations, as well as on the number of pulsars in the PTA. We showed that the cross-correlation uncertainty and the number of pulsars in a PTA can be combined into a single figure of merit for the PTA sensitivity, $N_{\mathrm{psrs}}/\sigma$, which succinctly maps the PTA configuration and noise specifications to the detectability of anisotropy. Our scaling relations show that increasing the number of pulsars in an array, while reducing the uncertainty on the cross correlation measurements, leads to higher total S/N and lower anisotropy decision thresholds. As shown in § 3, the cross-correlation uncertainty scales inversely with both the timing baseline and cadence of observation for each pulsar, with the power-law index dependent on the signal regime occupied by the PTA. The pulsar timing baseline is set to increase as PTAs continue operation into the future, as well as through IPTA data combinations (e.g., Perera et al. [50]). Improving observing cadence is more challenging due to constraints on the available telescope time for each PTA, though once again IPTA data combinations can help alleviate this problem. In addition to IPTA data combinations, the CHIME radio telescope (CHIME/Pulsar Collaboration et al. [17]) will offer $\sim$daily cadence to the NANOGrav PTA (McLaughlin [42]), which is an order of magnitude improvement over the current $\sim$monthly cadence employed by NANOGrav.

Finally, we examined the evolution of the S/N and anisotropy decision thresholds as a function of timing baseline using realistic NANOGrav data. Since we injected an isotropic GWB signal in these data, we found that the anisotropic S/N remains consistent with zero at all times and across all signal realizations. By contrast, the total and isotropic S/N increase with time, as expected (Siemens et al. [58]). We find that any anisotropic GWB power distribution with $C_{l=1} \gtrsim0.3C_{l=0}$ would be in tension with an isotropic model at the $p = 3 \times10^{-3}$ ($\sim3\sigma$) level. We note that these simulations held the number of pulsars in the array fixed to the 45 that were included in the NANOGrav 12.5 year dataset. However, this number will increase in future NANOGrav datasets. Based on results in § 4.3, this will lead to larger total S/N values and lower anisotropy decision thresholds, implying improved sensitivity to any anisotropy that might be present in the real GWB. IPTA data combinations will further increase the timing baseline and number of pulsars in the array allowing further improvements on the ability to detect anisotropy. Furthermore, new instruments such as ultra-wideband receivers, and new telescopes such as MeerKAT (Bailes et al. [11]), SKA (Bailes et al. [10]), and DSA-2000 (Hallinan et al. [28]) will also aid in the detection and characterization of anisotropy with PTAs.

The techniques and scaling relations that we have developed in this work are PTA-agnostic and can be projected onto any PTA specification, allowing for immediate usage by the broader PTA community. Yet we have made assumptions in our framework that can be generalized in the future. For example, while our techniques operate on PTA data at the level of cross correlations rather than TOAs, in order to get to that stage we have implicitly assumed that the GWB characteristic strain spectrum is well-described by a power-law model. This follows the same approach as Arzoumanian et al. [9], where an average power-law spectrum $h_c \propto f^{-2/3}$ is assumed for the GWB. For a GWB produced by a population of inspiraling SMBHBs, this power-law representation is an approximation to the true spectrum [51], where there are different SMBHBs contributing to the GWB at different frequencies [38]. Thus, a more appropriate way to characterize anisotropy in the SMBHB GW background would be to measure the cross correlations of pulsar TOAs as a function of GW frequency, rather than our current approach of computing cross correlations that are filtered against a power-law GWB spectral template (see Appendix § 8). We would then have a more general data structure that includes pulsar cross-correlations and uncertainties for each GW frequency, for which the methods developed here can be applied at each of those frequencies independently. We also note that the methods developed here can be modified to search for multiple backgrounds [59], where an astrophysical background (e.g. from SMBHBs, Kelley et al. 2017; Sesana et al. 2004; Burke-Spolaor et al. 2019) would be expected to be anisotropic, while a cosmological background (e.g. from cosmic strings, Ölmez et al. 2012) may be isotropic. We plan to explore these improvements and generalizations in future analyses in a bid to extract as much spatial and angular information as possible from the exciting new PTA datasets now under development. As mentioned, these techniques will not only aid in the detection of GWB anisotropy, but also in its characterization for the purposes of isolating regions of excess power that may be indicative of individually-resolvable GW sources, and as leverage for the separation of potentially multiple stochastic GW signals of astrophysical and cosmological origin.

NP was supported by the Vanderbilt Initiative in Data Intensive Astrophysics (VIDA) Fellowship. SRT acknowledges support from NSF AST-200793, PHY-2020265, PHY-2146016, and a Vanderbilt University College of Arts & Science Dean’s Faculty Fellowship. JDR was partially supported by start-up funds provided by Texas Tech University.

This work has been carried out by the NANOGrav collaboration, which is part of the International Pulsar Timing Array. The NANOGrav project receives support from National Science Foundation (NSF) Physics Frontiers Center award number #1430284 and #2020265.

This work was conducted using the resources of the Advanced Computing Center for Research and Education at Vanderbilt University, Nashville, TN.

SOFTWARE

We use lmfit (Newville et al. [45]) for the non-linear least-squares minimization in the square-root spherical harmonic basis. We use libstempo (Vallisneri et al. [66]) to generate our realistic PTA datasets and to inject the pulsar noise parameters and GWB signals in these datasets. We use the software packages Enterprise (Ellis et al. [20]) and enterprise extensions (Taylor et al. [60]) for model construction, along with PTMCMCSampler (Ellis & van Haasteren [21]) as the Markov Chain Monte Carlo sampler for our realistic PTA Bayesian analyses. We also extensively used Matplotlib (Hunter [34]), NumPy (Harris et al. [29]), Python (Oliphant [46]; ); Millman & Aivazis [43]), and SciPy (Virtanen et al. [68]).

Computing cross-correlations in realistic pulsar-timing datasets

Reiterating Equation (1), the cross-correlation value (and its uncertainty) measured between pulsar-$a$ and pulsar-$b$ are

where $\delta t_a$ is a vector of timing residuals for pulsar-$a$, $P_a = \langle\delta t_a\delta t_a^T\rangle$ is the measured autocovariance matrix of pulsar-$a$, and $\hat{S}_{ab}$ is the template scaled-covariance matrix between pulsar-$a$ and pulsar-$b$. This scaled-covariance matrix is a template for the spectral shape only, and is amplitude and ORF independent. It is related to the full covariance matrix by $S_{ab} = \langle\delta t_a\delta t_b^T\rangle= A^2\chi_{ab}\hat{S}_{ab}$.

Covariance matrices in the PTA likelihood are modeled using rank-reduced approximations, e.g., long-timescale stochastic processes are modeled on a truncated Fourier basis with a number of frequencies that is much smaller than the number of TOAs. The induced timing delays of such stochastic processes are modeled as $\mathbf{r} = \mathbf{T}\mathbf{b}$, where $\mathbf{b}$ is a vector of amplitude coefficients for the basis design matrix $\mathbf{T} = [\mathbf{M}\ \mathbf{F}]$, which itself is a concatenation of the stabilized timing-model design matrix $\mathbf{M}$ and Fourier design matrix $\mathbf{F}$. The former is a matrix of partial derivatives of TOAs with respect to timing model parameters, which is then column-normed or SVD-stabilized to reduce the dynamic range of the entries. The latter is a matrix of sine and cosine basis functions evaluated at all TOAs for each sampling frequency of the time series. The autocovariance matrix of pulsars is then modeled as $P_a = N_a + T_aB_{aa}T_a^T$, where $N_a$ is the white-noise covariance matrix of pulsar-$a$, and

is the covariance matrix of the $\mathbf{b}$ coefficients, $B = \langle\mathbf{b}\mathbf{b}^T\rangle$, with infinite variance in the timing model diagonal entries to approximate an unbounded uniform prior, and $\phi$ as the variance of the Fourier coefficients that is related to the power spectral density of the stochastic process. The inter-pulsar covariance matrix is treated similarly, except that the timing model entries can be completely ignored since they will be uncorrelated between different pulsars. Hence, $\hat{S}_{ab} = F_a\hat{\phi}F_b^T$, where

Cross-correlation measurement uncertainty from PTA specifications

$\hat{\phi}$ is the scaled Fourier covariance matrix of the GWB (i.e., amplitude and ORF independent).

We use the Woodbury matrix lemma to perform inversions, such that

Hence the numerator of Equation (1) can be written as

where

and the denominator is written as

where

Crucially, the $X$ and $Z$ matrices for each pulsar depend on measured noise characteristics (e.g., white noise and intrinsic red noise) or fixed aspects of modeling (e.g., timing-model and Fourier design matrices). They can thus be constructed and cached for follow-up analysis. The diagonal matrix $\hat{\phi}$ acts as a spectral template that we use in a procedure akin to matched filtering, where we assess the S/N of cross-correlated signals with, e.g., a power-law PSD that has a $-13/3$ exponent, matching expectations for a circular, GW-driven population of SMBHBs forming a stochastic GWB.

Thus, in production-level PTA searches for stochastic GWB, the cross-correlations are computed as

Cross-correlation measurement uncertainty from PTA specifications

As shown in § 3, the trace in Equation (1) can be written as

where $T$ is the timing baseline; $A$ is the amplitude of the GWB; $P_g(f)$ is the power spectral density of timing residuals induced by the GWB; $P_a(f)$ and $P_b(f)$ are the intrinsic power spectral density of timing residuals in pulsars $a$ and $b$; and $f_l = 1/T$ and $f_h$ are the low and high frequency cutoffs used in the GWB analysis. For a GWB with characteristic strain spectrum

and using the convention from Siemens et al. [58], the GWB spectrum can be written in the timing residual space as

where $\gamma= 3 - 2\alpha$, and $f_{\rm ref}$ is a reference frequency. Consequently, the uncertainty on the cross correlations can be written as

As described in Siemens et al. [58], assuming no intrinsic red noise in the pulsars, and all pulsars are identical, the intrinsic power can be written as the sum of the GWB power and the white noise, $w$, in each pulsar, $P_a(f) = P_b(f) = P_g(f) + 2w^2\Delta t$, where the cadence of pulsar observations is given by $c = 1/\Delta t$. Substituting the intrinsic and GWB power in Equation (34), the cross correlation uncertainty can be written as

The solution for the integral in Equation (35) is nontrivial, but as described in Siemens et al. [58] and Chamberlin et al. [15], it can be solved using hypergeometric functions. We can define three regimes in which PTAs can operate: (i) the weak-signal regime, where the intrinsic pulsar noise dominates over the GWB power, i.e. $2w^2\Delta t \gg b f^{-\gamma}$; (ii) the strong-signal regime, where the GWB power dominates the intrinsic white noise, i.e. $2w^2\Delta t \ll b f^{-\gamma}$; (iii) the intermediate-signal regime where only the lowest few frequency bins in the PTA are dominated by the GWB, while at higher frequencies, the white noise dominates the GWB power. Given the recent results from the regional PTAs ([9, 16, 25]) and the IPTA ([7]), currently PTAs are in the intermediate-signal regime.

In the weak-signal regime, Siemens et al. [58] showed that Equation (35) can be written as

which reduces to

Similarly, in the strong-signal regime, the same integral reduces to (Siemens et al. [58])

Note that the scaling expressions above include the GWB amplitude, while in the rest of this paper, we use amplitude-scaled cross correlation values and uncertainties.