Detecting Gravitational Wave Anisotropies from Supermassive Black Hole Binaries

Abstract

Several Pulsar Timing Array (PTA) collaborations have recently found evidence for a gravitational wave background (GWB) permeating our galaxy. The origin of this background is still unknown. Indeed, while the gravitational wave emission from inspiraling supermassive black hole binaries (SMBHB) is the primary candidate for its origin, several cosmological sources have also been proposed. One distinctive feature of SMBHB-generated backgrounds is the presence of GWB anisotropies stemming from the binaries distribution in the local Universe. However, none of the anisotropy searches performed to date reported a detection. In this work, we show that the lack of anisotropy detection is not currently in tension with a SMBHB origin of the background. We accomplish this by calculating the anisotropy detection probability of present and future PTAs. We find that a PTA with the noise characteristics of the NANOGrav 15-year data set had only a $~{}2\%-11\%$ probability of detecting SMBHB-generated anisotropies, depending on the properties of the SMBHB population. However, we estimate that for the IPTA DR3 data set these probabilities will increase to $~{}4\%-28\%$, putting more pressure on the SMBHB interpretation in case of a null detection. We also identify SMBHB populations that are more likely to produce detectable levels of anisotropies. This information could be used together with the spectral properties of the GWB to characterize the SMBHB population.

I Introduction

Several Pulsar Timing Array (PTA) collaborations have recently found strong evidence for the presence of a gravitational wave background (GWB) in the nHz frequency band [1, 2, 3]. The most plausible source for such a background is a population of supermassive black hole binaries (SMBHB) formed in galaxy mergers through cosmic history. However, cosmological sources (including cosmic inflation, scalar-induced GW, phase transitions, cosmic strings, and domain walls) have also been proposed as a possible source of the background (see, for example, Ref. [4]). Discriminating between these two possible interpretations is one of the main near-term goals of PTA experiments.

Searches for anisotropies in the GWB provide a promising path to establishing the origin of the GWB. Indeed, if the GWB is produced by a population of inspiraling SMBHB, clustering of host galaxies and Poisson fluctuations in binaries properties are expected to induce a detectable level of anisotropies in the GWB [5]. On the other hand, the level of anisotropies produced by most cosmological sources is well below the sensitivity of present and future PTAs [6, 7]. Therefore, detecting anisotropies in the GWB would provide strong evidence in favor of an astrophysical origin of the signal. At the same time, we might use the lack of GWB anisotropies to constrain the properties of the SMBHB population and eventually rule out the astrophysical interpretation in favor of a cosmological one.

The recent search for anisotropies in NANOGrav 15-year data set has reported null evidence, and placed an upper limit on the level of broadband anisotropies [8]. In view of these results a few questions arise: Is the lack of evidence for GWB anisotropies in tension with an SMBHB origin of the background? What can we learn about the SMBHB population given this null detection? Are the properties of the SMBHB population needed to reproduce the intensity of the observed signal compatible with a null anisotropy detection?

In this work, we address these questions by estimating the anisotropy detection probability for GWBs sourced by SMBHB populations whose demographic and evolution are consistent with the background spectral properties. We do this by running realistic anisotropy searches on mock PTA data in which we injected the GW signal generated from synthetic SMBHB populations produced using holodeck [9]. Specifically, we proceed as follows: for a given set of SMBHB population parameters, we generate several mock SMBHB populations and their associated PTA signal, specifically the pulsar-pair correlations that such populations would induce (see Sec. III). We then run a frequentist anisotropy search [10] on these mock pulsar-pairs correlations, and derive the anisotropy detection probability for a given set of population parameters (see Sec. II and Sec. IV for details). Finally, we discuss our results (see Sec. V), which consist in a mapping between sets of SMBHB population parameters and anisotropies detection probabilities. This mapping gives us the tools to interpret the results of present and future searches for anisotropies. Specifically, it allows us to assess if a null anisotropy detection can be considered in tension with an SMBHB origin of the background, and use the results of anisotropy searches to better characterize the properties of the SMBHB population.

II How to Search for Anisotropies

In this section, we provide a summary of how GWB anisotropies appear in PTA data and discuss some of the data analysis techniques commonly adopted to extract these signals.

II.1 Anisotropies and pulsar-pair correlations

PTAs track the arrival times of radio pulses emitted by a collection of galactic millisecond pulsars. The presence of a GWB will appear in PTA data as correlated perturbations, $\delta t$, in the arrival times of these pulses. Specifically, let’s consider a GWB whose plane-wave expansion is given by

$$h_{ij}(t,\bm{x})=\sum_{A}\int df\int_{S^{2}}d\hat{\Omega}\;h_{A}(f,\hat{\Omega})e^{i2\pi f(t-\hat{\Omega}\cdot\bm{x})}e_{ij}^{A}(\hat{\Omega})\,,$$ (1)

where $f$ is the GW frequency, $A=+,\times$ labels the two GW polarizations, $e_{ij}^{A}$ are the GW polarization tensors, and $\hat{\Omega}$ identifies the propagation direction for a single GW plane wave. For a stationary, Gaussian, and unpolarized background the polarization amplitudes, $h_{A}(f,\hat{\Omega})$, are complex random variables with vanishing expectation values and a two-point correlator given by:

$$\langle h_{A}^{*}(f,\hat{\Omega})h_{A^{\prime}}(f^{\prime},\hat{\Omega}^{\prime})\rangle=\frac{1}{16\pi}\delta_{AA^{\prime}}\delta(f-f^{\prime})\delta(\hat{\Omega},\hat{\Omega}^{\prime})S_{h}(\hat{\Omega},f)\,.$$ (2)

The GWB power spectrum can be factorized as $S_{h}(\hat{\Omega},f)=S_{h}(f)P(\hat{\Omega},f)$, where the function $S_{h}(f)$ describes the spectral content of the GWB, and $P(\hat{\Omega},f)$ describes the distribution of the GWB power in the sky.

The expected correlations between the perturbations in the pulse time of arrivals, usually referred to as timing residuals, for pulsar $a$ at time $t_{i}$ and the timing residials for pulsar $b$ at time $t_{j}$ induced by such a GWB are given by [11]:

$$\langle\delta t_{a}(t_{i})\delta t_{b}(t_{j})\rangle=\int\frac{df}{24\pi^{2}f^{2}}\,e^{2\pi if(t_{i}-t_{j})}\,S_{h}(f)\,\Gamma_{ab}(f)\,.$$ (3)

From this equation is clear that $S_{h}(f)$ encodes the time correlation of the residuals. The spatial correlations of the residuals are instead encoded in the overlap reduction function, $\Gamma_{ab}$, which is related to the sky-distribution of the GWB power as:

$$\Gamma_{ab}(f)=\frac{3}{2}\int_{S^{2}}\frac{d\hat{\Omega}}{4\pi}\,P(\hat{\Omega},f)\sum_{A=+,\times}\mathcal{F}^{A}(\hat{p}_{a},\hat{\Omega})\mathcal{F}^{A}(\hat{p}_{b},\hat{\Omega})$$ (4)

where $\mathcal{F}^{A}(\hat{p},\hat{\Omega})$ is the antenna response function for the $A$-th GW polarization:

$$\mathcal{F}^{A}(\hat{p},\hat{\Omega})=\frac{\hat{p}^{i}\hat{p}^{j}}{2(1-\hat{\Omega}\cdot\hat{p})}e_{ij}^{A}(\hat{\Omega})\,.$$ (5)

Working with a discrete set of frequencies, indexed by $n$, and rewriting the integral in Eq. (4) as a sum over equal-areal pixels, indexed by $k$, we get:

$$\Gamma_{ab,n}=\frac{3}{2}\sum_{k}P_{n,k}\left[\mathcal{F}_{a,k}^{+}\mathcal{F}_{b,k}^{+}+\mathcal{F}_{a,k}^{\times}\mathcal{F}_{b,k}^{\times}\right]\frac{\Delta\hat{\Omega}_{k}}{4\pi}\,.$$ (6)

In this work, we will normalize the GWB power such that $\int d\hat{\Omega}\,P(\hat{\Omega},f_{n})=\sum_{k}P_{n,k}\Delta\hat{\Omega}_{k}=4\pi$. With this normalization, a perfectly isotropic GWB is given by $P_{n,k}=1$, and when plugged into Eq. (6) it returns the familiar Hellings-Downs correlation pattern [13].

II.2 Anisotropy searches

Given a set of timing residuals, we can estimate the cross-correlation coefficients using the frequentist estimator developed in Refs. [14, 15, 16, 17]:

$$\rho_{ab}=\frac{\delta\mathbf{t}_{a}^{T}\mathbf{N}_{a}^{-1}\hat{\mathbf{S}}_{ab}\mathbf{N}_{b}^{-1}\delta\mathbf{t}_{b}^{T}}{\,{\rm tr}\left[\mathbf{N}_{a}^{-1}\hat{\mathbf{S}}_{ab}\mathbf{N}_{b}^{-1}\hat{\mathbf{S}}_{ba}\right]},$$ (7)

where $\rho_{ab}$ are the best estimates for the cross-correlation coefficients, $\delta\mathbf{t}_{a}$ contains the timing residuals for pulsar $a$, $\mathbf{N}_{a}=\langle\delta\mathbf{t}_{a}\delta\mathbf{t}_{a}^{T}\rangle$ is the $a$-th pulsar’s autocovariance matrix, and $\hat{\mathbf{S}}_{ab}$ is the template-scaled covariance matrix. This template-scaled covariance matrix encodes information about the spectral shape of the GWB but is agnostic about its overall amplitude and cross-correlation signature. It is related to the full covariance matrix, $\mathbf{S}_{ab}\equiv\langle\delta\mathbf{t}_{a}\delta\mathbf{t}_{b}^{T}\rangle$, as $\mathbf{S}_{ab}=A^{2}\Gamma_{ab}\hat{\mathbf{S}}_{ab}$, where $A$ is the GWB amplitude for a given spectral shape template. The values of $\mathbf{N}_{a}$ and $\hat{\mathbf{S}}$ are not known a priori, and need to be extracted from the data, usually using Bayesian inference techniques.

Equation (7) allows us to estimate broad-band correlations for the timing residuals. However, the GWB sky sourced by realistic SMBHB populations will not be constant across frequencies, and therefore induce correlations which are not constant in frequency space. These frequency-resolved correlations can be estimated using the recently developed per-frequency Optimal Statistic (PFOS) methods described in Ref. [18]. Assuming stationary Gaussian noise for the cross-correlation uncertainties, the likelihood function for these frequency-resolved correlations, $\bm{\rho}_{n}$, given a GWB sky, $\mathbf{P}_{n}$, can be written as:

$$p(\bm{\rho}_{n}|\mathbf{P}_{n})=\frac{\exp[-\frac{1}{2}(\bm{\rho}_{n}-\mathbf{R}\mathbf{P}_{n})\mathbf{\Sigma}^{-1}_{n}(\bm{\rho}-\mathbf{RP}_{n})]}{\sqrt{\det(2\pi\mathbf{\Sigma}_{n})}},$$ (8)

where $\bm{\Sigma}_{n}$ is the $N_{\text{cc}}\times N_{\text{cc}}$ diagonal noise covariance matrix for the cross-correlation measurements in the $n$-th frequency bin, and we have defined the PTA overlap response matrix

$$R_{ab,k}=\frac{3}{2N_{pix}}[\mathcal{F}_{a,k}^{+}\mathcal{F}_{b,k}^{+}+\mathcal{F}_{a,k}^{\times}\mathcal{F}_{a,k}^{\times}].$$ (9)

Therefore, given a set of measured cross-correlations, we can maximize this likelihood to obtain an estimate of the GWB sky, $\hat{\mathbf{P}}_{n}$. In this work, we choose to work in the pixel basis [20], in which the GWB sky is divided into independent equal-area pixels. In the pixel basis, reconstructing the GWB sky corresponds to reconstructing the GWB power in each of these equal area pixels. Other typical basis choices include decomposition of the sky into linear [21, 22, 23, 24] and square root spherical harmonics [25] or eigenmaps of the PTA Fisher matrix [26].

In the rest of this section, we discuss the three techniques used in this work to maximize the likelihood function in Eq. (8). Examples of sky maps reconstructed by using some of these techniques are shown in Fig. 1.

II.2.1 Analytical Solution

In the pixel basis, there exists an analytical maximum solution for the likelihood in Eq. (8), given by [27, 28, 29]:

$$\mathbf{\hat{P}}_{n}=\mathbf{M}_{n}^{-1}\mathbf{X}_{n},$$ (10)

where we have defined the Fisher information matrix, $\mathbf{M}_{n}=\mathbf{R}^{T}\mathbf{\Sigma}^{-1}_{n}\mathbf{R}$, and the “dirty map”, $\mathbf{X}_{n}=\mathbf{R}^{T}\mathbf{\Sigma}^{-1}_{n}\bm{\rho}_{n}$. The diagonal elements of the Fisher matrix provide an estimate of the uncertainty on the recovered pixel values, while the dirty map provides an inverse noise-weighted representation of the power as seen by the PTA.

The inversion of the Fisher matrix typically limits the resolution of our maximum likelihood sky maps beyond the fundamental PTA limit. The spatial resolution of a PTA is approximately set by requiring $N_{\text{pix}}\leq N_{\text{cc}}$ [28], where $N_{cc}$ is the number of cross-correlations in the data set, and $N_{\rm pix}$ is the number of pixels in the sky-map tessellation. For a HEALpix [30] sky map, $N_{\rm pix}$ is related to the $N_{\text{side}}$ parameter as $N_{\text{pix}}=12\cdot N_{\text{side}}$. Given the NANOGrav 15yr configuration with 67 pulsars and $N_{\text{cc}}=2211$, this translates to $N_{\text{side}}\leq 8$. However, we find the inverse of the Fisher matrix to be numerically ill-defined for $N_{\text{side}}=8$, which limits the resolution of the NANOGrav 15-year data set to $N_{\textrm{side}}\leq 4$, corresponding to 12 ($N_{\textrm{side}}=1)$, 48 ($N_{\text{side}}=2$) and 192 ($N_{\text{side}}=4$) pixels for the sky maps.

Another downside of this map-reconstruction method is that it usually produces maps with negative power in some regions of the sky. These negative-power regions are –of course– unphysical and prevent us from normalizing the recovered maps.

II.2.2 Radiometer Search

If we assume that most of the GW power is coming from a single bright pixel, instead of trying to reconstruct the entire GW sky, we can focus on reconstructing the power in each pixel independently, neglecting correlations between neighbouring pixels. The power in each pixel can be derived by replacing the Fisher matrix with its diagonal components in Eq. (10).

This reconstruction method, usually known as radiometer [31, 32], is best suited for GWB skies that contain a single, very bright hot spot instead of a collection of sources with similar signal amplitudes. While some SMBHB populations might generate GW skies with very bright hot spots, it is not clear a priori if a radiometer search is best suited to recover generic SMBHB skies. As part of our analysis, we will assess how well the radiometer search performs on realistic SMBHB skies.

II.2.3 Numerical Solution

Finally, we can maximize the cross-correlation likelihood numerically. Specifically, we can reconstruct the GW sky by numerically solving the following equation

$$\mathbf{M}_{n}\mathbf{\hat{P}}_{n}=\mathbf{X}_{n}\,.$$ (11)

The advantage of this method is that it allows us to impose positivity on the recovered sky maps and therefore interpret each pixel value as physical power. Additionally, this approach does not rely on our ability to calculate the inverse of $\mathbf{M}_{n}$.

III Mock Data

In this section we describe the procedure we used to generate mock PTA data. While, in principle, we could produce a set of mock timing residuals for each synthetic SMBHB population and then derive the corresponding pulsar cross-correlations, we instead directly generate the cross-correlation coefficients, $\rho_{ab,n}$, to maximize the computational efficiency of our procedure. Specifically, we proceed by following these three main steps:

Population Synthesis (Sec. III.1) –

We explore the parameter space describing formation and evolution of SMBHBs. For each point in this space, we generate several realizations of the cosmic SMBHB population using the holodeck package [9].

Sky Creation (Sec. III.2) –

For each SMBHB population generated in the previous step, we create a GWB sky, $P_{n,k}$, by randomly placing the population binaries in the sky.

Correlation Creation (Sec. III.3) –

For each GWB sky, we generate theoretical pulsar-pair correlation according to Eq. (6). To simulate measurement noise, on top of these theoretical correlations, we add zero-mean Gaussian noise whose standard deviation is given by the cross-correlations uncertainty measured from the NANOGrav 15-year data set [8].

III.1 Population Synthesis with holodeck

To generate SMBHB populations, we used the same procedure as in Ref. [35], which we briefly summarize in this section. For more details, we refer the reader to Refs. [35, 9].

Model Component	Parameter	Fiducial value
GSMF [36]	$\psi_{z}$	-0.60
	$m_{\psi,z}$	+0.11
	$\alpha_{\psi,0}$	-1.21
	$\alpha_{\psi,z}$	-0.03
$M_{\text{BH}}-M_{\text{bulge}}$ [37, 38, 39, 40, 41]	$\alpha_{\mu}$	1.10
	$f_{\star,\text{bulge}}$	0.615
$\dfrac{da}{dt}$ [42, 35]	$a_{c}/{\rm pc}$	$10^{2}$
	$a_{\text{init}}/{\rm pc}$	$10^{4}$
	$\nu_{\text{outer}}$	2.5

Table 1: Fiducial model parameters as presented in table B1 of Ref. [35].

III.1.1 Galaxy Merger Rate

As SMBHB are formed in galaxy mergers, one of the key ingredients needed to model SMBHB populations is the rate at which such mergers take place. The comoving number density of galaxies mergers, $\eta_{\rm gal-gal}$, as function of the stellar mass of the primary galaxy, $m_{1\star}$, the stellar mass ratio, $q_{\star}\equiv m_{2\star}/m_{1\star}$, and redshift, $z$, can be written as [36]:

$$\frac{\partial^{3}\eta_{\text{gal-gal}}}{\partial m_{\star 1}\partial q_{\star}\partial z}=\frac{\Psi(m_{\star 1},z^{\prime})}{m_{\star 1}\ln(10)}\frac{P(m_{\star 1},q_{\star},z^{\prime})}{T_{\text{gal-gal}}(m_{\star 1},q_{\star},z^{\prime})}\frac{\partial t}{\partial z^{\prime}},$$ (12)

where $\Psi$ denotes the galaxy stellar-mass function (GSMF), $P$ the galaxy pair fraction and $T_{\text{gal-gal}}$ the galaxy merger time. Since the galaxy merger spans a finite timescale ($T_{\rm gal-gal}$), in Eq. (12), we distinguish between the redshift at which the galaxy pair forms ($z^{\prime}=z^{\prime}[t]$ at some initial time $t$) and the redshift at which the system becomes a post-merger galaxy remnant ($z=z[t+T_{\rm gal-gal}]$).

The GSMF provides the differential number-density of galaxies per decade of stellar mass. Following Ref. [35], we parametrize the GSMF as a single Schechter function [43]:

$$\Psi(m_{\star 1},z^{\prime})=\ln(10)\Psi_{0}\cdot\bigg{(}\frac{m_{\star 1}}{M_{\psi}}\bigg{)}^{\alpha_{\psi}}\exp\bigg{(}-\frac{m_{\star 1}}{M_{\psi}}\bigg{)},$$ (13)

where the phenomenological parameters $\Psi_{0}$, $M_{\psi}$ and $\alpha_{\psi}$ follow the redshift parametrization introduced in [36]:

$$\begin{split}\log_{10}(\Psi_{0}/\,{\rm Mpc}^{-3})&=\psi_{0}+\psi_{z}\cdot z\\ \log_{10}(M_{\psi}/M_{\odot})&=m_{\psi,0}+m_{\psi,z}\cdot z\\ \alpha_{\psi}&=1+\alpha_{\psi,0}+\alpha_{\psi,z}\cdot z\end{split}$$ (14)

This parametrization introduces six free parameters for the GSMF. Of these six parameters, the mass normalization $\psi_{0}$ and the characteristic mass $m_{\psi,0}$ are allowed to vary in our analysis, while the remaining four are fixed to the fiducial values provided in Table 1.

The galaxy pair fraction and merger time are kept fixed to fiducial power-law functions of the galaxy stellar mass, galaxy mass ratio, and initial redshift. For a full description of these functions we direct the reader to Ref. [35].

III.1.2 SMBH-Host Relation

By populating galaxies with supermassive black holes (SMBHs), we can relate the galaxies merger rate in Eq. (12) with the SMBHs merger rate. We populate the galaxies assuming a correspondence between the host galaxy and the SMBH at its center. Specifically, we assume that the SMBH mass, $M_{\rm BH}$, is related to the bulge mass of the host galaxy, $M_{\rm bulge}$, according to the relation [44]:

$$\log_{10}(M_{\text{BH}}/M_{\odot})=\mu+\alpha_{\mu}\log_{10}\bigg{(}\frac{M_{\text{bulge}}}{10^{11}M_{\odot}}\bigg{)}+\mathcal{N}(0,\epsilon_{\mu})\,$$ (15)

where $\mathcal{N}(0,\epsilon_{\mu})$ is a normally distributed random scatter with zero mean and standard deviation $\epsilon_{\mu}$, and $M_{\rm bulge}=f_{\star,{\rm bulge}}\cdot m_{\star}$ is the fraction of the galaxy stellar mass contained in the bulge. In our analysis, we vary the mass normalization $\mu$ and the scattering width $\epsilon_{\mu}$, while keeping the other parameters fixed to the fiducial values reported in Table 1.

Using the $M_{\rm BH}-M_{\rm bulge}$ relation we can then translate the number density of galaxy mergers into a number density for SMBHB mergers, $\eta$, as:

$$\frac{\partial^{3}\eta}{\partial M\partial q\partial z}=\frac{\partial^{3}\eta_{\text{gal-gal}}}{\partial m_{\star 1}\partial q_{\star}\partial z}\frac{\partial m_{\star 1}}{\partial M}\frac{\partial q_{\star}}{\partial q}\,,$$ (16)

where $M=m_{1}+m_{2}$ is the total binary mass, $q\equiv m_{2}/m_{1}<1$ is the binary mass ratio, and $m_{1,2}$ are the masses of the SMBHs in the binary system.

Binary property	Range	Number of bins
$M/M_{\odot}$	$[10^{4},10^{12}]$	91
$q$	$[10^{-3},1]$	81
$z$	$[10^{-3},10]$	101
$f/\rm{nHz}$	$[1.85,5.93]$	3

Table 2: Details of the SMBHB parameter space binning. Frequency bins are linearly spaced, while the bins for all the other dimensions are logarithmically spaced.

III.1.3 Binary Evolution

For binaries on circular orbits, the observer-frame GW frequency, $f$, is related to the binary’s rest frame orbital frequency, $f_{p}$, as $f=2f_{p}/(1+z)$. The orbital frequency, in turn, is uniquely related to the binary orbital separation, $a$, according to $f_{p}^{2}=4\pi GM/a^{3}$.

The evolution of the binary system’s orbital separation, $a$, is driven by a combination of GW emission and interactions with the astrophysical environment. The GW contribution to the binary hardening is given by [45]:

$$\frac{da}{dt}\bigg{|}_{\text{GW}}=-\frac{64G^{3}}{5c^{5}}\frac{m_{1}m_{2}M}{a^{3}}\,,$$ (17)

where $G$ is the gravitational constant. We describe environmental effects using a phenomenological model which aims to capture the overall effects of more complicated explicit binary evolution models, while introducing a small number of free parameters:

$$\frac{da}{dt}\bigg{|}_{\text{phenom}}=H_{a}\bigg{(}\frac{a}{a_{c}}\bigg{)}^{1-\nu_{\text{inner}}}\bigg{(}1+\frac{a}{a_{c}}\bigg{)}^{\nu_{\text{inner}}-\nu_{\text{outer}}}.$$ (18)

Here, the power law indices $\nu_{\text{inner}}$ and $\nu_{\text{outer}}$ are introduced, dividing the hardening process into a small-separation and a large-separation regime with different power-law behaviours and a turnover at a critical separation $a_{c}$. While $\nu_{\text{inner}}$ is allowed to vary in our model, $\nu_{\text{outer}}$ is fixed to a standard literature value of $+2.5$ [42]. The total rate of binary evolution is obtained by summing the two contributions:

$$\frac{da}{dt}=\frac{da}{dt}\bigg{|}_{\text{GW}}+\frac{da}{dt}\bigg{|}_{\text{phenom}}\,.$$ (19)

The normalization of the phenomenological hardening rate, $H_{a}$, is determined by requiring that the binary lifetime, from the initial separation $a_{\rm init}$ to the innermost stable circular orbit $a_{\rm isco}=6GM/c^{2}$, matches a target value given by:

$$\tau_{f}=\int_{a_{\text{init}}}^{a_{\text{isco}}}\bigg{(}\frac{da}{dt}\bigg{)}^{-1}da\,.$$ (20)

The value of $\tau_{f}$ is a free parameter that we allow to vary in our analysis.

Given a model for the binary orbital evolution, we can write the number of SMBHB per given mass, mass-ratio, redshift, and log orbital frequency, as [46, 47, 48]:

$$\frac{\partial^{4}N}{\partial M\partial q\partial z\partial\ln f_{p}}=\frac{\partial^{3}\eta}{\partial M\partial q\partial z}\cdot\tau(f_{p})\cdot 4\pi c(1+z)d_{c}^{2}\,,$$ (21)

where $d_{c}$ is the comoving distance of a source at redshift $z$, and $\tau(f_{p})$ is the binary hardening timescale, which characterizes the amount of time a SMBHB spends in a given frequency bin:

$$\tau(f_{p})\equiv f_{p}/(df_{p}/dt)=-\frac{2a}{3}\frac{da}{dt}\,.$$ (22)

III.1.4 Population Synthesis

Given a set of astrophysical parameters, $\bm{\theta}=\{\psi_{0},\,m_{\psi,0},\,\mu,\,\epsilon_{\mu},\,\tau_{f},\,\nu_{\rm inner}\}$, we generate random SMBHB populations by drawing the number of binaries in a given $(M,\,q,\,z,\,f)$ bin from a Poisson distribution centered around the expectation value

$$N(M,q,z,f_{p})=\frac{\partial^{4}N}{\partial M\partial q\partial z\partial\ln f_{p}}\Delta M\,\Delta q\,\Delta z\,\Delta\ln f_{p},$$ (23)

where $(\Delta M,\,\Delta q,\,\Delta z,\,\Delta f_{p})$ are the bin sizes for each SMBHB parameter. Details about the binning of the SMBHB parameter space are reported in Tab. 2.

To ensure that the SMBHB populations considered in our analysis are consistent with the spectral properties of the observed GWB, we sample the astrophysical parameters from the posterior distributions derived by fitting the SMBHB signal to the observed GWB spectrum. Specifically, we sample from the Monte Carlo chains obtained for the Phenom+Astro model in Ref. [35]. For each $\bm{\theta}$ in these chains, we then generate 50 independent SMBHB populations.

III.2 GWB Sky Maps

Given a mock SMBHB population, the next step consists in generating an associated GW sky. We start by assuming that each of the SMBHB is on a circular orbit, and assigning them the $M$, $q$, $z$, and $f$ values corresponding to their bin centers. We can then derive their sky and polarization-averaged GW amplitude as [49]:

$$h_{s}^{2}(f)=\frac{32}{5c^{8}}\frac{(G\mathcal{M})^{10/3}}{d_{c}^{2}}(2\pi f_{p})^{4/3}\,,$$ (24)

where $\mathcal{M}=Mq^{3/5}/(1+q)^{6/5}$, and $f_{p}=f(1+z)/2$. This amplitude can then be translated into a characteristic strain as [50]:

$$h_{c}^{2}(f)=h_{s}^{2}(f)\frac{f}{\Delta f}\,,$$ (25)

where the frequency bin width is given by $\Delta f=1/T_{\text{obs}}$.

For each frequency bin, we then place the loudest 20 sources on random pixels of a sky-tessellation derived using HEALpix [30], and add their squared characteristic strain to the corresponding pixels. The remaining binaries are added to the sky map as an isotropic background (equally divided among all pixels) with a characteristic strain given by

$$h^{2}_{c,\scriptscriptstyle{\rm BG}}(f)=\sum_{M,q,z,f_{p}}N(M,q,z,f_{p})h_{c}^{2}(f)\,,$$ (26)

where the sum runs over all binaries expect the loudest 20. The maps are then normalized such that $\sum_{k}P_{n,k}\Delta\hat{\Omega}_{k}=4\pi$, such that $P_{n}=1$ describes an isotropic sky and a theoretical upper limit of $P_{n,k,\rm max}=4\pi/\Delta\hat{\Omega}=N_{pix}$ is imposed for the power in each pixel, which corresponds to a singular GW source being present.

To marginalize over the sky locations of the loudest sources, for each SMBHB population, we generate 30 independent sky maps per frequency bin. Given that we generate 50 SMBHB populations for each set of astrophysical parameters (see the previous section for details), we produce a total of 1500 sky maps per frequency bin for each set of astrophysical parameters. Examples of mock sky maps are provided in the upper panels of Fig. 1.

Figure 1: Examples of GWB power distributions induced by SMBHB populations (upper panel), and the reconstructed sky-maps obtained using the analytical (middle panel) and numerical (lower panel) reconstruction methods. The recovered pixel values are rescaled by their corresponding uncertainties derived from a null distribution of recovered sky maps from isotropic GW skies, $\sigma_{P}$, indicating the significance of deviation from isotropy. The left panels illustrate the reconstruction of a GW sky with a hot spot that leads to an anisotropy detection in our analysis. The right panels show a GW sky that does not produce a detection.

Figure 2: Pixel upper limits derived for $f_{\rm GW}=3.95$ nHz and $N_{\text{side}}=2$ using the analytical (left panel) and numerical (right panel) reconstruction methods. The red stars indicate the position of the pulsars contained in the NANOGrav 15-year data set. The higher pixel values for the analytical upper limits reflect our inability to normalize the reconstructed sky maps due to the presence of negative pixel values.

III.3 Mock Correlations

The final step of the procedure consists of generating mock cross-correlations for each of the mock GWB skies. For a given sky map, $P_{i,k}$, the expected cross-correlations values are given by Eq. (6). However, any realistic measurement will not return these theoretical values due to finite experimental precision. Therefore, to simulate a realistic anisotropy search, we generate mock correlation values by sampling from a Gaussian distribution centered around the expected value given by Eq. (6), and with a standard deviation given by the cross-correlation uncertainties measured from the NANOGrav 15-year data set using the PFOS [18].

We use this procedure to generate two mock PTA data sets. The first data set consists of the cross-correlations values for the $67$ pulsars monitored for the NANOGrav 15-year data release. For the second data set, we try to mimic the capabilities of a future IPTA and generate mock cross-correlations for the 115 pulsars currently observed by NANOGrav [52] together with EPTA+InPTA [53], PPTA [54], and MeerKat [55]. Not having a realistic noise estimate for this data set, we use the cross-correlation uncertainties derived from the NANOGrav 15-year data set for the subset of pulsar pairs already included in this data set, and draw the noise values randomly from this set of cross-correlation uncertainties for pulsar pairs not included in the current NANOGrav observations.

It is important to note that, even for a PTA with infinite resolution, cosmic variance, polarization of the GW emission, and contributions from the pulsar term will induce deviations from the expected cross-correlations in Eq. (6) (see, for example, Ref. [56]). While these deviations are below the current precision at which we can reconstruct the cross-correlations, their impact on our ability to reconstruct the underlying GWB sky has not yet been investigated. In this work, we will follow standard conventions of anisotropy searches, and ignore the impact of these effects. We leave a dedicated study of their impact to a follow-up analysis [12].

IV Detection Statistics

In the final step of our procedure, we take the mock cross-correlations generated for each SMBHB population and try to reconstruct the underlying GW sky by using the techniques discussed in Sec. II. For each of these reconstructed skies, we want to establish if they show any evidence for anisotropies. We do this by constructing an appropriate detection statistic that we then calibrate against a null distribution.

When searching for anisotropies, we assume isotropy as the null hypothesis. Therefore, we calibrate all our detection statistics against a null distribution derived by re-running our pipeline on mock cross-correlations generated from perfectly isotropic skies (instead of the SMBHB populated skies used in Sec. III). Specifically, we generate null distributions by proceeding as follows:

• For each pulsar pair, we generate a cross-correlation value drawing from a normal distribution with a mean given by the Hellings & Downs (HD) value and a standard deviation given by the cross-correlation uncertainties measured from the NANOGrav 15-year data set using the PFOS estimator [18]. We note that, even for isotropic skies, the cosmic variance effect already discussed in the previous section will induce deviations from HD correlations. In our analysis, following standard conventions (see for example Refs. [21, 24, 8]), we ignore this effect. The impact of this assumption will be investigated in a follow-up analysis [12].

• By using the methods discusses in Sec. II, we recover the GWB sky associated to these mock correlations.

• We compute the value of the detection statistic on the recovered GWB sky.

• We repeat this procedure $10^{6}$ times to derive a null distribution for the detection statistic.

By comparing the detection statistic values of a reconstructed sky with this null distribution we can assess the significance of the recovered anisotropies.

In the rest of this section, we discuss three detection statistics used in this work to define anisotropy detection.

IV.1 Pixel upper limits

An intuitive detection statistic is provided by the reconstructed GWB power in each pixel. By comparing the reconstructed pixel power with the pixel-specific null distribution, we can claim that a reconstructed sky map provides a detection with a $\sim 3\sigma$ global significance if any of the reconstructed pixel power has a $p$-value smaller than $3\times 10^{-3}/N_{\rm pix}$. Here, the trial factor $1/N_{\rm pix}$ converts the local significance of each pixel to a global significance.

The pixel threshold values for maps with $N_{\rm side}=2$, and for a GW frequency of $f_{\rm GW}=3.95\,{\rm nHz}$ are reported in Fig. 2. As expected, the lowest threshold values are found in the region of the sky with the highest density of pulsars.

IV.2 Anisotropic signal-to-noise ratio

Instead of defining a detection statistic based on the individual pixels, we can also use one that takes the full sky map into account. One such statistic is the anisotropic signal-to-noise ratio (SNR) [10], which is given by the maximum likelihood ratio between an anisotropic and an isotropic sky:

$${\rm SNR}=\sqrt{2\ln\bigg{[}\frac{p(\bm{\rho}|\mathbf{\hat{P}})}{p(\bm{\rho}|\mathbf{P}_{\text{iso}})}\bigg{]}}\,,$$ (27)

where $p$ is the likelihood function given in Eq. (8), $\mathbf{\hat{P}}$ is the recovered anisotropic sky map, and $\mathbf{P}_{\text{iso}}$ is the perfectly isotropic sky map. This detection statistics takes advantage of the full reconstructed sky, such that sources located in different regions of the sky can contribute to a classification as a detection (even if below the individual pixel thresholds discussed in the previous section).

Due to its global nature, this detection statistic is not suitable for assessing the significance of GW skies reconstructed using the radiometer approach. Indeed, the radiometer method only reconstructs the power of individual pixels while assuming that all the others do not contribute to the observed correlation pattern. Because of this, we will not use the SNR detection statistic in conjunction with the radiometer reconstruction method.

As for the pixel-based statistic, the SNR values are calibrated against null distributions to identify sky maps that give us a detection with a $p$-value $p<3\times 10^{-3}$ (corresponding to a $\sim 3\sigma$ Gaussian-equivalent significance).

IV.3 $\chi^{2}$-Statistic

Finally, we consider a third test statistic based on the “distance” between the reconstructed sky map and the “typical” recovered map for an underlying isotropic sky. Specifically, assuming that the uncertainties on the pulsar-pair correlations are normally distributed, the null distribution of the recovered pixel values also follows a multivariate normal distribution, $\mathcal{N}(\mathbf{\mu},\mathbf{\Sigma})$, where $\mathbf{\mu}$ and $\mathbf{\Sigma}$ are the mean and covariance of the recovered pixel values for the sky maps of the null distribution. We can than calculate a probability distance estimator, the Mahalanobis distance [57], for a specific sky map, $\mathbf{P}$, as

$$D^{2}=(\mathbf{P}-\mu)^{T}\mathbf{\Sigma}^{-1}(\mathbf{P}-\mu).$$ (28)

As for the previous test statistics, we compare the Mahalanobis distance for any recovered map with its null distribution. In this case the null distribution can be derived analytically, it follows a $\chi^{2}$ distribution with the number of degrees of freedom given by the number of pixels. Recovered skies that produce a Mahalanobis distance with a $p$-value smaller than $p<3\times 10^{-3}$ are considered as detections.

Our application of the Mahalanobis distance as a detection statistic assumes that the null distribution is given by a multivariate normal distribution. While this is the case for both the radiometer pixel basis and the analytical map reconstruction, the positivity constraint that we impose during the numerical reconstruction changes the shape of the null distribution such that we can not use the Mahalanobis distance as an estimator for anisotropy detection in combination with the numerical approach.

V Results

The main result of our analysis consists of an estimate of the anisotropy detection probability for different combinations of SMBHB population parameters. For each set of these parameters, $\bm{\theta}$, the anisotropy detection probability, $p_{\bm{\theta}}$, is estimated as

$$p_{\bm{\theta}}=\frac{N_{\bm{\theta},\rm det}}{N_{\rm tot}}\,,$$ (29)

where $N_{\rm tot}=1500$ is the total number of mock GWB skies generated for each set of SMBHB population parameters, and $N_{\bm{\theta},\rm det}$ is the fraction of those skies that gave us an anisotropy detection.

In Fig. 3, we report the posterior distributions of the SMBHB parameters used in this work. These distributions were derived in Ref. [35] by fitting the spectrum of SMBHB-sourced GWBs against the NANOGrav 15-year data. In this figure, the individual samples are colored according to their anisotropy detection probability, $p_{\bm{\theta}}$, for a NANOGrav-like PTA at $f_{\rm GW}=3.95\,{\rm nHz}$ (the most sensitive frequency bin). From this plot, we can identify some general trends in the detection probability values. The parameter that mostly impacts the detection probability is $\nu_{\rm inner}$, with larger values of this parameter leading to larger anisotropy detection probabilities (up to $p_{\bm{\theta}}\sim 0.45$ for $\nu_{\rm inner}\sim 0)$. This is consistent with the results of Ref. [51], where it was shown that larger $\nu_{\rm inner}$ lead to a faster evolution of the binaries at small separation and fewer sources contributing to the background in the lowest frequency bins. Averaging over the SMBHB population parameters, we find an average detection probability of $\bar{p}_{\bm{\theta}}\simeq 0.07$. This suggests that the null-detection for anisotropies reported in Ref. [8] is still perfectly compatible with an SMBHB interpretation of the GWB.

Figure 3: The off-diagonal panels show the 2D posterior distribution for the SMBHB population parameters derived in Ref. [35] by fitting the SMBHB signal to the observed GWB spectrum. The black contours represent the 68%, 95% and 99.7% Bayesian credible regions for these posterior distributions, while the dots show the MCMC samples from which the posterior distributions are derived. The individual MCMC samples are colored according to their estimated anisotropy detection probability, $p_{\bm{\theta}}$. The diagonal panels show the marginalized 1D posterior distribution, the lines are colored according to the marginalized anisotropy detection probability.

The observed trends in detection probabilities suggest that anisotropy detection (or lack thereof) provides additional information about the SMBHB properties in addition to the ones already encoded in the background spectral properties. This extra information could be used in future analyses to break the degeneracy between SMBHB populations leading to similar GWB spectra. Specifically, we envision an analysis in which the information provided by a positive (or negative) anisotropy detection is included in the likelihood used in Monte Carlo explorations of the SMBHB parameter space.

As expected, we find that most of the GW skies leading to a detection present a hot spot in the region of the sky in which we have the highest density of pulsars, and consequently a higher sensitivity to GWB anisotropies. This fact is clearly illustrated by Fig. 4, which reports the average of the detected GW skies, and shows how the typical detectable sky has more power in the region of the sky that is more densely populated by pulsars included in the data set. This fact also suggests that searches for anisotropies will greatly benefit from the more uniform sky coverage that will be provided by the upcoming IPTA DR3. Estimates for anisotropy detection probabilities in an IPTA-like data set are reported in Fig. 5. We find that, for this PTA, detection probabilities can reach up to $p_{\theta}\sim 0.85$ for $\nu_{\rm inner}\sim 0$, while the average detection probability grows to $\bar{p}_{\bm{\theta}}\sim 0.16$. This result illustrates how future anisotropy searches could provide the first evidence in support of SMBHB as the origin of the background, or –in case of a null detection– put this interpretation under pressure.

Figure 4: Average of the GWB skies (at $f_{\rm GW}=3.95\,{\rm Hz}$) that would lead to a detection of anisotropies in a PTA dataset with the noise properties of NANOGrav 15-year data set (using numerical map reconstruction and SNR-based detection statistic). The red stars indicate the position of the pulsars contained in the NANOGrav 15-year data set.

Figure 5: Same for Fig. 3 but for a PTA containing all the pulsars currently observed by EPTA+InPTA, MeerKAT, NANOGrav, and PPTA.

Figure 6: The left panel shows the average (histograms) and $68\%$ intervals (bars) for anisotropy detection probability (for $N_{\rm side}=2$) in the first three frequency bins and for different combinations of the map reconstruction method (numerical , analytical , or radiometer ) and detection statistic (SNR ,pixel-based , or $\chi^{2}$-based statistic ). The right panel shows the same quantities for the second frequency bin and different values of $N_{\rm side}$.

While in Fig. 3 and 5 we only report detection probabilities for the second frequency bin, we have also derived analogous results for the remaining lowest five frequency bins. We find that detection probabilities at other frequencies show similar trends in the SMBHB parameter space, even if with different overall normalizations. The left panel of Figure 6 reports the average detection probabilities in the first three frequency bins, and shows how detection probabilities are the largest in the second (we have also checked detection probabilities at higher frequencies, finding that sensitivity to anisotropies keeps degrading going to higher frequencies). The second frequency bin results the most likely to detect GW anisotropies from SMBHB due to the balancing of competing effects: the level of anisotropies from SMBHBs is expected to increase at higher frequencies due to decreasing numbers of binaries contributing to the signal at those frequencies [51]; however, the red-tilt of the GWB spectrum causes the cross-correlation uncertainties to strongly increase at higher frequencies, limiting the sensitivity to anisotropies to the lower frequency range. The inclusion of pulsars with a timing baseline significantly shorter than 15 years introduces additional uncertainties in the first frequency bin, which leaves the second bin as the one most sensitive to searches for anisotropies.

While in Fig. 3 and Fig. 5 we only report the results derived with a combination of numerical map reconstruction and SNR-based detection statistic, we also derived these detection probabilities using all possible combinations of sky reconstruction techniques (discussed in Sec. II) and detection statistics (discussed in Sec. IV). We find that for all these possible map reconstruction-detection statistic combinations we observe the same detection probability trends in the SMBHB parameter space, but with different overall normalization. Specifically, we find that numerical map reconstruction combined with the SNR-based detection statistics leads to the best detection prospects when coupled with a map resolution of $N_{\rm side}=2$ or $N_{\rm side}=4$ (see the right panel of Fig. 6). Except for the lower resolution maps ($N_{\rm side}=1)$, the analytical reconstruction tends to underperform the numerical reconstruction methods, while the radiometer basis delivers similar results. This shows how, despite realistic SMBHB skies presenting more than a single bright pixel, the radiometer basis is still able to resolve them.

Since most of the GW skies that lead to a positive detection for anisotropies are characterized by a loud binary, it is interesting to consider the interplay with continuous wave (CW) searches. By comparing the CW upper limits provided in Ref. [58] with the characteristic strain of all binaries in our mock GW skies, we find that roughly $35\%$ of the anisotropy detections in the second frequency bin would be accompanied by a positive detection in CW searches at the $2\sigma$-level (see the left panel of Fig. 7). Similarly, we find that around $40\%$ of the skies with a CW-detectable binary in the second frequency bin would also lead to a positive anisotropy detection (see the right panel of Fig. 7). Of course, these results are only rough estimates based on marginalized CW upper limits. A more detailed analysis would require performing a full CW search on the SMBHB populations used to derive the mock GW skies. We leave this analysis to future studies [59].

Figure 7: Left: Fraction of anisotropy detections above the CW upper limits. Right: Fraction of CW detectable skies that also led to an anisotropy detection.

Finally, we want to emphasize that in deriving our results we ignored the impact that interference between the binaries signals can have on the pulsars cross correlations [56, 19]. This effect –which is a combination of what is commonly known as pulsar variance and cosmic variance– can induce significant deviations from the expected cross-correlation values given in Eq. (4) and make anisotropy detection more challenging. The precise impact of these effects, together with possible mitigation strategies, will be discussed in future work [12]. However, because of this approximation, the detection probabilities reported herein should be regarded as optimistic estimates, contingent upon improvements in detection methods to account for such interference effects.

VI Discussion

In this work, we developed a pipeline to perform realistic anisotropy searches on GWBs produced by synthetic SMBHB populations. We then used this framework to study the prospects of anisotropy detection by current and future PTAs, under the assumption that the observed GWB is produced by a population of SMBHB. The main results of our analysis are the following:

1. Based on our analysis, the lack of evidence for GWB anisotropies in a PTA with the noise characteristics of the NANOGrav 15-year data set is not in tension with assuming SMBHBs as a source of the observed background. Indeed, we find that for a SMBHB-generated GWB there was only a $\sim 7\%$ probability of detecting anisotropies in such a data set.

2. We estimate that the probability of detecting GWB anisotropies in a SMBHB-generated background will increase to $\sim 16\%$ for the upcoming IPTA DR3 data set. This indicates that, as we keep improving the sensitivity of our observations, a null anisotropy detection will start to put pressure on an SMBHB interpretation of the background.

3. We find that the anisotropy detection probability strongly depends on the properties of the SMBHB population. For example, while the average detection probability is around $\sim 7\%$ for a PTA with the noise properties of NANOGrav 15-year data set, probabilities can reach up to $\sim 45\%$ for SMBHB populations where the phenomenological hardening processes speed up at smaller separations (i.e., when $\nu_{\rm inner}\sim 0)$. This indicate that an anisotropy detection (or lack of thereof) can convey information about the SMBHB population properties, and break the degeneracy between SMBHB populations with similar GWB spectral properties.

4. Among the considered reconstruction methods and detection statistics, the combination that performs the best uses a numerical map reconstruction method and an SNR-based detection statistic.

5. We find that only $\sim 40\%$ of the GWB skies leading to detection of anisotropies contain a binary with a CW signal above current CW upper limits. This emphasizes the complementarity between CW and anisotropy searches in the effort to detect the first direct evidence for SMBHB being the source of the GWB.