# Forecasting the sensitivity of Pulsar Timing Arrays to gravitational wave backgrounds

Stanislav Babak, Mikel Falxa, Gabriele Franciolini, Mauro Pieroni

► To cite this version:

Stanislav Babak, Mikel Falxa, Gabriele Franciolini, Mauro Pieroni. Forecasting the sensitivity of Pulsar Timing Arrays to gravitational wave backgrounds. Physical Review D, 2024, 110 (6), pp.063022. 10.1103/PhysRevD.110.063022. [hal-04557205]

HAL Id: hal-04557205

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

Submitted on 27 Jan 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Forecasting the sensitivity of pulsar timing arrays
to gravitational wave backgrounds

Stanislav Babak,<sup>1,*</sup> Mikel Falxa,<sup>2,1,†</sup> Gabriele Franciolini,<sup>3,‡</sup> and Mauro Pieroni<sup>3,§</sup>

<sup>1</sup>Astroparticule et Cosmologie, CNRS, Université Paris Cité, F-75013 Paris, France

<sup>2</sup>LPC2E, CNRS, 3 Av. de la Recherche Scientifique, 45071 Orléans, France

<sup>3</sup>CERN, Theoretical Physics Department, Esplanade des Particules 1, Geneva 1211, Switzerland

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

(Received 20 April 2024; accepted 30 July 2024; published 11 September 2024)

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 pulsar timing array 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.

DOI: 10.1103/PhysRevD.110.063022

# I. INTRODUCTION

Pulsar timing array (PTA) 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 [1,2], European Pulsar Timing Array (EPTA) [in combination with Indian Pulsar Timing Array (InPTA)] [3–5], Parkes Pulsar Timing Array (PPTA) [6–8], and Chinese Pulsar Timing Array (CPTA) [9] 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 [10].

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_{GW}(f) \propto f^{n_r}$, with $n_T \approx 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 inspiraling supermassive black hole binaries (SMBHBs) on circular orbits, whose frequency spectrum is characterized by a scaling law

$\Omega_{GW} \propto f^{2/3}$ [11]. Nevertheless, orbital eccentricity and stellar environment could affect the binary evolution that is coupled to the statistical fluctuation in the distribution of SMBHBs, which may play an important role and can lead to modified spectra (see, e.g., [12–21]). 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 [22–37], cosmic strings and domain walls [38–57], scalar-induced GWs generated from primordial fluctuations [58–82], tensor perturbations generated during inflation (see, e.g., [83,84]), and many others (see also [85–91]).

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

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Open access publication funded by CERN.

<footnote>*Contact author: stas@apc.in2p3.fr</footnote>

<footnote>†Contact author: falxa@apc.in2p3.fr</footnote>

<footnote>‡Contact author: gabriele.franciolini@cern.ch</footnote>

<footnote>§Contact author: mauro.pieroni@cern.ch</footnote>
resulting from the analysis of the recent PTA datasets [4], 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. [92] 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. [93], 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 Refs. [92,94–96], 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, [3]) and forecasted from mock datasets 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<fnref n="1" /> and define the Fourier domain data as

$$ \tilde{d}_i^k = \int_{-T_{\text{obs}}/2}^{T_{\text{obs}}/2} d_I(t) e^{-i2\pi f_k t} dt. \quad (1) $$

Given the finite duration of the data stream, we adopt a finite set of frequencies for the Fourier basis, $f_k = k/T_{\text{obs}}$, starting from the inverse of the observation time $T_{\text{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}. \quad (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 HD correlation pattern to be discussed in Sec. II B.

**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.<fnref n="2" /> Crucially, the drastic decrease only happens around the frequency $1/T_{\text{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_{\text{obs}}$ [92]. We model this transfer function as

$$ \mathcal{T}(f) \approx [1 + 1/(f T_{\text{obs}})]^{-6}. \quad (3) $$

An additional loss of the sensitivity is induced around $f = 1/\text{yr}$ from fitting the sky position and the proper motion and around $f = 2/\text{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.

### A. Noise model

Following the discussion in Ref. [92], we isolate a few key features of the noise budget. For each pulsar, the noise contribution comes schematically from

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

defined as follows:

<footnote n="1">¹In reality, we only need stationarity in the frequency band relevant to the targeted stochastic GW signal. See conclusions for a more comprehensive discussion of this point.</footnote>

<footnote n="2">²Sensitivity to the GWs with frequencies below $1/T_{\text{obs}}$ can be retained in delays compared to higher order spin-down terms of the pulsar timing model (see, e.g., [97,98]).</footnote>
1. White noise

The time of arrival of the signal shows an overall white noise (WN) 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^{\text{WN}}$, where $I$ is the index associated with each pulsar. This is typically controlled by the parameters EQUAD, ECORR, and EFAC in PTA analyses [2,4]. This can be described by a time delay power spectrum

$$P_I^{\text{WN}} = 2\sigma^2\Delta t, \qquad (5)$$

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

2. Red noise

The stochasticity in the pulsar's rotation causes an "achromatic" red noise (RN), 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^{\text{RN}} = A_I^{\text{RN}} \left(\frac{f}{f_r}\right)^{\gamma_I^{\text{RN}}}; \qquad (6)$$

where $f_r$ is a reference frequency often assumed to be 1/yr, and the amplitude $A_I^{\text{RN}}$ and the spectral index $\gamma_I^{\text{RN}}$ is intrinsic to each pulsar $I$.

3. 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 [99] 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. [100]) of the form

$$P_{n,I} = P_I^{\text{WN}} + A_I^{\text{RN}} \left(\frac{f}{f_r}\right)^{\gamma_I^{\text{RN}}} + P_I^{\text{DM,SV}}, \qquad (7)$$

where we arbitrarily fix the reference frequency at $f_r = f_{\text{yr}}$ and $\gamma_I^{\text{RN}} < 0$. $P_I^{\text{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 rediscuss this point in the concluding section.

As we will see, the simplified modeling summarized 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.

B. 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 d^2\Omega_{\hat{k}} \int_{-\infty}^{\infty} df [\tilde{h}_{+}(f, \hat{k}) e_{ab}^{+}(\hat{k}) + \tilde{h}_{\times}(f, \hat{k}) e_{ab}^{\times}(\hat{k})] e^{i2\pi f(t - \hat{k} \cdot \vec{x}/c)}, \qquad (8)$$

where we introduced the polarization tensors $e_{ab}^{+, \times}(\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'}^{*}(f', \hat{k}') \rangle = \frac{1}{16\pi} S_h(f) \delta(f-f') \delta_{PP'} \delta^2(\hat{k}-\hat{k}'), \qquad (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., [101])

$$\Omega_{\text{GW}} h^2 = \frac{h^2}{\rho_c} \frac{d\rho_{\text{GW}}}{d\log f} \equiv \frac{2\pi^2 f^3}{3H_0^2/h^2} S_h, \qquad (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 \text{ Gyr})$ is the Hubble parameter today. In the last equality, we have associated $\Omega_{\text{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 [92]

$$\tilde{\Delta}t_I(f) = \int d^2\Omega_{\hat{k}} \sum_P R_I^P(f, \hat{k}) \tilde{h}_P(f; \hat{k}), \qquad (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}) = \frac{e_{ab}^P(\hat{k})}{i4\pi f} \frac{\hat{p}_I^a \hat{p}_I^b}{1 + \hat{p}_I \cdot \hat{k}} (1 - e^{-i2\pi fD(1+\hat{k} \cdot \hat{p}_I)/c}). \quad (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_* \approx 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_{\text{obs}} \sim \text{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

$$\langle \tilde{\Delta}t_I \tilde{\Delta}t_J^* \rangle = \frac{1}{2} \delta(f - f') C_{h,IJ} = \frac{1}{2} \delta(f - f') R_{IJ} S_h(f), \quad (13)$$

where we introduced [92]

$$R_{IJ} = \chi_{IJ} \cdot \mathcal{R}(f) [T_I(f)T_J(f)T_{IJ}/T_{\text{obs}}]^{1/2}, \quad (14)$$

and we conveniently defined the frequency dependent, sky-averaged, quadratic response function $\mathcal{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., [102]), 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} = \min[T_I, T_J]$.

## 1. The Hellings-Downs correlation

We introduced the HD function [10] 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}, \quad (15)$$

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

To estimate the sensitivity to the HD angular correlation function, it is convenient to decompose it on a suitable basis. Following Ref. [3], we adopt the following templates:

(1) 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}), \quad (16)$$

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

(2) HD expansion in Legendre polynomials: following [3,103], we expand the HD function as

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

Using the standard normalization of the Legendre polynomials, the coefficients are found by [103,104]

$$\begin{aligned} a_{\ell} &= \frac{2\ell + 1}{2} \int_{-1}^{1} \chi_{\text{HD}}(x) P_{\ell}(x) \mathrm{d}x \\ &= \frac{3}{2(\ell + 2)(\ell + 1)\ell(\ell - 1)}. \end{aligned} \quad (18)$$

for $\ell \ge 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$.

## 2. Effective sensitivity

We compute the SNR by

$$\text{SNR}^2 = \sum_{f_k} [C_{IJ}^{-1} C_{KL}^{-1} R_{JK} R_{LI}] S_h^2, \quad (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 [92]

$$S_{\text{eff}}(f) = (C_{IJ}^{-1} C_{KL}^{-1} R_{JK} R_{LI})^{-1/2}, \quad (20)$$

in such a way that the SNR computation takes the familiar form $\text{SNR}^2 = \sum_k (S_h/S_{\text{eff}})^2 = T_{\text{obs}} \int df (S_h/S_{\text{eff}})^2$.

In the weak signal limit (which is in practice the limit in which $RS_h \ll P_n$), one can expand the correlation matrix above to find

$$S_{\text{eff}}(f) = \left( \sum_{I \neq J} \frac{T_{IJ}}{T_{\text{obs}}} \frac{T_I(f)T_J(f)\chi_{IJ}^2}{P_{n,I}(f)P_{n,J}(f)/\mathcal{R}^2(f)} \right)^{-1/2}, \quad (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. [92] and shown in Refs. [105,106] to report current PTA sensitivities. As we will see, when considering future

<footnote n="3">³Notice that for a small number of pulsars, a sensible choice would be to divide the angular range unevenly and choose bins containing an equal number of pulsar pairs [3].</footnote>
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.

### C. 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 [107,108]

$$ -\ln \mathcal{L}(\tilde{d}|\theta) = \text{const.} + \sum_{k,IJ} [\ln[C_{IJ}(f_k, \theta)] + \tilde{d}_I^k C_{IJ}^{-1}(f_k, \theta) \tilde{d}_J^{k*} / T_{\text{obs}}], \quad (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. We note the absence of a factor of 2 from the likelihood definition coming from our adoption of positive frequencies only.

### 1. Fisher information matrix estimates

Here we briefly summarize 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} = -\left. \frac{\partial^2 \log \mathcal{L}}{\partial \theta^\alpha \partial \theta^\beta} \right|_{\theta=\theta_0} = \sum_k \text{Tr} \left[ C^{-1} \frac{\partial C}{\partial \theta^\alpha} C^{-1} \frac{\partial C}{\partial \theta^\beta} \right]_{\theta=\theta_0}, \quad (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} [C^{-1} - C^{-1} \tilde{d}^k \tilde{d}^{k*} C^{-1} / T_{\text{obs}}] = 0, \quad (24) $$

which is solved by $C(f_k, \theta_0) = \tilde{d}^k \tilde{d}^{k*} / T_{\text{obs}}$. In practice, the discrete sum over finite frequencies can be replaced with a continuous integral over the frequency range. Keeping all indices fully expressed, one obtains

$$ F_{\alpha\beta} = \sum_{f_k} C_{IJ}^{-1} C_{KL}^{-1} \frac{\partial(R_{JK}S_h)}{\partial\theta^\alpha} \frac{\partial(R_{LI}S_h)}{\partial\theta^\beta}. \quad (25) $$

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

### 2. Monte Carlo Markov chain 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 using EMCEE [109] 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 [110]. At the moment, PTA datasets are limited by radiometer noise. The SKA will provide high-precision pulsar timing measurement with uncertainties below ~100 ns [110], making it roughly 10 times better than current generation telescopes [111]. 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:

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

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

### A. Fake PTA data and time-domain Bayesian analysis

#### 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 $\tilde{d} = [d_1(t), d_2(t), \dots, d_{N_{\text{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 \rightarrow d_I - \sum_k \alpha_k \vec{M}_{I,k}, \quad (26) $$

where $M_{I,k}$ is a design matrix (representing the derivative of the timing model with respect to 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 [112] 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), \quad (27)$$

where $X_{i,I}$, $Y_{i,I} \sim 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_{\text{pulsars}}$ dimensional multivariate normal distribution $\vec{X}_i, \vec{Y}_i \sim 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.

## 2. Generating fake data

We simulate the data in the time domain.<fnref n="4" /> 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 (27) 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:

(1) Offset: $\vec{M}_0 \propto \vec{1}$.

(2) Spin rate: $\vec{M}_1 \propto \vec{t}_I$.

(3) Quadratic spin-down rate: $\vec{M}_2 \propto \vec{t}_I^2$.

(4) RA: $\vec{M}_3 \propto \sin(2\pi f_{yr} \vec{t}_I)$.

(5) 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 Right

Ascension (RA) and Declination (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 noninformative improper prior. This marginalization is responsible for the transmission function described in (3) [92].

## B. Mock datasets

We simulated two PTA mock datasets:

(1) EPTA20: built from the 25 pulsars of the EPTA DR2new 10.4 yr dataset with the same noise properties (as in [113]) 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 reinjected according to their maximum likelihood values in [113].

(2) SKA10: SKA-like PTA with 50 pulsars and 10 yr observations. Every pulsar is simulated with a 2-week cadence of observation, $\sigma = 100$ ns 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 summarized in Table I. 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 [1,3], potentially originating from the population of SMBHBs in the Universe. This signal will act as foreground, and we

<table><caption>TABLE I. Summary of mock datasets properties. $N_{\text{pulsars}}$ is the number of pulsars, $\nu$ are the observing frequency subbands, $T_{\text{obs}}$ is the time of observation, $\sigma$ is the level of white noise, $\Delta t$ is the cadence, TM is the timing model where $\text{TM}_5$ is the five-component model described in Sec. III A 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. [113] correspond to the choice of parameters consistent with ERPTA DR2new analysis.</caption><thead><tr><th></th><th>EPTA20</th><th>SKA10</th></tr></thead><tbody><tr><td>$N_{\text{pulsars}}$</td><td>25 + 25</td><td>50</td></tr><tr><td>$\nu$ [GHz]</td><td>[113]</td><td>1.4–3</td></tr><tr><td>$T_{\text{obs}}$ [yr]</td><td>20.8 + 10.4</td><td>10</td></tr><tr><td>$\sigma$ [s]</td><td>~10<sup>-6</sup></td><td>10<sup>-7</sup></td></tr><tr><td>$\Delta t$ [days]</td><td>~3</td><td>14</td></tr><tr><td>TM</td><td>TM<sub>5</sub> + [113]</td><td>TM<sub>5</sub></td></tr><tr><td>$\log_{10} A_{\text{RN}}$</td><td>[113]</td><td>[-17, -13]</td></tr><tr><td>$\gamma_{\text{RN}}$</td><td>[113]</td><td>[1, 5]</td></tr><tr><td>$\log_{10} A_{\text{DM}}$</td><td>[113]</td><td>[-17, -13]</td></tr><tr><td>$\gamma_{\text{DM}}$</td><td>[113]</td><td>[1, 5]</td></tr><tr><td>$\log_{10} A_{\text{SV}}$</td><td>[113]</td><td>–</td></tr><tr><td>$\gamma_{\text{SV}}$</td><td>[113]</td><td>–</td></tr></tbody></table>

<footnote n="4">⁴The time-domain simulations were made using https://github.com/mfalxa/fakepta.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Comparison between the measurement uncertainties estimated using FIM (blue) and time-domain Bayesian parameter estimation (red, see Sec. III A 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.</figcaption></figure>

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 D).

### C. Time-domain Bayesian analysis

The timing model marginalized Gaussian likelihood used in time-domain PTA analysis [112] can be expressed as [analogously to the previously defined likelihood (22)]

$$ -\ln \mathcal{L}(d|\theta) = \text{const.} + \sum_{IJ} [\ln(\det\{C_{IJ}\}) + d_I(t)C_{IJ}^{-1}(\theta)d_J^T(t)], \quad (28) $$

where $C_{IJ} = C_I \delta_{IJ} + C_{h,IJ}$ has a block structure of size $N_{\text{pulsars}} \times N_{\text{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^{\text{TM}} + C_I^{\text{RN}} + C_I^{\text{DM}} + C_I^{\text{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 [114,115] and sampler PTMCMC [116]. ENTERPRISE provides computation of the likelihood and prior based on the chosen model. All noise components are modeled as Gaussian processes with incomplete Fourier basis as eigenfunctions [112]. 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_{\text{GW}}^{\text{PL}}(f) = A_{\text{GWB}} \left( \frac{f}{f_{\text{yr}}} \right)^{n_T}, \quad (29) $$

with $A_{\text{GWB}} \approx 6.3 \times 10^{-8}$ and $n_T \approx 2$. For definiteness, these parameters are fixed to the maximum likelihood values obtained by the EPTA [3], which is in good agreement with the other observations [106]. 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} \approx 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_{\text{GWB}})$.

Assuming this SGWB, we obtain SNR = 4.5 with the current EPTA (DRNew) pulsar configuration, which is compatible with the observations [3]. On the other hand, we obtain SNR = 184 for EPTA20 and SNR = 292 for SKA10.

To validate our framework based on the approximated likelihood introduced in Sec. II C, we compare our results to the ones obtained with state-of-the-art techniques described in Sec. III C. 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 [4]. 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_{\text{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.
<figure><img src="image_4.png" /><figcaption>FIG. 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 (29).</figcaption></figure>

### A. Future effective sensitivity

Using the metrics defined in Sec. II B 2, we can forecast the effective sensitivity achieved by current and future experiments. In particular, we plot $S_{\text{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 \text{ 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_{\text{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. (29). 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.

### B. 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.<fnref n="5" /> The parameters determining the noise of each pulsar are sampled from distributions built out of the currently observed ones [4].

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.

Focusing on the evolution of SNR, we observe that at small observation times, it grows rapidly as SNR $\sim T_{\text{obs}}^3$. This can be explained as follows (see analytical considerations reported in Ref. [117]). Assuming WN dominates at the initially relatively high frequencies $f_{\text{min}} \sim 1/T_{\text{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 prefactor 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_{\text{obs}}}$

<footnote n="5">⁵For both Figs. 3 and 4, we generated 30 different realizations for each value of $T_{\text{obs}}$ and $N_{\text{pulsars}}$, to estimate the scatter induced by randomly choosing noise parameter for each pulsar.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Left panel: scaling of the SNR as a function of the observation time $T_{\text{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_{\text{obs}}$ regime is dictated by our choice of the spectrum with $n_T = 2$, while the high-$T_{\text{obs}}$ regime is universal. Right panel: scaling of the uncertainties on the amplitude and tilt of the SGWB as a function of observation time.</figcaption></figure>

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_{\text{obs}} \sim 10$ yrs (see also [118]). 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 autocorrelation 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_{\text{pulsars}}}$.

Had we only included the off-diagonal components, featuring the HD correlation, we would have obtained a scaling SNR $\sim (\sum_{I \neq J} \chi_{IJ}^2)^{1/2}$, which behaves as SNR $\sim [N_{\text{pulsars}}(N_{\text{pulsars}} - 1)]^{1/2} \sim N_{\text{pulsars}}$.

### C. 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., [119]). For this purpose, we allow $\chi_{IJ}$ to vary, and, adopting the HD curve parametrizations described in Sec. II B 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 parametrization (16) found

<figure><img src="image_6.png" /><figcaption>FIG. 4. Scaling of the SNR (left panel) and uncertainties (right panel) as a function of the number of pulsars, assuming fixed $T_{\text{obs}} = 20$ yrs. The error band indicates the standard deviation of both results induced by producing multiple realizations of the pulsar sample.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 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 seven equally spaced bins. Uncertainties are comparable with the one observed found in EPTADR2 analyses [3]. 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).</figcaption></figure>

using FIM. For this analysis, we considered seven 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% CL in the central bins. These results are in good agreement with the posterior distributions reported by the EPTA collaboration [3], 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_{II}$ 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% CL.

## 1. 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. [120,121], 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}}. \qquad (30) $$

<figure><img src="image_8.png" /><figcaption>FIG. 6. Same as Fig. 5, but assuming the SKA10 future configuration.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 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_{\text{pulsars}}$. We adopt $N_{\text{modes}} = 30$. Left panel: zero noise limit. The zeroth and first orders are zero, and the corresponding uncertainty scales as $\propto 1/N_{\text{pulsars}}$ (indicated by the dashed black line). The higher-order pieces reach the plateau predicted by cosmic variance in Eq. (30). See the discussion in the main text. Right panel: same as the left panel, but includes noise. The plateau at small values of $N_{\text{pulsars}}$ is reached when information is small enough that our prior on $\sigma_{a_\ell} = 1$ becomes dominant.</figcaption></figure>

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 (N_{\text{modes}})^{-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_{\text{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_{\text{pulsars}}$, without showing any significant departure at high $N_{\text{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 (30). The plateau is reached at an increasing number of pulsars for higher $\ell \ge 2$ due to the increasingly larger $N_{\text{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_{\text{pulsars}}$, where their posterior distribution is dominated by the prior (assumed to be Gaussian with $\sigma_{a_\ell} = 1$). When $N_{\text{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_{\text{pulsars}}$ scaling, while starting from the lowest orders in $\ell$ (especially $\ell = 2, 3$) we notice a flattening of the curve already from $N_{\text{pulsars}} \sim \mathcal{O}(\text{few} \times 10^2)$. We do not push the computations to even larger $N_{\text{pulsars}}$ for two main reasons: (1) this configuration becomes increasingly optimistic for any future scenario and (2) to retain a reasonable computational cost.

## D. 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 log-normal (LN) shape

$$h^2 \Omega_{\text{CGW}}(f) = 10^{\alpha_{\text{LN}}} e^{[-\ln^2(f/(10^{\alpha_{\text{LN}}}) \text{ Hz}))/2 \times 10^{2\beta_{\text{LN}}}}. \quad (31)$$

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

<figure><img src="image_10.png" /><figcaption>FIG. 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_{\text{LN}} = -7$, $\beta_{\text{LN}} = -1$, and $\gamma_{\text{LN}} = -7.5$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIG. 9. Left panel: injected SGWB including PL [purple, Eq. (29)] and LN [cyan, Eq. (31)]. The gray solid lines show the effective sensitivity obtained with the future SKA10 configuration (with $N_{\text{pulsars}} = 100$ and $T_{\text{obs}} = 10$ yrs). Different values of the amplitude are shown, varying $\alpha_{\text{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_{\text{LN}}$.</figcaption></figure>

spectrum which reaches a peak close to the nHz frequencies, also to comply with $\Delta N_{\text{eff}}$ upper bounds on cosmological (i.e., from the early Universe) GWs, forcing the total primordial GW abundance $\Omega_{\text{GW}}h^2 < 1.6 \times 10^{-6}$ ($\Delta N_{\text{eff}}/0.28$) at 95% CL [122]. 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 that assume specific spectral shapes motivated by new physics cosmological scenarios (see, e.g., [5, 18, 86, 87, 123]) will be presented elsewhere [124].

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

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

Collectively, the injected signal features five parameters, which are fixed to $\{\alpha_*, n_T, \alpha_{\text{LN}}, \beta_{\text{LN}}, \gamma_{\text{LN}}\} = \{-7.2, 2, -7, -1, -7.5\}$. In Fig. 8, we compare the FIM estimates to the state-of-the-art multicomponents 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_{\text{LN}}, \gamma_{\text{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_{\text{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_{\text{LN}}$. Interestingly, the relative uncertainty on the signal amplitude falls below ~30% around $\alpha_{\text{LN}} \approx -7.4$, indicating detection capability, when the LN amplitude crosses the effective sensitivity.

## V. CONCLUSIONS

PTA 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_{\text{GW}}h^2 \sim 10^{-9}$ at the best-constrained frequencies $f \sim 1/(10 \text{ 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 PTA. 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.<fnref n="6" /> We can also study the spatial correlation in the pulsar's pairs either using a binned consideration or Legendre decomposition. Finally, the

<footnote n="6">⁶For previous studies forecasting the scaling of the uncertainties with future PTA configurations, see [118,125].</footnote>
formalism allows us to investigate multicomponent 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 inspiraling 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_{\text{CGW}}, \text{sky})$ (see [126,127]) which is similar to the timing (linearized) model and will give a contribution to the $G$ matrix [92] and, therefore, to the transmission $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 [92] using the piecewise constant model for dispersion measure (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. Finally, we stress again that in the FIM we make two assumptions that render our approach insensitive to potential noise mismodeling. First, we describe the pulsar noises using their best parameters measured by the EPTA collaboration or from time-domain inference applied to mock data simulations. Moreover, we do not include noise parameters in the FIM and therefore do not capture potential degeneracies between signal and noise. The comparison in Sec. III shows that our results are sufficiently robust despite neglecting these issues. Furthermore, these effects become less and less relevant as observations will enter in the signal dominated regime.

We conclude by mentioning possible further future directions:

(1) 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.

(2) 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.

(3) With our code, one could also test the sensitivity to different correlation patterns beyond HD, e.g., motivated by modified theories of gravity [119,128–130].

(4) It would be interesting to extend this framework to include SGWB anisotropies (see, e.g., [131–133] for work in this direction), going beyond the weak signal approximation.

(5) Finally, one could include information from astrometry [134–137], to forecast how this would improve the constraining power of PTA observations (see [138]).

## ACKNOWLEDGMENTS

We thank Bruce Allen, Paul Frederik Depta, and Valerie Domcke for valuable discussions on the HD cosmic variance. M. P. thanks Rutger van Haasteren for the useful discussion on PTA data analysis. S. B. acknowledges support from funding from the French National Research Agency (Grant No. ANR21-CE31-0026, project MBH waves), the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 101007855. M. P. acknowledges the hospitality of Imperial College London, which provided office space during some parts of this project.

[1] G. Agazie et al. (NANOGrav Collaboration), Astrophys. J. Lett. **951**, L8 (2023).

[2] G. Agazie et al. (NANOGrav Collaboration), Astrophys. J. Lett. **951**, L9 (2023).

[3] J. Antoniadis et al. (EPTA, InPTA: Collaborations), Astron. Astrophys. **678**, A50 (2023).

[4] J. Antoniadis et al. (EPTA Collaboration), Astron. Astrophys. **678**, A48 (2023).

[5] J. Antoniadis et al. (EPTA Collaboration), Astron. Astrophys. **685**, A94 (2024).

[6] D. J. Reardon et al., Astrophys. J. Lett. **951**, L6 (2023).

[7] A. Zic et al., Pub. Astron. Soc. Aust. **40**, e049 (2023).
[8] D. J. Reardon et al., Astrophys. J. Lett. **951**, L7 (2023).

[9] H. Xu et al., Res. Astron. Astrophys. **23**, 075024 (2023).

[10] R. w. Hellings and G. s. Downs, Astrophys. J. Lett. **265**, L39 (1983).

[11] E. S. Phinney, arXiv:astro-ph/0108028.

[12] A. Sesana, A. Vecchio, and C. N. Colacino, Mon. Not. R. Astron. Soc. **390**, 192 (2008).

[13] B. Kocsis and A. Sesana, Mon. Not. R. Astron. Soc. **411**, 1467 (2011).

[14] L. Z. Kelley, L. Blecha, and L. Hernquist, Mon. Not. R. Astron. Soc. **464**, 3131 (2017).

[15] D. Perrodin and A. Sesana, Astrophysics and Space Science Library **457**, 95 (2018).

[16] J. Ellis, M. Fairbairn, G. Hütsi, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe, Astron. Astrophys. **676**, A38 (2023).

[17] G. Agazie et al. (NANOGrav Collaboration), Astrophys. J. Lett. **952**, L37 (2023).

[18] A. Afzal et al. (NANOGrav Collaboration), Astrophys. J. Lett. **951**, L11 (2023).

[19] A. Ghoshal and A. Strumia, J. Cosmol. Astropart. Phys. 02 (2024) 054.

[20] M. Bonetti, A. Sesana, E. Barausse, and F. Haardt, Mon. Not. R. Astron. Soc. **477**, 2599 (2018).

[21] M. Aghaie, G. Armando, A. Dondarini, and P. Panci, Phys. Rev. D **109**, 103030 (2024).

[22] Z. Arzoumanian et al. (NANOGrav Collaboration), Phys. Rev. Lett. **127**, 251302 (2021).

[23] X. Xue et al., Phys. Rev. Lett. **127**, 251303 (2021).

[24] Y. Nakai, M. Suzuki, F. Takahashi, and M. Yamada, Phys. Lett. B **816**, 136238 (2021).

[25] P. Di Bari, D. Marfatia, and Y.-L. Zhou, J. High Energy Phys. 10 (2021) 193.

[26] A. S. Sakharov, Y. N. Eroshenko, and S. G. Rubin, Phys. Rev. D **104**, 043005 (2021).

[27] S.-L. Li, L. Shao, P. Wu, and H. Yu, Phys. Rev. D **104**, 043510 (2021).

[28] A. Ashoorioon, K. Rezazadeh, and A. Rostami, Phys. Lett. B **835**, 137542 (2022).

[29] M. Benetti, L. L. Graef, and S. Vagnozzi, Phys. Rev. D **105**, 043520 (2022).

[30] J. Barir, M. Geller, C. Sun, and T. Volansky, Phys. Rev. D **108**, 115016 (2023).

[31] M. Hindmarsh and J. Kume, J. Cosmol. Astropart. Phys. 04 (2023) 045.

[32] P. Athron, A. Fowlie, C.-T. Lu, L. Morris, L. Wu, Y. Wu, and Z. Xu, Phys. Rev. Lett. **132**, 221001 (2024).

[33] Y. Gouttenoire and T. Volansky, arXiv:2305.04942.

[34] I. Baldes, M. Dichtl, Y. Gouttenoire, and F. Sala, arXiv:2306.15555.

[35] S.-P. Li and K.-P. Xie, Phys. Rev. D **108**, 055018 (2023).

[36] P. Di Bari and M. H. Rahat, arXiv:2307.03184.

[37] Y. Gouttenoire, Phys. Rev. Lett. **131**, 171404 (2023).

[38] J. Ellis and M. Lewicki, Phys. Rev. Lett. **126**, 041304 (2021).

[39] S. Datta, A. Ghosal, and R. Samanta, J. Cosmol. Astropart. Phys. 08 (2021) 021.

[40] R. Samanta and S. Datta, J. High Energy Phys. 05 (2021) 211.

[41] W. Buchmuller, V. Domcke, and K. Schmitz, Phys. Lett. B **811**, 135914 (2020).

[42] S. Blasi, V. Brdar, and K. Schmitz, Phys. Rev. Lett. **126**, 041305 (2021).

[43] S. Ramazanov, E. Babichev, D. Gorbunov, and A. Vikman, Phys. Rev. D **105**, 063530 (2022).

[44] E. Babichev, D. Gorbunov, S. Ramazanov, and A. Vikman, J. Cosmol. Astropart. Phys. 04 (2022) 028.

[45] M. Gorghetto, E. Hardy, and H. Nicolaescu, J. Cosmol. Astropart. Phys. 06 (2021) 034.

[46] W. Buchmuller, V. Domcke, and K. Schmitz, J. Cosmol. Astropart. Phys. 12 (2021) 006.

[47] J. J. Blanco-Pillado, K. D. Olum, and J. M. Wachter, Phys. Rev. D **103**, 103512 (2021).

[48] R. Z. Ferreira, A. Notari, O. Pujolas, and F. Rompineve, J. Cosmol. Astropart. Phys. 02 (2023) 001.

[49] H. An and C. Yang, Phys. Rev. D **109**, 123508 (2024).

[50] Z.-Y. Qiu and Z.-H. Yu, Chin. Phys. C **47**, 085104 (2023).

[51] Z.-M. Zeng, J. Liu, and Z.-K. Guo, Phys. Rev. D **108**, 063005 (2023).

[52] S. F. King, D. Marfatia, and M. H. Rahat, Phys. Rev. D **109**, 035014 (2024).

[53] E. Babichev, D. Gorbunov, S. Ramazanov, R. Samanta, and A. Vikman, Phys. Rev. D **108**, 123529 (2023).

[54] N. Kitajima, J. Lee, K. Murai, F. Takahashi, and W. Yin, Phys. Lett. B **851**, 138586 (2024).

[55] B. Barman, D. Borah, S. Jyoti Das, and I. Saha, J. Cosmol. Astropart. Phys. 10 (2023) 053.

[56] G. Servant and P. Simakachorn, Phys. Rev. D **108**, 123516 (2023).

[57] Y. Gouttenoire and E. Vitagliano, arXiv:2306.17841.

[58] V. Vaskonen and H. Veermäe, Phys. Rev. Lett. **126**, 051303 (2021).

[59] V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. Lett. **126**, 041303 (2021).

[60] N. Bhaumik and R. K. Jain, Phys. Rev. D **104**, 023531 (2021).

[61] K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. Lett. **126**, 131301 (2021).

[62] K. Kohri and T. Terada, Phys. Lett. B **813**, 136040 (2021).

[63] G. Domènech and S. Pi, Sci. China Phys. Mech. Astron. **65**, 230411 (2022).

[64] R. Namba and M. Suzuki, Phys. Rev. D **102**, 123527 (2020).

[65] S. Sugiyama, V. Takhistov, E. Vitagliano, A. Kusenko, M. Sasaki, and M. Takada, Phys. Lett. B **814**, 136097 (2021).

[66] Z. Zhou, J. Jiang, Y.-F. Cai, M. Sasaki, and S. Pi, Phys. Rev. D **102**, 103527 (2020).

[67] J. Lin, S. Gao, Y. Gong, Y. Lu, Z. Wang, and F. Zhang, Phys. Rev. D **107**, 043517 (2023).

[68] K. Rezazadeh, Z. Teimoori, S. Karimi, and K. Karami, Eur. Phys. J. C **82**, 758 (2022).

[69] M. Kawasaki and H. Nakatsuka, J. Cosmol. Astropart. Phys. 05 (2021) 023.

[70] W. Ahmed, M. Junaid, and U. Zubair, Nucl. Phys. B**984**, 115968 (2022).

[71] Z. Yi and Q. Fei, Eur. Phys. J. C **83**, 82 (2023).

[72] Z. Yi, J. Cosmol. Astropart. Phys. 03 (2023) 048.

[73] V. Dandoy, V. Domcke, and F. Rompineve, SciPost Phys. Core **6**, 060 (2023).
[74] J.-X. Zhao, X.-H. Liu, and N. Li, Phys. Rev. D **107**, 043515 (2023).

[75] L. Liu, Z.-C. Chen, and Q.-G. Huang, J. Cosmol. Astropart. Phys. **11** (2023) 071.

[76] G. Ferrante, G. Franciolini, A. Iovino, Jr., and A. Urbano, J. Cosmol. Astropart. Phys. **06** (2023) 057.

[77] S. Wang, Z.-C. Zhao, J.-P. Li, and Q.-H. Zhu, Phys. Rev. Res. **6**, L012060 (2024).

[78] S. Wang, Z.-C. Zhao, and Q.-H. Zhu, Phys. Rev. Res. **6**, 013207 (2024).

[79] Y. Cai, M. Zhu, and Y.-S. Piao, Phys. Rev. Lett. **133**, 021001 (2024).

[80] G. Franciolini, A. Iovino, Junior., V. Vaskonen, and H. Veermae, Phys. Rev. Lett. **131**, 201401 (2023).

[81] S. Balaji, G. Domènech, and G. Franciolini, J. Cosmol. Astropart. Phys. **10** (2023) 041.

[82] L. Liu, Z.-C. Chen, and Q.-G. Huang, Phys. Rev. D **109**, L061301 (2024).

[83] X. Niu and M. H. Rahat, Phys. Rev. D **108**, 115023 (2023).

[84] S. Vagnozzi, J. High Energy Astrophys. **39**, 81 (2023).

[85] G. Franciolini, D. Racco, and F. Rompineve, Phys. Rev. Lett. **132**, 081001 (2024).

[86] E. Madge, E. Morgante, C. Puchades-Ibáñez, N. Ramberg, W. Ratzinger, S. Schenk, and P. Schwaller, J. High Energy Phys. **10** (2023) 171.

[87] D. G. Figueroa, M. Pieroni, A. Ricciardone, and P. Simakachorn, Phys. Rev. Lett. **132**, 171002 (2024).

[88] J. Garcia-Bellido, A. Papageorgiou, M. Peloso, and L. Sorbo, J. Cosmol. Astropart. Phys. **01** (2024) 034.

[89] K. Murai and W. Yin, J. High Energy Phys. **10** (2023) 062.

[90] R. A. Konoplya and A. Zhidenko, Phys. Lett. B **853**, 138685 (2024).

[91] C. Smarra et al. (EPTA Collaboration), Phys. Rev. Lett. **131**, 171001 (2023).

[92] J. S. Hazboun, J. D. Romano, and T. L. Smith, Phys. Rev. D **100**, 104028 (2019).

[93] https://github.com/Mauropieroni/fastPTA/.

[94] M. Anholm, S. Ballmer, J. D. E. Creighton, L. R. Price, and X. Siemens, Phys. Rev. D **79**, 084030 (2009).

[95] S. J. Chamberlin, J. D. E. Creighton, X. Siemens, P. Demorest, J. Ellis, L. R. Price, and J. D. Romano, Phys. Rev. D **91**, 044048 (2015).

[96] P. A. Rosado, A. Sesana, and J. Gair, Mon. Not. R. Astron. Soc. **451**, 2417 (2015).

[97] W. DeRocco and J. A. Dror, Phys. Rev. Lett. **132**, 101403 (2024).

[98] W. DeRocco and J. A. Dror, Phys. Rev. D **108**, 103011 (2023).

[99] A. Chalumeau et al. (EPTA Collaboration), Mon. Not. R. Astron. Soc. **509**, 5538 (2021).

[100] R. N. Caballero et al. (EPTA Collaboration), Mon. Not. R. Astron. Soc. **457**, 4421 (2016).

[101] C. Caprini and D. G. Figueroa, Classical Quantum Gravity **35**, 163001 (2018).

[102] A. Kehagias and A. Riotto, J. Cosmol. Astropart. Phys. **06** (2024) 059.

[103] J. Gair, J. D. Romano, S. Taylor, and C. M. F. Mingarelli, Phys. Rev. D **90**, 082001 (2014).

[104] J. D. Romano and B. Allen, arXiv:2308.05847.

[105] G. Agazie et al. (NANOGrav Collaboration), Astrophys. J. Lett. **951**, L10 (2023).

[106] G. Agazie et al. (International Pulsar Timing Array Collaboration), Astrophys. J. **966**, 105 (2024).

[107] C. R. Contaldi, M. Pieroni, A. I. Renzini, G. Cusin, N. Karnesis, M. Peloso, A. Ricciardone, and G. Tasinato, Phys. Rev. D **102**, 043502 (2020).

[108] J. R. Bond, A. H. Jaffe, and L. Knox, Phys. Rev. D **57**, 2117 (1998).

[109] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publ. Astron. Soc. Pac. **125**, 306 (2013).

[110] G. Janssen et al., Proc. Sci., AASKA14 (2015) 037 [arXiv: 1501.00127].

[111] T. J. W. Lazio, Classical Quantum Gravity **30**, 224011 (2013).

[112] R. van Haasteren and M. Vallisneri, Phys. Rev. D **90**, 104012 (2014).

[113] J. Antoniadis et al. (EPTA, InPTA Collaboration), Astron. Astrophys. **678**, A49 (2023).

[114] J. A. Ellis, M. Vallisneri, S. R. Taylor, and P. T. Baker, Enterprise: Enhanced numerical toolbox enabling a robust pulsar inference suite, Zenodo, 10.5281/zenodo.4059815 (2020).

[115] S. R. Taylor, P. T. Baker, J. S. Hazboun, J. Simon, and S. J. Vigeland, enterprise_extensions, https://github.com/nanograv/enterprise-extensions (2021), v2.3.3.

[116] J. Ellis and R. van Haasteren, jellis18/ptmcmcsampler: Official release, 10.5281/zenodo.1037579 (2017).

[117] X. Siemens, J. Ellis, F. Jenet, and J. D. Romano, Classical Quantum Gravity **30**, 224015 (2013).

[118] N. S. Pol et al. (NANOGrav Collaboration), Astrophys. J. Lett. **911**, L34 (2021).

[119] Q. Liang, M.-X. Lin, and M. Trodden, J. Cosmol. Astropart. Phys. **11** (2023) 042.

[120] E. Roebber and G. Holder, Astrophys. J. **835**, 21 (2017).

[121] B. Allen, Phys. Rev. D **107**, 043018 (2023).

[122] N. Aghanim et al. (Planck Collaboration), Astron. Astrophys. **641**, A6 (2020); **652**, C4(E) (2021).

[123] J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe, Phys. Rev. D **109**, 023522 (2024).

[124] C. C. et al. (to be published).

[125] A. R. Kaiser, N. S. Pol, M. A. McLaughlin, S. Chen, J. S. Hazboun, L. Z. Kelley, J. Simon, S. R. Taylor, S. J. Vigeland, and C. A. Witt, Astrophys. J. **938**, 115 (2022).

[126] S. Babak and A. Sesana, Phys. Rev. D **85**, 044034 (2012).

[127] J. A. Ellis, X. Siemens, and J. D. E. Creighton, Astrophys. J. **756**, 175 (2012).

[128] R. C. Bernardo and K.-W. Ng, Phys. Rev. D **107**, L101502 (2023).

[129] R. C. Bernardo and K.-W. Ng, Phys. Rev. D **109**, L101502 (2024).

[130] R. C. Bernardo and K.-W. Ng, arXiv:2310.07537.

[131] Y. Ali-Haïmoud, T. L. Smith, and C. M. F. Mingarelli, Phys. Rev. D **102**, 122005 (2020).

[132] Y. Ali-Haïmoud, T. L. Smith, and C. M. F. Mingarelli, Phys. Rev. D **103**, 042009 (2021).

[133] N. M. J. Cruz, A. Malhotra, G. Tasinato, and I. Zavala, arXiv:2402.17312.
[134] E. V. Linder, Phys. Rev. D **34**, 1759 (1986).

[135] C. R. Gwinn, T. M. Eubanks, T. Pyne, M. Birkinshaw, and D. N. Matsakis, Astrophys. J. **485**, 87 (1997).

[136] J. Garcia-Bellido, H. Murayama, and G. White, J. Cosmol. Astropart. Phys. **12** (2021) 023.

[137] S. Jaraba, J. García-Bellido, S. Kuroyanagi, S. Ferraiuolo, and M. Braglia, Mon. Not. R. Astron. Soc. **524**, 3609 (2023).

[138] M. Çalışkan, Y. Chen, L. Dai, N. Anil Kumar, I. Stomberg, and X. Xue, J. Cosmol. Astropart. Phys. **05** (2024) 030.