Forecasting the sensitivity of Pulsar Timing Arrays to gravitational wave backgrounds

Abstract

Pulsar Timing Array (PTA) observations hinted towards the existence of a stochastic gravitational wave background (SGWB) in the nHz frequency band. Still, the nature of the SGWB signal cannot be confidently inferred from current data, and the leading explanation invokes mergers of supermassive black holes. If confirmed, such discovery would not only represent a turning point in our understanding of astrophysics, but it may severely limit the capability of searching for additional cosmological sources in the nHz frequency range. In this work, we build a simple framework to forecast the sensitivity of future PTA configurations and assess the parameter estimation of SGWB, which could consist of several contributions. We release the python code fastPTA implementing this framework and ready to use.

I Introduction

Pulsar timing array (PTAs) experiments are unique probes to study the gravitational-wave (GW) spectrum in the nano-Hertz (nHz) frequency range. Recently, the results reported by the NANOGrav Agazie et al. (2023a, b), EPTA (in combination with InPTA) Antoniadis et al. (2023a, b, c), PPTA Reardon et al. (2023a); Zic et al. (2023); Reardon et al. (2023b) and CPTA Xu et al. (2023) collaborations hinted towards the existence of a common spectrum of a stochastic nature, with around $(2\div 4)\sigma$ significance for a Hellings-Down (HD) angular correlation, consistent with the quadrupolar nature of GWs in general relativity Hellings and Downs (1983).

All collaborations reported the signal could be explained by a stochastic GW background (SGWB) with a preference for blue-tilted spectrum of energy density $\Omega_{\rm GW}(f)\propto f^{n_{T}}$, with $n_{T}\simeq 2$, and large uncertainties varying among different collaborations. The astrophysical explanation of this signal would rely on the incoherent superposition of GW signals from a population of inspiralling supermassive black hole (SMBH) binaries on circular orbits, whose frequency spectrum is characterised by a scaling law $\Omega_{\rm GW}\propto f^{2/3}$ Phinney (2001). Nevertheless, orbital eccentricity and stellar environment could affect the binary evolution which is coupled to the statistical fluctuation in the distribution of SMBHBs may play an important role and can lead to modified spectra (see, e.g., Sesana et al. (2008); Kocsis and Sesana (2011); Kelley et al. (2017); Perrodin and Sesana (2018); Ellis et al. (2023); Agazie et al. (2023c); Afzal et al. (2023); Ghoshal and Strumia (2024); Bonetti et al. (2018)). On the other hand, with current data, it is still impossible to rule out a cosmological explanation for the observed signal. The most studied scenarios capable of producing a SGWB in the nHz frequency range would be: first-order phase transitions Arzoumanian et al. (2021); Xue et al. (2021); Nakai et al. (2021); Di Bari et al. (2021); Sakharov et al. (2021); Li et al. (2021); Ashoorioon et al. (2022); Benetti et al. (2022); Barir et al. (2023); Hindmarsh and Kume (2023); Gouttenoire and Volansky (2023); Baldes et al. (2023); Li and Xie (2023); cosmic strings and domain walls Ellis and Lewicki (2021); Datta et al. (2021); Samanta and Datta (2021); Buchmuller et al. (2020); Blasi et al. (2021); Ramazanov et al. (2022); Babichev et al. (2022); Gorghetto et al. (2021); Buchmuller et al. (2021); Blanco-Pillado et al. (2021); Ferreira et al. (2023); An and Yang (2023); Qiu and Yu (2023); Zeng et al. (2023); King et al. (2024); Babichev et al. (2023); Kitajima et al. (2024); Barman et al. (2023); scalar-induced GWs generated from primordial fluctuations Vaskonen and Veermäe (2021); De Luca et al. (2021); Bhaumik and Jain (2021); Inomata et al. (2021); Kohri and Terada (2021); Domènech and Pi (2022); Namba and Suzuki (2020); Sugiyama et al. (2021); Zhou et al. (2020); Lin et al. (2023); Rezazadeh et al. (2022); Kawasaki and Nakatsuka (2021); Ahmed et al. (2022); Yi and Fei (2023); Yi (2023); Dandoy et al. (2023); Zhao et al. (2023); Ferrante et al. (2023); Cai et al. (2023); Franciolini et al. (2023); Balaji et al. (2023); Liu et al. (2024); tensor perturbations generated during inflation (see, e.g., Vagnozzi (2023)) and many others (see also Franciolini et al. (2024); Madge et al. (2023); Figueroa et al. (2023); Garcia-Bellido et al. (2024); Murai and Yin (2023); Konoplya and Zhidenko (2023); Smarra et al. (2023)).

Future observations will narrow down the range of models compatible with the data. For example, the forecasted sensitivity of future PTA experiments may improve by various orders of magnitude compared to the ones achieved by current experiments. This will allow us to better understand the nature of the dominant source of GWB in the nHz while constraining signals from new physics appearing with amplitude below current sensitivity. Crucially, the presence of the currently observed foreground will greatly limit the constraining power of future observations on subdominant SGWB. In other words, regardless of the signal origin, currently hinted SGWB will constitute an unavoidable source of foreground noise when searching for subdominant contributions.

In this work, we aim to build a simplified framework capable of reliably estimating future PTA sensitivity without the need for expensive data simulation and analyses. We will remain agnostic on the nature of the signal, which we assume to be an isotropic, stationary SGWB characterized by the maximum likelihood power-law model resulting from the analysis of the recent PTA datasets Antoniadis et al. (2023b), and we will model noises based on the currently observed pulsars. We then forecast future PTA uncertainties and describe how they will shrink both with a longer observation time and a larger number of pulsars. Finally, we will forecast the sensitivity to subdominant contribution to the SGWB, considering simple models as a proof of concept. Detailed analysis considering motivated cosmological SGWB, coming from new physics would be presented elsewhere.

The paper is organized as follows. In Sec. II, we discuss how we build the PTA noise sensitivity, following Ref. Hazboun et al. (2019) and with some additional simplifying assumptions. In Sec. III we describe the PTA data simulations we adopt to validate our proposed methodology against state-of-the-art pipelines for the PTA data analysis. In sec. IV we report the results while presenting conclusions and outlook in Sec. V.

At the repository linked in Ref. cod , we release the code fastPTA readily usable to estimate the sensitivities and the measurement uncertainties with future PTA configuration, with longer observation times and a simulated larger set of pulsars, which also features the possibility of including additional subdominant GWB contributions.

II Building the Pulsar timing array sensitivity

In this section, we describe how, following Ref. Anholm et al. (2009); Chamberlin et al. (2015); Rosado et al. (2015); Hazboun et al. (2019), we build the PTA sensitivity curves based on a few simplifying assumptions. In the following sections, we show that this framework reproduces measurement uncertainties estimated using current data (EPTA DR2, in particular, Antoniadis et al. (2023a)) and forecasted from mock data sets generated with state-of-the-art techniques.

As a starting point, we assume the combined PTA data consists of residuals of time of arrivals (TOAs) $d_{I}(t)$, where the index $I$ runs up the total number of pulsars observed $1,\dots,N_{p}$. Moreover, we will assume the GW signal and noise to be stationary and define the Fourier domain data as

$$\tilde{d}_{i}^{k}=\int_{-T_{\rm obs}/2}^{T_{\rm obs}/2}d_{I}(t)e^{-i2\pi f_{k}t}{\rm d}t.$$ (1)

Given the finite duration of the data stream, we adopt a finite set of frequencies for the Fourier basis, $f_{k}=k/T_{\rm obs}$, starting from the inverse of the observation time $T_{\rm obs}$ and going up to the Nyquist frequency. We assume that the data contains the stochastic GW signal $\tilde{s}_{I}^{k}$ and noise $\tilde{n}_{I}^{k}$, as $\tilde{d}_{I}^{k}=\tilde{s}_{I}^{k}+\tilde{n}_{I}^{k}$, which we will further assume to be Gaussian with zero mean $\langle\tilde{s}_{I}^{k}\rangle=\langle\tilde{n}_{I}^{k}\rangle=0$.

Given the aforementioned assumptions, one can write down the full covariance matrix, including contributions from both intrinsic pulsar noise, measurement process, and GWB as

$$C_{IJ}=C_{n,IJ}+C_{h,IJ}\,.$$ (2)

In the following, we will assume the noise contribution to the covariant matrix to be dominated by the diagonal terms in the pulsar indices, i.e., to be uncorrelated among different pulsars. This allows us to approximate $C_{n,IJ}=\delta_{IJ}P_{n,I}$. Additionally, the GW contribution presents the unique Hellings and Downs (HD) correlation pattern to be discussed in Sec. II.2.

Transmission function.

The signal is extracted from timing residuals built by subtracting the expected time of arrival computed using the timing model. Fitting for the parameters entering the timing model (i.e., describing each pulsar intrinsic rotation frequency, its derivative, proper motion, etc), results in a polynomial suppression of sensitivity at low frequencies. Crucially, the drastic decrease only happens around the frequency $1/T_{\rm obs}$. This suppression of signal can effectively be described by a transmission function scaling as $1/f^{6}$ (for a quadratic spin-down model) below the frequency of $1/T_{\rm obs}$ Hazboun et al. (2019). We model this transfer function as

$${\cal T}(f)\simeq\left[1+1/\left(fT_{\text{\tiny obs}}\right)\right]^{-6}.$$ (3)

An additional loss of the sensitivity is induced around $f=1/{\rm yr}$ from fitting the sky position and the proper motion and around $f=2/{\rm yr}$ due to parallax. This generates large spikes in the sensitivity curve. Furthermore, if the pulsar is in a binary system with a period that falls within the frequency range, an additional dip in the transmission function will appear due to the fitting of the orbital parameters. As these losses typically appear at relatively high frequencies, we will neglect them in our estimates.

II.1 Noise model

Following the discussion in Ref. Hazboun et al. (2019), we isolate a few key features of the noise budget. For each pulsar, the noise contribution comes schematically from

$$\text{Noise budget}=\text{WN + RN + DM + SV},$$ (4)

defined as follows:

White noise (WN).

The time of arrival of the signal shows an overall white noise due to the finite signal-to-noise ratio (SNR) of the match filtering process adopted to extract them. It appears as a flat contribution to the power spectra of timing residuals. We characterize the white noise contribution for each pulsar as $P_{I}^{\rm WN}$, where $I$ is the index associated with each pulsar. This is typically controlled by the parameters EQUAD, ECORR, and EFAC in PTA analyses Antoniadis et al. (2023b); Agazie et al. (2023b). This can be described by a time delay power spectrum

$$P_{I}^{\rm WN}=2\sigma^{2}\Delta t,$$ (5)

where $\Delta t$ is the inverse of the observing cadence and $\sigma$ is the rms timing uncertainty.

Red noise (RN).

The stochasticity in the pulsar’s rotation causes an “achromatic" red noise, which we model with a power-law power spectral density (PSD). This noise is independent of the frequency of the radio observation (hence “achromatic”) and often dominates at low frequencies:

$$P_{I}^{\rm RN}=A^{\rm RN}_{I}\left(\frac{f}{f_{r}}\right)^{\gamma^{\rm RN}_{I}};$$ (6)

where $f_{r}$ is a reference frequency often assumed to be 1/yr, and the amplitude $A^{\rm RN}_{I}$ and the spectral index $\gamma^{\rm RN}_{I}$ is intrinsic to each pulsar $I$.

Chromatic noise components (DM, SV).

Other important noise sources are temporal variations in dispersion measure (DM) and scattering variations (SV). Those components are chromatic and add time delays to the ToAs as $\propto\nu^{-2}$ (DM) and $\propto\nu^{-4}$ (SV), where $\nu$ is the observed radio frequency, and both are caused by the time-varying electron column density in the interstellar medium along the line of sight (see Chalumeau et al. (2021) and references therein). The PSD for both noises is usually described by power-law similar to Eq. (6) but with added chromaticity.

All in all, we can approximate the noise model as a white noise contribution with the addition of a low-frequency red power law (see e.g., Ref. Caballero et al. (2016)) of the form

$$P_{n,I}=P_{I}^{\rm WN}+A^{\rm RN}_{I}\left(\frac{f}{f_{r}}\right)^{\gamma_{I}^{\rm RN}}+P_{I}^{\rm DM,SV},$$ (7)

where we arbitrarily fix the reference frequency at $f_{r}=f_{\rm yr}$ and $\gamma_{I}^{\rm RN}<0$. $P_{I}^{\rm DM,\,SV}$ in the previous equation indicate the possible presence of DM and SV terms, modeled as power-laws analogously to eq. (6), only present for some of the pulsars. In this work, we have neglected the chromaticity and treated DM and SV on an equal footing as achromatic noise. The main reason is that there is still a significant coupling between chromatic and achromatic noise components in the currently observed data. We will re-discuss this point in the concluding section.

As we will see, the simplified modeling summarised above provides a sufficiently accurate description of current sensitivity, allowing the reproduction of the uncertainties obtained in current analyses. After validating this construction, we will use this framework to forecast future detector performance. In particular, we will investigate the limit in which the currently observed SGWB dominates over the experimental noise, in what is referred to as the signal-dominated regime. This further motivates the simplifying assumptions we made to describe the noise model.

II.2 Stochastic GW background model

The spin-2 perturbations of the metric can be expressed in terms of plane GWs with frequency $f$, polarizations $\{+,\times\}$, and propagation directions $\hat{k}$ as

$$h_{ab}(t,\vec{x})=\int{\rm d}^{2}\Omega_{\hat{k}}\>\int_{-\infty}^{\infty}{\rm d}f\>\left[\tilde{h}_{+}(f,\hat{k})e^{+}_{ab}(\hat{k})\right.\\ \left.+\tilde{h}_{\times}(f,\hat{k})e^{\times}_{ab}(\hat{k})\right]e^{i2\pi f(t-\hat{k}\cdot\vec{x}/c)}\,,$$ (8)

where we introduced the polarization tensors $e^{+,\times}_{ab}(\hat{k})$, with the two polarizations denoted $P=\{+,\times\}$ in the following. We assume that the SGWB is stationary, unpolarized, and isotropic, which means

$$\langle\tilde{h}_{P}(f,\hat{k})\tilde{h}_{P^{\prime}}^{*}(f^{\prime},\hat{k}^{\prime})\rangle=\frac{1}{16\pi}S_{h}(f)\delta(f-f^{\prime})\delta_{PP^{\prime}}\delta^{2}(\hat{k},\hat{k}^{\prime})\,,$$ (9)

where $S_{h}(f)$ is the (one-sided) strain power spectral density of the GWB. We define the SGWB energy density (per logarithmic frequency interval) as (see e.g., Caprini and Figueroa (2018))

$$\Omega_{\rm GW}h^{2}=\frac{h^{2}}{\rho_{c}}\frac{d\rho_{\rm GW}}{d\log f}\equiv\frac{2\pi^{2}f^{3}}{3H_{0}^{2}/h^{2}}S_{h},$$ (10)

where $\rho_{c}/h^{2}=3(H_{0}/h)^{2}/(8\pi G)$ is the Universe critical energy density and $H_{0}/h=1/(9.78{\rm Gyr})$ is the Hubble parameter today. In the last equality, we have associated $\Omega_{\rm GW}h^{2}$ to the strain power spectral density $S_{h}$.

The timing residual response of the signal coming from a pulsar $I$ to the SGWB can be expressed in Fourier space as Hazboun et al. (2019)

$$\tilde{h}_{I}(f)=\int{\rm d}^{2}\Omega_{\hat{k}}\>\sum_{P}R_{I}^{P}(f,\hat{k})\tilde{h}_{P}(f;\hat{k}),$$ (11)

which is integrated over all propagation directions $\hat{k}$. The response function for a pulsar located at a distance $D$ along the direction $\hat{p}_{I}$ is

$$R_{I}^{P}(f,\hat{k})\equiv\frac{\epsilon^{P}_{ab}(\hat{k})}{i4\pi f}\frac{\hat{p_{I}}^{a}\hat{p_{I}}^{b}}{1+\hat{p_{I}}\cdot\hat{k}}\left(1-e^{-i2\pi fD(1+\hat{k}\cdot\hat{p}_{I})/c}\right)\,.$$ (12)

When dealing with PTA observations, it is easy to see that the characteristic frequency $f_{*}=(2\pi D/c)^{-1}$ associated to a distance $D$ (which is of order kpc) turns out to be $f_{*}\simeq 2\times 10^{-12}$Hz. On the other hand, the minimum frequency accessible in our analysis is limited by the observation time $f\sim 1/T_{\rm obs}\sim$nHz. Therefore, one finds that for PTA experiments $fD\gg 1$, and the frequency-dependent term in the response function (12) is well approximated by $R^{P}(f)\propto 1/f$, while the rapidly oscillating piece is negligible when averaged over the directions $\hat{k}$ in (11).

It follows that the SGWB signal covariance matrix is

$$C_{h,IJ}\equiv\langle\tilde{h}_{I}\tilde{h}_{J}^{*T}\rangle=R_{IJ}S_{h}(f),$$ (13)

where we introduced Hazboun et al. (2019)

$$R_{IJ}=\chi_{IJ}\cdot{\cal R}(f)\left[{\cal T}_{I}(f){\cal T}_{J}(f){T_{IJ}}/{T_{\text{obs}}}\right]^{1/2},$$ (14)

and we conveniently defined the frequency dependent, sky-averaged, quadratic response function ${\cal R}(f)\equiv 1/12\pi^{2}f^{2}$. The first geometrical factor is the well-known HD correlation pattern as a function only of the angular separation between pulsars as required by symmetries (see e.g. Kehagias and Riotto (2024)), discussed below. In contrast, the time-dependent factor is introduced to account for each pulsar transmission function as well as the individual observation time of each pulsar $T_{I}$. This is because we can only include off-diagonal components correlating signals between different pulsars for an effective overlapping time defined as $T_{IJ}={\rm min}[T_{I},T_{J}]$.

II.2.1 The Hellings-Downs correlation

We introduced the HD function Hellings and Downs (1983) for a pair of pulsars separated by an angle $\zeta_{IJ}\equiv\arccos(\hat{p}_{I}\cdot\hat{p}_{J})$, which takes the analytic expression

$$\chi_{IJ}=\frac{1}{2}+\frac{3}{2}\xi_{IJ}\left[\ln\xi_{IJ}-\frac{1}{6}\right]+\frac{1}{2}\,\delta_{IJ}\,,$$ (15)

where $\xi_{IJ}\equiv\left({1-\cos(\zeta_{IJ})}\right)/{2}$.

To estimate the sensitivity to the HD angular correlation function, it is convenient to decompose it on a suitable basis. Following Ref. Antoniadis et al. (2023a), we adopt the following templates:

• Binned HD function: we model the correlation function as a step-wise constant over $n$ bins in the angular variable $\xi_{IJ}$. This takes the form

  $$\chi_{IJ}=\sum_{i=0}^{n-1}b_{i}\Theta(\zeta_{IJ}-\zeta^{i})\Theta(\zeta^{i+1}-\zeta_{IJ}),$$ (16)

  where the bins in the angular variables are equally spaced $\zeta^{i}=i\pi/n$. We also introduced the Heaviside theta function $\Theta$. We impose Gaussian priors on the coefficients $b_{i}$ with width $\sigma_{b_{i}}=1$.

• HD expansion in Legendre polynomials: following Gair et al. (2014); Antoniadis et al. (2023a), we expand the HD function as

  $$\chi_{IJ}=\sum_{\ell=0}^{n}a_{\ell}P_{\ell}(\cos\zeta_{IJ}).$$ (17)

  Using the standard normalization of the Legendre polynomials, the coefficients are found by Gair et al. (2014); Romano and Allen (2023)

  $$a_{\ell} =\frac{2\ell+1}{2}\int_{-1}^{1}\chi_{\text{\tiny HD}}(x)P_{\ell}(x)\textrm{d}x$$

  $$=\frac{3}{2}\frac{(2\ell+1)}{(\ell+2)(\ell+1)\ell(\ell-1)}.$$ (18)

  for $\ell\geq 2$, and $a_{0}=a_{1}=0$. We will estimate the sensitivity to the HD curve by constraining the expansion parameters $a_{\ell}$. We impose priors on $a_{\ell}$ to be Gaussian with width $\sigma_{a_{\ell}}=1$.

II.2.2 Effective sensitivity

We compute the SNR by

$${\rm SNR}^{2}\equiv\sum_{f_{k}}\left[C^{-1}_{IJ}C^{-1}_{KL}R_{JK}R_{LI}\right]S_{h}^{2},$$ (19)

where, here and throughout the work, we implicitly assume Einstein’s summation convention on the pulsar indices in capital letters. It is instructive to define an effective sensitivity, as Hazboun et al. (2019)

$$S_{\rm eff}(f)=\left(C^{-1}_{IJ}C^{-1}_{KL}R_{JK}R_{LI}\right)^{-1/2},$$ (20)

in such a way that the SNR computation takes the familiar form ${\rm SNR}^{2}=T_{\text{obs}}\sum_{k}(S_{h}/S_{\text{eff}})^{2}$.

In the weak signal limit (which is in practice the limit in which ${\cal R}S_{h}\ll P_{n}$), one can expand the correlation matrix above to find

$$S_{\rm eff}(f)=\left(\sum_{I\neq J}\frac{T_{IJ}}{T_{\text{obs}}}\frac{\mathcal{T}_{I}(f)\mathcal{T}_{J}(f)\chi^{2}_{IJ}}{P_{n,I}(f)P_{n,J}(f)/\mathcal{R}^{2}(f)}\right)^{-1/2}\,,$$ (21)

which includes contributions from the Hellings and Downs factors $\chi_{IJ}$ and the individual pulsar strain-noise power spectral densities. This corresponds to the result of Ref. Hazboun et al. (2019) and shown in Refs. Agazie et al. (2023d, e) to report current PTA sensitivities. As we will see, when considering future PTA datasets (including longer observation times and a larger number of pulsars), it is important to adopt the full expression for the effective sensitivity to go beyond the weak signal approximation. We report examples of future effective sensitivity, in the presence of the currently observed GW background in Fig. 2 in the next section.

II.3 Likelihood function and parameter estimation

The log-likelihood for Gaussian and zero mean data $\tilde{d}^{k}$, with $k$ running over frequencies $f_{k}$, described only by their variance, dubbed $C_{IJ}(f_{k},{\theta})$, can be written as Contaldi et al. (2020); Bond et al. (1998)

$$-\ln\mathcal{L}(\tilde{d}|{\theta})\propto\sum_{k,IJ}\ln\left[C_{IJ}(f_{k},{\theta})\right]+\tilde{d}_{I}^{k}C_{IJ}^{-1}(f_{k},{\theta})\tilde{d}_{J}^{k*},$$ (22)

which is also known as Whittle likelihood. This likelihood function for ${\theta}$ corresponds to the probability of the data $\tilde{d}$ given ${\theta}$, the signal model parameters.

Fisher information matrix estimates.

Here we briefly summarise the main ingredients of the Fisher Information Matrix (FIM) formalism, which is often used in the limit of large SNR to assess parameter uncertainties.

The FIM $F_{\alpha\beta}$ is defined as

$$F_{\alpha\beta}\equiv-\left.\frac{\partial^{2}\log\mathcal{L}}{\partial{\theta}^{\alpha}\partial{\theta}^{\beta}}\right|_{{\theta}={\theta}_{0}}=\sum_{k}\textrm{Tr}\left[C^{-1}\frac{\partial C}{\partial{\theta}^{\alpha}}C^{-1}\frac{\partial C}{\partial{\theta}^{\beta}}\right]_{{\theta}={\theta}_{0}}\;,$$ (23)

where, ${\theta}_{0}$ represents the maximum likelihood estimator of model parameters, determined by imposing

$$\left.\frac{\partial\log\mathcal{L}}{\partial{\theta}^{\alpha}}\right|_{{\theta}={\theta}_{0}}\propto\sum_{k}\frac{\partial C}{\partial{\theta}^{\alpha}}\left[C^{-1}-C^{-1}\;\tilde{d}^{k}\tilde{d}^{k*}\;C^{-1}\right]=0\;,$$ (24)

with $C(f_{k},{\theta}_{0})=\tilde{d}^{k}\tilde{d}^{k*}$. In practice, the discrete sum over finite frequencies can be replaced with a continuous integral over the frequency range as

$$F_{\alpha\beta}\equiv\sum_{f_{k}}\frac{\partial\log C}{\partial{\theta}^{\alpha}}\frac{\partial\log C}{\partial{\theta}^{\beta}}\,.$$ (25)

Keeping all indices fully expressed, one obtains

$$F_{\alpha\beta}\equiv\sum_{f_{k}}C^{-1}_{IJ}C^{-1}_{KL}\frac{\partial(R_{JK}S_{h})}{\partial{\theta}^{\alpha}}\frac{\partial(R_{LI}S_{h})}{\partial{\theta}^{\beta}}.$$ (26)

Finally, the covariance matrix $C_{\alpha\beta}$, is obtained by inverting the FIM, from which one can estimate uncertainties as $\sigma_{\alpha}\equiv\sqrt{F_{\alpha\alpha}^{-1}}$.

MCMC analyses.

Given that the FIM is an unreliable estimator in the low SNR limit, where the validity of the Gaussian approximation is expected to degrade, we will test the FIM results with full-fledged Bayesian parameter estimation. For this purpose, we simulate frequency domain data $\tilde{d}^{k}\equiv\tilde{d}(f_{k})$, where $k$ indexes the frequencies within the detector sensitivity range. We generate Gaussian realizations for the signal and all noise components, with zero mean and variances defined by their respective power spectral densities. The parameter estimation is then performed by running a Monte Carlo Markov Chain (MCMC) using emcee Foreman-Mackey et al. (2013) with uninformative priors on signal parameters.

III Simulating future pulsar timing array networks

To test our method, we compare to currently available datasets, while we also propose to simulate future generation datasets. In the upcoming years, the Square Kilometer Array (SKA) will be a central protagonist among PTA collaborations Janssen et al. (2015). At the moment, PTA datasets are limited by radiometer noise. The SKA will provide high-precision pulsar timing measurement with uncertainties below $\sim$100ns Janssen et al. (2015), making it roughly 10 times better than current generation telescopes Lazio (2013). With a reduction of the white noise at relatively high frequencies, it also comes with a better determination of the RN, relevant at low frequencies. For our purposes, we will simulate two mock datasets:

• An EPTA-like with 50 pulsars and 20 yr observations. We will denote this as EPTA20 in the rest of the paper.

• A SKA-like with 50 pulsars and 10 yr observations. We will denote this configuration as SKA10.

III.1 Fake PTA data and time-domain Bayesian analysis

III.1.1 Time domain representation

As mentioned in Sec. II, we are dealing with the timing residuals obtained by fitting and subtracting the timing model from observed pulses’ ToAs. The residuals from all pulsars are concatenated in a single array forming PTA data set $\vec{d}=[d_{1}(t),d_{2}(t),...,d_{N_{\rm pulsars}}(t)]$. We expect that the residual errors in fitting the timing model parameters are small and could be approximated by a linear model

$$d_{I}^{\prime}\rightarrow d_{I}-\sum_{k}\alpha_{k}\vec{M}_{Ik},$$ (27)

where $M_{Ik}$ is a design matrix (representing the derivative of the timing model w.r.t timing model parameters), and $\alpha_{k}$ are residual timing model errors, which we consider as random variables. We have already used the design matrix in the frequency domain while building the transmission function.

For the description of the noise components, we use Gaussian process representation (see van Haasteren and Vallisneri (2014) for a detailed description) with a discrete Fourier basis $f_{i}=i/T_{I}$ with $i=[1,...,N_{f}]$ as the eigenfunctions:

$$n_{I}(t)=\sum_{i}X_{i,I}\cos(2\pi f_{i}t)+Y_{i,I}\sin(2\pi f_{i}t),$$ (28)

where $X_{i,I},Y_{i,I}\sim\mathcal{N}(0,S(f_{i})\Delta f_{i})$, $S(f)$ the one-sided PSD of the noise and $\Delta f_{i}=f_{i+1}-f_{i}$. The weights $X_{i,I},Y_{i,I}$ for the pulsar intrinsic noise components are uncorrelated among the pulsars. For GWB signals, they are spatially correlated between pulsars $I$ and $J$ and are distributed as a $N_{\rm pulsars}$ dimensional multivariate normal distribution $\vec{X}_{i},\vec{Y}_{i}\sim\mathcal{N}(0,\chi_{IJ}S(f_{i})\Delta f_{i})$ where $\chi_{IJ}$ are the introduced above HD correlation coefficients. In general, the noise covariance matrix is given as

$$C_{I}^{n} = \langle n_{I}(t)n_{I}(t)\rangle$$

$$= \begin{bmatrix}\cos(2\pi f_{1}\vec{t}_{I})\\ \sin(2\pi f_{1}\vec{t}_{I})\\ \vdots\\ \cos(2\pi f_{N_{f}}\vec{t}_{I})\\ \sin(2\pi f_{N_{f}}\vec{t}_{I})\end{bmatrix}\cdot\begin{bmatrix}\langle X_{1,I}^{2}\rangle\\ \langle Y_{1,I}^{2}\rangle\\ \vdots\\ \langle X_{N_{f},I}^{2}\rangle\\ \langle Y_{N_{f},I}^{2}\rangle\end{bmatrix}\cdot\begin{bmatrix}\cos(2\pi f_{1}\vec{t}_{I})\\ \sin(2\pi f_{1}\vec{t}_{I})\\ \vdots\\ \cos(2\pi f_{N_{f}}\vec{t}_{I})\\ \sin(2\pi f_{N_{f}}\vec{t}_{I})\end{bmatrix},$$

where $\langle X_{1,I}^{2}\rangle$ being the variance of the $X_{i,I},Y_{i,I}$ as defined by their Gaussian distribution.

III.1.2 Generating fake data

We simulate the data in the time domain. We start with choosing the pulsar array and epochs of observation, which we populate with the white noise with dispersion $\sigma$. We simulated the red (spin) noise of each pulsar using Gaussian Process according to (28) where the weights $\vec{X}_{i,I}$ and $\vec{Y}_{i,I}$ are drawn from their respective normal probability distributions with variance defined by PSD, $S(f)$. The simulated timing model includes only five components with the following design matrix:

• Offset : $\vec{M}_{0}\propto\vec{1}$,

• Spin rate : $\vec{M}_{1}\propto\vec{t}_{I}$,

• Quadratic spin-down rate : $\vec{M}_{2}\propto\vec{t}_{I}^{2}$,

• RA : $\vec{M}_{3}\propto\sin(2\pi f_{yr}\vec{t}_{I})$,

• DEC : $\vec{M}_{4}\propto\cos(2\pi f_{yr}\vec{t}_{I})$.

The overall offset, spin rate, and quadratic spindown rate correspond to the inaccuracy in the pulsar rotation. The offset in the pulsar sky position is modeled by RA and DEC terms, giving annual modulation due to the Earth’s orbital motion. During the Bayesian PTA analysis, we marginalize over the coefficients $\alpha_{k}$ of the linear timing model, assuming non-informative improper prior. This marginalization is responsible for the transmission function described in (3) Hazboun et al. (2019).

III.2 Mock datasets

We simulated two PTA mock datasets:

• EPTA20 : built from the 25 pulsars of the EPTA DR2new 10.4 yr dataset with the same noise properties (as in Antoniadis et al. (2023d)) but with a doubled observation time, to which we append 25 additional pulsars with the same noise properties and observation time as EPTA DR2new but randomized positions drawn uniformly in the sky. The noises are re-injected according to their maximum likelihood values in Antoniadis et al. (2023d).

• SKA10 : SKA-like PTA with 50 pulsars and 10 yr observations. Every pulsar is simulated with a 2-week cadence of observation, $\sigma=$100ns timing uncertainty, and random position uniformly drawn in the sky. Time-correlated noises with power-law spectra are included (red noise and dispersion measure noise) with log-amplitude $\log_{10}A$ and spectral index $\gamma$ drawn uniformly between $\log_{10}A=[-17,-13]$ and $\gamma=[1,5]$.

These properties are summarised in Tab. 1.

	EPTA20	SKA10
$N_{\rm pulsars}$	25 + 25	50
$\nu$ [GHz]	Antoniadis et al. (2023d)	1.4 - 3
$T_{\rm obs}$ [yr]	20.8 + 10.4	10
$\sigma$ [s]	$\sim 10^{-6}$	$10^{-7}$
$\Delta t$ [days]	$\sim 3$	14
TM	TM5 + Antoniadis et al. (2023d)	TM5
$\log_{10}A_{\rm RN}$	Antoniadis et al. (2023d)	[-17, -13]
$\gamma_{\rm RN}$	Antoniadis et al. (2023d)	[1, 5]
$\log_{10}A_{\rm DM}$	Antoniadis et al. (2023d)	[-17, -13]
$\gamma_{\rm DM}$	Antoniadis et al. (2023d)	[1, 5]
$\log_{10}A_{\rm SV}$	Antoniadis et al. (2023d)	-
$\gamma_{\rm SV}$	Antoniadis et al. (2023d)	-

Table 1: Summary of mock datasets properties. $N_{\rm pulsars}$ is the number of pulsars, $\nu$ are the observing frequency sub-bands, $T_{\rm obs}$ is the time of observation, $\sigma$ is the level of white noise, $\Delta t$ is the cadence, TM is the timing model where TM₅ is the five-component model described in III.1.2, $\log_{10}A$ and $\gamma$ give the range of the uniform distribution used to draw the red noise parameters. The entries marked by Ref. Antoniadis et al. (2023d) correspond to the choice of parameters consistent with ERPTA DR2new analysis.

Figure 1: Comparison between the measurement uncertainties estimated using FIM (blue) and time-domain Bayesian parameter estimation (red, see Sec. III.1 for more details). For visualization purposes, the FIM result is shown by building samples drawn from a multivariate Gaussian distribution with covariance $FIM^{-1}$. For both parameters, we remove the mean observed values to compare with the FIM which does not provide information on it. In the two-dimensional panels, we indicate 68%, 95%, and 99% C.I. Left panel: current posterior obtained with EPTADR2new datasets (EPTA10). Center panel: forecasted EPTA 20 yr dataset. Right panel: forecasted SKA 10 yr dataset.

In both datasets, we inject a stochastic GW signal with HD correlations and a power-law spectrum with $\log_{10}A=-14$ and $\gamma=3$ to reproduce the signal observed in Antoniadis et al. (2023a); Agazie et al. (2023a), potentially originating from the population of SMBHB in the Universe. This signal will act as foreground, and we will also inject additional stochastic GW signals into the datasets to test our ability to detect them in the presence of a foreground (see Sec. IV.4).

III.3 Time-domain Bayesian analysis

The timing model marginalized Gaussian likelihood used in time-domain PTA analysis van Haasteren and Vallisneri (2014) can be expressed as (analogously to the previously defined likelihood (22))

$$-\ln\mathcal{L}(d|{\theta})\propto\textrm{det}\{C_{IJ}\}+\sum_{IJ}d_{I}(t)C_{IJ}^{-1}({\theta})d_{J}^{\intercal}(t),$$ (29)

where $C_{IJ}=C_{I}\delta_{IJ}+C_{h,IJ}$ has a block structure of size $N_{\rm pulsars}\times N_{\rm pulsars}$ in which the GW signal $C_{h,IJ}=\chi_{IJ}C_{h}$ ($I\neq J$) appears off diagonal. The individual pulsar’ noise components form the diagonal $I=J$ where $C_{I}=C_{I}^{TM}+C_{I}^{RN}+C_{I}^{DM}+C_{I}^{SV}$ (corresponding to the timing model, the red noise, the dispersion measure noise, and the chromatic noise respectively).

We perform Bayesian analysis using the software ENTERPRISE Ellis et al. (2020); Taylor et al. (2021) and sampler PTMCMC Ellis and van Haasteren (2017). ENTERPRISE provides computation of the likelihood and prior based on the chosen model. All noise components are modelled as Gaussian Processes with incomplete Fourier basis as eigenfunctions van Haasteren and Vallisneri (2014). We are using PTMCMC without a parallel tempering option, which is sufficient for our purpose. The resulting posterior distributions for both EPTA20 and SKA10 are shown in Fig. 1 (central and right panels, red color).

IV Results

Throughout this section, we will assume that the SGWB observed by the various collaborations is described by a power-law (PL)

$$h^{2}\Omega_{\rm GW}^{\rm PL}(f)=A_{\rm GWB}\left(\frac{f}{f_{\rm yr}}\right)^{n_{T}}$$ (30)

with $A_{\rm GWB}\simeq 6.3\times 10^{-8}$ and $n_{T}\simeq 2$. For definiteness, these parameters are fixed to the maximum likelihood values obtained by the EPTA Antoniadis et al. (2023a), which is in good agreement with the other observations Agazie et al. (2023e). To make contact with other quantities used in the literature, this corresponds to a strain amplitude ($h_{c}\equiv\sqrt{fS_{h}}$) at $f=f_{\text{yr}}$ of $A_{h_{c}}\simeq 10^{-14}$ and pulsar timing residuals spectral index $\gamma=5-n_{T}=3$. When performing parameter estimation, it is convenient to define the log of the amplitude $\alpha_{*}\equiv\log_{10}(A_{\rm GWB})$.

Assuming this SGWB, we obtain SNR = 4.5 with the current EPTA (DRNew) pulsar configuration, which is compatible with the observations Antoniadis et al. (2023a). On the other hand, we obtain ${\rm SNR}=184$ for EPTA20 and ${\rm SNR}=292$ for SKA10.

To validate our framework based on the approximated likelihood introduced in section II.3, we compare our results to the ones obtained with state-of-the-art techniques described in III.3. In Fig. 1, we compare the uncertainties obtained with our simplified FIM approach to actual PTA parameter estimation. We remove the mean from the posterior distribution of both amplitude and tilt showing the deviations from the mean values. In the left panel, we compare results for the current EPTA DR2new configuration, assuming $T_{\text{obs}}=10.33$yrs of data with the 25 pulsars adopted by EPTA in their analysis Antoniadis et al. (2023b). In the center and right panels, we compare our forecasted sensitivity to signal parameters with the simulations for EPTA20 and SKA10, with $T_{\text{obs}}=20$yrs and $T_{\text{obs}}=10$yrs, respectively and $N_{\rm pulsars}=50$, built as described in Sec. III. In all cases, we find very good agreement between both results, showing the validity of our assumptions in Sec. II for our purposes.

IV.1 Future effective sensitivity

Using the metrics defined in sec. II.2.2, we can forecast the effective sensitivity achieved by current and future experiments. In particular, we plot $S_{\rm eff}(f)$ Fig. 2 for different detector designs, i.e., EPTA10 and SKA10, translated in the GW energy density using Eq. (10). Notice we plot the effective sensitivity as a function of frequencies also below $f\sim 1/(10{\rm yrs})$ for presentation purposes. In the same plot, we show the maximum likelihood power-law SGWB for comparison. These results assume the HD function to be fixed and described by Eq. (15) (i.e., we do not infer it from the data).

In both cases, the dashed lines report the putative sensitivity obtained without injecting an SGWB. The effective sensitivity reaches a minimum of around $f\sim 3$nHz, with a steep high-frequency slope mostly induced by the large WN components. Indeed, in the high-frequency part of the sensitivity curve, we observe $\Omega_{\rm GW}(f)\sim f^{5}$, as expected from a WN-dominated regime, including frequency-dependent factors coming from the definition of energy density and the response function. Instead, at low frequencies, the sensitivity is degraded with a tilt steeper than $\propto 1/f$, expected from the behavior of transmission functions alone, showing the relevance of the red-noise components included in the pulsar noise budget. SKA10 sensitivity greatly surpasses EPTA10 due to the much smaller expected WN (with an improvement that is more than an order of magnitude at high frequencies).

The solid lines report the effective sensitivity obtained including the SGWB background as in Eq. (30). While including the SGWB does not significantly affect the sensitivity of EPTA10, showing that one still falls in the weak signal regime, the SKA10 effective sensitivity strongly degrades, driven by the large SGWB foreground compared to the pulsar noises.

In practice, the take-home message is that searches for subdominant contributions to the SGWB below the solid lines would be inconclusive, at odds with the conclusion drawn considering the effective sensitivity with neglected foregrounds. This is one of the main results of this paper. We will show examples of multi-SGWB analyses in the following sections, corroborating the qualitative conclusions shown here.

Figure 2: Current sensitivity of the EPTADR2new dataset, denoted as EPTA10, (blue lines), and forecasted future sensitivity achieved by the SKA 10 yr (yellow lines). Both scenarios are reported assuming the presence (absence) of an SGWB in a solid (dashed) line. The black solid line corresponds to the maximum likelihood power-law signal (30).

Figure 3: Left panel: scaling of the SNR as a function of the observation time $T_{\text{\tiny obs}}$ assuming a PTA network composed of 30 pulsars. The vertical error bars indicate the standard deviation of the result obtained by varying over a random realization of the PTA array (both for pulsar positions and noises). The black dashed lines indicate the two limiting scaling regimes expected for the SNR. The low-$T_{\rm obs}$ regime is dictated by our choice of the spectrum with $n_{T}=2$, while the high-$T_{\rm obs}$ regime is universal. Right panel: scaling of the uncertainties on the amplitude and tilt of the SGWB as a function of observation time.

Figure 4: Scaling of the SNR (left panel) and uncertainties (right panel) as a function of the number of pulsars, assuming fixed $T_{\text{\tiny obs}}=20$yrs. The error band indicates the standard deviation of both results induced by producing multiple realizations of the pulsar sample.

IV.2 Scaling with time and number of pulsars

With the proposed framework, we can forecast the evolution with time and the number of pulsars of both the SNR and measurement precision.

In Fig. 3, we report the results of different random realizations of $30$ pulsars with variable observation time $T_{\text{obs}}$. We indicate with the uncertainty bands the corresponding spread obtained from randomly selecting pulsars with different sky localization (assuming isotropic distribution) and noises. The parameters determining the noise of each pulsar are sampled from distributions built out of the currently observed ones Antoniadis et al. (2023b).

In these plots, we aim to test the scaling of the relevant quantities. For this reason, we show FIM results also for low-SNR values where FIM predictions are unreliable. It should be kept in mind that results obtained in this regime only serve the purpose of testing the scalings and should not be used to deduce the detectability of the signal. Also, in this case, these results assume the HD function to be fixed and described by Eq. (15) (i.e., we do not infer it from the data). We do not expect this assumption to modify the scaling obtained.

Figure 5: Uncertainties on the HD curve were obtained using the Fisher matrix analysis and assuming a power-law signal. Left panel: constraints on the coefficients of the binned HD function. We divide the angular direction into 7 equally spaced bins. Uncertainties are comparable with the one observed found in EPTADR2 analyses Antoniadis et al. (2023a). By construction, within the FIM formalism, the central value corresponds to the injection. Right panel: Posterior predictive distribution of the HD function obtained from the uncertainties on the coefficients of the Legendre polynomial expansion in Eq. (17).

Figure 6: Same as Fig. 5, but assuming the SKA10 future configuration.

Focusing on the evolution of SNR, we observe that at small observation times, it grows rapidly as SNR$\sim T_{\text{\tiny obs}}^{3}$. This can be explained as follows (see analytical considerations reported in Ref. Siemens et al. (2013)). Assuming WN dominates at the initially relatively high frequencies $f_{\text{min}}\sim 1/T_{\text{\tiny obs}}$, and that $\gamma>0$, the integral for the SNR is dominated by the contribution close to $f_{\text{min}}$. Longer observation times allow to probe lower frequencies, enhancing the ratio $S_{h}/S_{\text{eff}}\sim f_{\text{min}}^{-\gamma}\sim T_{\text{obs}}^{\gamma}$. Including the additional pre-factor and frequency summation, this gives SNR$\propto T_{\text{obs}}^{\gamma}$. When the effective sensitivity falls below the GW amplitude at the minimum frequency, one transitions to the intermediate regime, and the above scaling is broken. Eventually, the SNR growth with time converges to $\propto\sqrt{T_{\text{obs}}}$ as predicted in the strong signal regime, which we observe on the right side of the plot.

An analogous behavior is observed in the evolution of uncertainties on signal parameters $\alpha_{*}$ and $n_{T}$. By construction, they are expected to scale inversely with the SNR, in the high SNR regime. We observe that both feature an initially steep decline that converges towards a $1/\sqrt{T_{\rm obs}}$ in the far future. We stress that current PTA experiments are currently probing the intermediate regime between the two asymptotic scalings, seen in the plot around $T_{\rm obs}\sim 10$yrs (see also Pol et al. (2021)). Additionally, we see that the uncertainty on the tilt reaches the scaling regime more slowly because more frequency modes need to be observed to better pin down the spectral tilt compared to the overall amplitude.

In Fig. 4 we show the corresponding scaling as a function of the number of pulsars. In the left panel, we see that the SNR grows as $\sqrt{N}$ as a function of number of pulsars. This scaling originates from including auto-correlation terms in our analysis (i.e., diagonal terms in the pulsar indices of the covariant matrix). Correspondingly, the uncertainties on SGWB spectral amplitude and tilt decrease as $1/\sqrt{N_{\rm pulsars}}$. Had we only included the off-diagonal components, featuring the HD correlation, we would have obtained a scaling ${\rm SNR}\sim\left(\sum_{I\neq J}\chi_{IJ}^{2}\right)^{1/2}$ which behaves as ${\rm SNR}\sim\left[N_{\rm pulsars}(N_{\rm pulsars}-1)\right]^{1/2}\sim N_{\rm pulsars}$.

Figure 7: Uncertainties on the fifth-order expansion coefficients of the HD curve in Legendre’s polynomials $a_{\ell}$ as a function of the number of pulsars $N_{\rm pulsars}$. We adopt $N_{\rm modes}=30$. Left panel: zero noise limit. The $0^{\rm th}$ and $1^{\rm st}$ orders are zero, and the corresponding uncertainty scales as $\propto 1/N_{\rm pulsars}$ (indicated by the dashed black line). The higher-order pieces reach the plateau predicted by cosmic variance in Eq. (31). See the discussion in the main text. Right panel: Same as the left panel, but includes noise. The plateau at small values of $N_{\rm pulsars}$ is reached when information is small enough that our prior on $\sigma_{a_{\ell}}=1$ becomes dominant.

IV.3 Sensitivity to the HD correlation pattern

In this section, we show how the formalism developed in the current work can also forecast the sensitivity to the HD correlation. This is a crucial observable to identify GW signals and could be used to test modified theories of gravity (e.g., Liang et al. (2023)). For this purpose, we allow $\chi_{IJ}$ to vary, and, adopting the HD curve parameterizations described in Sec. II.2.1, we assess the ability of different PTA configurations to constrain the expansion coefficients.

First of all, we consider the case of the EPTA10 configuration. In Fig. 5 (left panel), we report the uncertainties on the coefficient $b_{i}$ of the binned parameterization (16) found using FIM. For this analysis, we considered 7 bins as well as the signal amplitude and slope as free parameters. As one can notice, the uncertainties are large both at narrow and wide angular separations, where fewer pulsar pairs are present, while reaching down to values around $\Delta b_{i}\sim 0.3$ at $68\%$ C.L. in the central bins. These results are in good agreement with the posterior distributions reported by the EPTA collaboration Antoniadis et al. (2023a), although they choose unequally spaced bins to evenly distribute the number of pulsar pairs as a function of angular separation.

In the right panel of Fig. 5, we show the posterior predictive distribution of $\chi_{IJ}$ obtained using the uncertainties estimated with FIM for the current EPTA10 configuration. In this case, the FIM analysis adopts the coefficients of the Legendre polynomial expansion (17) up to the sixth order. As expected, relatively tight bounds are obtained in the central portion of the angular separation variable (which contains a larger number of pulsar pairs), with constraints that degrade at the edges of the plot.

Having checked the validity of our framework to estimate the sensitivity to HD correlations, we forecast the corresponding uncertainties obtained with the future SKA10 dataset. In Fig. 6, we show the corresponding posterior on binned expansion parameters (left panel) as well as posterior predictive distribution reconstructed using Legendre expansion (right panel). As expected, the improved SNR ratio in the SKA10 forecast leads to significantly narrower uncertainties on the HD correlation function. In the central bins, one finds uncertainties around $\Delta b_{i}\sim 0.04$ at $68\%$ C.L.

Cosmic variance.

For a Gaussian background, the multipole moments are drawn from a Gaussian distribution with variance $a_{\ell}$ defined in Eq. (17). Even neglecting noise, and assuming perfect measurement precision, the estimate of the variance of the distribution is limited by the number of modes observed for each $\ell$. As described in Refs. Roebber and Holder (2017); Allen (2023), assuming single mode observation and no noise, the uncertainty on $a_{\ell}$ should plateau towards the value

$$\Delta a_{\ell}=\frac{a_{\ell}}{\sqrt{2\ell+1}}\;.$$ (31)

In reality, multiple frequency modes are observed by PTA, and one expects the limit on the uncertainties due to cosmic variance to scale as $\Delta a_{\ell}\sim\left(N_{\rm modes}\right)^{-1/2}$.

In Fig. 7, we show how the uncertainty on $a_{\ell}$ scales as a function of the number of pulsars. To better show the qualitative behavior, we first compute the results assuming no noise and $30$ observed frequencies (left plot) and $T_{\rm obs}=10$yr, only including terms up to $\ell=6$ in the FIM. As one can see, the uncertainty for $\ell=0,1$ scales as $1/N_{\rm pulsars}$, without showing any significant departure at high $N_{\rm pulsars}$. This is because, $a_{0,1}=0$ due to the quadruple nature of GR, and their inference does not suffer from self-noise. On the other hand, the uncertainty on high multiples converges towards a plateau, which is dictated by cosmic variance (31). The plateau is reached at an increasing number of pulsars for higher $\ell\geq 2$ due to the increasingly larger $N_{\rm pulsars}$ needed to cover the sky positions with sufficient finer angular resolution to saturate the information available.

In the right panel of Fig. 7, we show analogous results assuming the noise to be present. In this case, we simulate 30 pulsars observed for $T_{\text{obs}}=10$yrs. We see that uncertainties on $a_{\ell}$ feature a plateau at low values of $N_{\rm pulsars}$, where their posterior distribution is dominated by the prior (assumed to be Gaussian with $\sigma_{a_{\ell}}=1$). When $N_{\rm pulsars}\gtrsim 10$ (incidentally close to the currently analyzed sample of pulsars by EPTA), uncertainties start to decrease, especially at lower multipoles. Qualitatively, we observe the same behavior obtained in the left panel for the idealized case of no noise. The $\ell=0,1$ uncertainties keep the $1/N_{\rm pulsars}$ scaling, while starting from the lowest orders in $\ell$ (especially $\ell=2,3$) we notice a flattening of the curve already from $N_{\rm pulsars}\sim{\cal O}({\rm few}\cdot 10^{2})$. We do not push the computations to even larger $N_{\rm pulsars}$ for two main reasons: 1) this configuration becomes increasingly optimistic for any future scenario, and 2) to retain a reasonable computational cost.

IV.4 Sensitivity to subdominant SGWBs

In this section, we forecast future sensitivity to subdominant contribution to the nHz SGWB. For presentation purposes, we consider the scenario in which the additional subdominant GWB contribution is described by a lognormal (LN) shape

$$h^{2}\Omega_{\rm CGW}(f)=10^{\alpha_{\rm LN}}e^{\left[-\ln^{2}(f/(10^{\gamma_{\rm LN}}{\rm Hz}))/(2\cdot 10^{2\beta_{\rm LN}})\right]}.$$ (32)

This should serve as a toy model we adopt to test our framework. Typical new physics scenarios of interest for PTA observations predict the existence of a blue tilted spectrum which reaches a peak close to the nHz frequencies, also to comply with $\Delta N_{\rm eff}$ upper bounds on cosmological (i.e., from the early universe) GWs, forcing the total primordial GW abundance $\Omega_{\rm GW}h^{2}<1.6\cdot 10^{-6}\left(\Delta N_{\rm eff}/0.28\right)$ at 95% C.L. Aghanim et al. (2020). Nonetheless, flatter spectra in the nHz range are predicted by cosmic string scenarios, as well as second order induced GWs from enhanced and flat curvature power spectra. In-depth analyses which assume specific spectral shapes motivated by new physics cosmological scenarios (see e.g., Afzal et al. (2023); Antoniadis et al. (2023c); Madge et al. (2023); Figueroa et al. (2023); Ellis et al. (2024)) will be presented elsewhere et al. (2024).

Again, to validate our framework in the case of multiple SGWB contributions, we generate a synthetic dataset as described in Sec. III.2 with the SKA10 configuration (50 pulsars and $T_{\rm obs}=10$yrs) of observation of an SGWB which takes contributions from two sources

$$h^{2}\Omega_{\rm GW}=h^{2}\Omega_{\rm GW}^{\text{PL}}+h^{2}\Omega_{\rm GW}^{\text{LN}}.$$ (33)

Collectively, the injected signal features five parameters, which are fixed to $\{\alpha_{*},n_{T},\alpha_{\rm LN},\beta_{\rm LN},\gamma_{\rm LN}\}=\{-7.2,2,-7,-1,-7.5\}$. In Fig. 8, we compare the FIM estimates to the state-of-the-art multi-components injection and inference. As in Fig. 1, one observes excellent agreement between both results, validating our framework even in the presence of multiple contributions to the SGWB.

We further explore the sensitivity to an LN subdominant contribution to the SGWB in Fig. 9. We simulate different scenarios in which the amplitude of the LN signal is varied, keeping fixed remaining parameters $\{\alpha_{*},n_{T},\beta_{\rm LN},\gamma_{\rm LN}\}=\{-7.2,2,-1,-7.5\}$. Uncertainties on PL parameters remain practically constant. This is because most of the injections feature a negligible contribution from the LN bump. Only large values of $\alpha_{\rm LN}$ the uncertainty grow slightly, due to the partial contamination around LN peak frequencies. On the other hand, we observe increasing precision on the LN parameters with larger $\alpha_{\rm LN}$. Interestingly, the relative uncertainty on the signal amplitude falls below $\sim 30\%$ around $\alpha_{\rm LN}\simeq-7.4$, indicating detection capability, when the LN amplitude crosses the effective sensitivity.

Figure 8: Same as Fig. 1, but injecting a SGWB from multiple sources: i) the dominant GW signal parametrized as a PL with $\alpha_{*}=-7.2$ and $n_{T}=2$ and ii) subdominant LN bump controlled by $\alpha_{\rm LN}=-7$, $\beta_{\rm LN}=-1$ and $\gamma_{\rm LN}=-7.5$.

Figure 9: Left panel: Injected SGWB including PL (purple, Eq. (30)) and LN (cyan, Eq. (32)). The grey solid lines show the effective sensitivity obtained with the future SKA10 configuration (with $N_{\rm pulsars}=100$ and $T_{\rm obs}=10$yrs). Different values of the amplitude are shown, varying $\alpha_{\rm LN}\in[-9,-7]$ with regular spacing. Right panel: Relative uncertainties on the model parameters, injecting the same values shown in the left panel. While uncertainties on PL parameters remain practically constant, we see improved accuracy on the LN parameters with larger $\alpha_{\rm LN}$.

V Conclusions

Pulsar timing array experiments provide an unprecedented opportunity to search for the existence of GWs in the nHZ frequency range. This will allow to constrain both astrophysical and cosmological sources of such signatures, with increasing precision. Current data report growing evidence for the existence of a signal with an amplitude $\Omega_{\rm GW}h^{2}\sim 10^{-9}$ at the best constrained frequencies $f\sim 1/(10{\rm yrs})$. In the future, longer observation time and the growing number of pulsars available will drastically expand PTA sensitivity.

In this work, we have demonstrated how to achieve fast estimation of detectability and parameter estimation for a stochastic GW signal which one could expect for a particular Pulsar Timing Array. Despite the simplistic assumptions adopted, we find that, compared to a full Bayesian analysis using the software ENTERPRISE on the EPTA DR2new and simulated data, this approach gives reliable results. Moreover, we have demonstrated that this method recovers an expected scaling of the PTA sensitivity with observation time and the number of pulsars. We can also study the spatial correlation in the pulsar’s pairs either using a binned consideration or Legendre decomposition. Finally, the formalism allows us to investigate multi-component stochastic GW signals of different strengths, going beyond the weak signal approximation so far adopted in the literature, providing fast and reliable forecasts on the constraining power of PTA configurations on both astrophysical and cosmological SGWBs. We release the code fastPTA which implements this framework.

If we neglect a pulsar term, extending this work to include the resolvable GW signals from individual sources (e.g., continuous GW (CGW) from inspiralling SMBHBs) would be quite simple. In this case, assuming that the CGW candidates are identified, we can include CGW in the transmission function, which will appear as a notch at the CGW frequency. To demonstrate this, we can write the CGW in the form $s=\sum_{k}a_{k},h_{k}(t,f_{\rm{CGW}},\rm{sky})$ (see Babak and Sesana (2012); Ellis et al. (2012)) which is similar to the timing (linearised) model and will give a contribution to the $G$ matrix Hazboun et al. (2019) and, therefore, to the transmission ${\cal T}$. Extension of this approach to CGWs with pulsar terms is somewhat more complicated and will be discussed elsewhere.

As promised, we want to say a few words about mimicking the chromatic noises in our approach. We can restore the dependence on the radio frequency, as it will become important for the simultaneous ultra-broad band observations. One possibility on how it can be done is presented in Hazboun et al. (2019) using the DMX model. In the case of Gaussian Process representation, we need to translate it to DMX-compatible representation. The DM variations through the DMX model appear in the transmission function.

We conclude by mentioning possible further future directions:

• Astrophysical models predict the SGWB comes from the superposition of unresolved signals. In this case, spectral fluctuations due to finite source numbers within frequency bins may appear in the observable band at a relatively "high" frequency. It would be interesting to estimate the future sensitivity as well as frequency resolutions to such spikes, which are not expected in cosmological scenarios.

• Throughout our work, we have assumed the dominant GWB to be stationary and isotropic. In particular, the effective sensitivity reported in our work assumes this background cannot be reduced. If the signal is of astrophysical origin, with sufficient sensitivity one should be able to resolve individual SMBH mergers and remove some CGW. This could improve the effective sensitivity to subdominant GWBs.

• With our code, one could also test the sensitivity to different correlation patterns beyond HD (e.g. motivated by modified theories of gravity).

• Finally, it would be interesting to extend this framework to include SGWB anisotropies (see e.g., Ali-Haïmoud et al. (2020, 2021); Cruz et al. (2024) for work in this direction), going beyond the weak signal approximation.