Gravitational Wave Background search with Pulsar Timing Arrays

Tutorial given at the IPTA Student Week 2026 by Aurélien Chalumeau

Big thanks to Daniel Reardon, Serena Valtolina and Valentina DiMarco for their help !

Table of Content:

• 1. Introduction: The GWB in the context of PTAs
• 2. The simulated IPTA data sets
• 3. Modeling the noise and GWB
• 4. Bayesian Recovery of the GWB
• 5. Optimal Statistics (OS)
• 6. Bonus: Sensitivity as a Function of Pulsar Number

# Use this to allow tab-completion
%config Completer.use_jedi = False
%config IPCompleter.omit__names = 0

from __future__ import division

print('\nLoading modules ...', end='')

# System
import shutil, os, glob, sys
import multiprocessing as mp

# Maths
import numpy as np
#import numpy.ma as ma
import scipy.linalg as sl
from scipy import stats
from scipy.integrate import quad
from scipy.stats import norm
from scipy.optimize import curve_fit
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

# Plots
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
import matplotlib.colors as mcolors
from matplotlib.colors import LogNorm
import corner

# Astro general
from astropy.coordinates import SkyCoord
import astropy.units as u

# PTA specifics

## libstempo
import libstempo as LT

## enterprise & enterprise_extensions
import enterprise
from enterprise.pulsar import Pulsar
from enterprise.signals import parameter, gp_priors, gp_signals, white_signals, signal_base
from enterprise.signals.gp_bases import createfourierdesignmatrix_dm
from enterprise_extensions.sampler import JumpProposal

from enterprise_extensions.blocks import white_noise_block, red_noise_block, dm_noise_block, common_red_noise_block
from enterprise_extensions import hypermodel
from enterprise_extensions.sampler import save_runtime_info
from enterprise_extensions.frequentist import optimal_statistic as ostat

# la_forge
from la_forge.core import Core

# defiant
import defiant
from defiant import OptimalStatistic
from defiant import utils, orf_functions
from defiant import plotting as defplot
from defiant.null_distribution import phase_shift_OS, sky_scramble_OS

# PTMCMC sampler
from PTMCMCSampler.PTMCMCSampler import PTSampler as ptmcmc

print("OK !")

# Matplotlib settings
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'STIXGeneral'
linestyles = "-"
linewidths = 1.0
shade = True
shade_gradient = 2.0
shade_alpha = 0.2
colors = mcolors.CSS4_COLORS['mediumblue']
sigmas = [1, 2]
tfs = 18 #tick
lfs = 27 #label

1

1. Introduction: The Gravitational Wave Background in the context of Pulsar Timing Arrays

Useful functions

def get_HD_curve(zeta):
    coszeta = np.cos(zeta * np.pi / 180)
    x = np.where(zeta == 0, 1., (1 - coszeta) / 2)  # avoid log(0)
    HD = 3/2 * (1/3 + x * (np.log(x) - 1/6))
    return np.where(zeta == 0, 1., HD)

1.1 What is the Gravitational Wave Background (GWB)?

A GWB is a stochastic superposition of gravitational wave signals from a large number of unresolved sources. Pulsar Timing Arrays (PTAs) are sensitive to the GWB in the nanohertz band ($10^{-9}$ – $10^{-6}$ Hz), where the main expected sources are:

1. The population of supermassive black hole binaries (SMBHBs)
2. Topological defects:
   1. Cosmic strings: 1D topological defects from phase transitions in the early Universe
   2. Domain walls: 2D topological defects from discrete symmetry breaking
3. Cosmological first order phase transitions, PTAs are sensitive at the QCD scale
4. Primordial GWBs:
   1. Inflationary GWs. from the amplification of quantum vacuum fluctuations of the gravitational field during inflation
   2. Scalar-induced GWs, second-order effect, arising from enhanced primordial scalar perturbations at second order in perturbation theory

PTA data sets might contain all (or a sub-set) of these astrophysical (for the SMBHB) and cosmological (for the others) signals, thus separating and constraining them will likely become one of the main challenges.

1.2 More about the GW spectrum from SMBHBs

The GW signal from a population of SMBHBs is made of a stochastic background (the incoherent superposition of unresolved binaries) and bright single sources (individually resolvable binaries above the background).

The figure below shows a simulated GW spectrum (orange line) made of thousands of single sources (small black points), while few of them could rise over the global population at their emitted GW frequency (thick black points). The black solid line displays a power-law that might fit well with such a GWB, and the red dotted line shows a typical PTA sensitivity level.

1.3. Standard approximations

The GWB searched by PTAs is often modeled under three standard assumptions:

Assumption	Meaning
Isotropic	The signal power is uniformly distributed across the sky
Stationary	Statistical properties do not evolve over the observation time
Unpolarized	No preferred polarization state; equal power in $+$ and $\times$ modes

We will follow these assumptions in this tutorial.

Key consequence 1:

Under these three approximations, the spatial correlations between any two pulsars i and j depend only on their angular separation $\zeta_{ij}$ and is defined by the famous Hellings-Downs relation $\Gamma(\zeta_{ij})$, which describes the level of correlations within the pulses arrival times between two pulsars, as function of their angular separation, as

$$\Gamma(\zeta_{ij}) = \frac{3}{2} \left[ \frac{1}{3} + x_{ij}\left(\ln x_{ij} - \frac{1}{6}\right) \right] + \frac{1}{2}\delta_{ij}, \quad x_{ij} = \frac{1 - \cos\zeta_{ij}}{2},$$

Here, we use a normalisation so that $\Gamma(0) = 0.5$.

Here is how it actually looks like

zeta = np.linspace(1e-5,180,1000)
HD = get_HD_curve(zeta)

plt.figure(figsize=(10,5))
plt.plot(zeta, HD, lw=3)
plt.axhline(0, c='k', ls='--', zorder=0, alpha=.3)
plt.xlabel("Angle $\zeta_{ij}$ between Earth-pulsar baseline pairs [deg.]", fontsize=18)
plt.ylabel("Arrival time correlation $\Gamma(\zeta_{ij})$", fontsize=18)
plt.xticks(fontsize=14)
plt.yticks(fontsize=14)
plt.tight_layout()
plt.grid(alpha=.3)
plt.show()

The Hellings-Downs (HD) correlation is the smoking gun of the GWB: detecting it in PTA data is what distinguishes a true GW signal from any other noise or signal.

Key consequence 2:

The GWB is fully characterized by its one-sided power spectral density: $S_h(f)$

We typically model a power-law and mention the predicted gammas

$$S_h(f) = \frac{A^2}{12\pi^2} \left(\frac{f}{f_{\rm ref}}\right)^{-\gamma} f^{-3} \quad \text{[s}^3\text{]}$$

where $A$ is the GWB characteristic strain amplitude at $f_{\rm ref}$, $\gamma$ the spectral index, which is equal to $13/3$ for a population of circular and GW-driven SMBHBs, and $f_{\rm ref}$ is typically $1\,\text{yr}^{-1}$.

With this tutorial, we will

• Vizualize a simulated PTA data set made of 133 pulsars
• See how we could model the GWB signal in PTA data sets
• See how to run a Bayesian analysis to (1) constrain the GWB spectral properties.
• See how to run the Optimal Statistic approach (OS) to (1) recover the GWB amplitude, (2) recover the HD curve and (3) compute the "frequentist" significance of the signal
• A Bonus part will follow right after that

Some selected references related to Section 1

• The Nanohertz Gravitational Wave Astronomer (Taylor, 2021), very useful book about PTA.
• Upper limits on the isotropic gravitational radiation background from pulsar timing analysis (Hellings & Downs, 1983), reference for the Hellings & Downs correlations
• A Practical Theorem on Gravitational Wave Backgrounds (Phinney, 2001), reference for the predicted power-law slope for a GW-driven and circular population of supermassive black hole binaries

2

2. The simulated IPTA data sets (2 versions here)

</h1>

Useful functions

#############
### Plot Pulsar positions in a sky map projection
#############

def plot_SkyMap(psrs, coordtype, psrcolor='w', show_mw=True):
    # Initiate and set Figure
    fig = plt.figure(figsize=(12, 6))
    fig.patch.set_facecolor('k')
    ax = fig.add_subplot(111, projection="mollweide") # projections could be "mollweide", "hammer", "aitoff" or "lambert"
    ax.set_facecolor('k')
    
    # Milky Way disk
    if show_mw:
        Plot_MW(ax, coordtype=coordtype, n_stars=5e5, angle_mask=10)

    # Pulsars
    for psr in psrs:
        if coordtype=="gal":
            coord = SkyCoord(ra=psr._raj * u.rad, dec=psr._decj * u.rad, frame='icrs')
            gal = coord.galactic
            x = gal.l.rad
            if x > np.pi:
                x -= 2 * np.pi
            y = gal.b.rad
        elif coordtype=="eq":
            x  = psr._raj - np.pi   # shift [0,2π] → [-π,π]
            y = psr._decj
        ax.plot(x, y, marker='*', markerfacecolor=psrcolor, markeredgecolor=psrcolor, markeredgewidth=0.4, markersize=8, zorder=10)

    # Set Legends
    legend_elements = [
        Line2D([0], [0], marker='*', color=psrcolor, label='Pulsar', markerfacecolor='w', markersize=8, lw=0),
        Patch(facecolor='darkorange', edgecolor='k', lw=1, linestyle=':', label='Milky Way')]
    leg = plt.legend(handles=legend_elements, fontsize=16, loc='upper right', facecolor='k', edgecolor='white', labelcolor='white')

    # Set Axes
    xticks = np.linspace(-np.pi, np.pi, 13)
    if coordtype=="gal":
        xtick_labels = [f"{int(lon)}°" for lon in np.linspace(180, -180, 13, endpoint=True)]
        ax.set_xlabel('Galactic Longitude $l$', fontsize=20, color='white')
        ax.set_ylabel('Galactic Latitude $b$', fontsize=20, color='white')
    elif coordtype=="eq":
        xtick_labels = [f"{int(h)}h" for h in np.linspace(0, 24, 13, endpoint=True)]
        ax.set_xlabel('RA [hr]', fontsize=18, color='white')
        ax.set_ylabel('DEC [deg]', fontsize=18, color='white')
    ax.set_xticks(xticks)
    ax.set_xticklabels(xtick_labels, fontsize=20, color='white')
    ax.tick_params(colors='white', labelsize=20)

    # Set Mollweide frame/grid
    ax.grid(alpha=0.4, color='white', zorder=1)
    for spine in ax.spines.values():
        spine.set_edgecolor('white')

    plt.tight_layout()
    plt.show()

def Plot_MW(ax, coordtype="eq", n_stars=500000, angle_mask=4):
    rng = np.random.default_rng(42)
    b_thin  = rng.normal(0, np.deg2rad(3),  int(n_stars * 0.80))
    b_thick = rng.normal(0, np.deg2rad(15), int(n_stars * 0.15))
    b_halo  = rng.normal(0, np.deg2rad(40), int(n_stars * 0.05))
    b_fake  = np.clip(np.concatenate([b_thin, b_thick, b_halo]), -np.pi/2, np.pi/2)
    
    # --- Longitude with galactic-center concentration ---
    # Sample l from a mixture: a narrow Gaussian at l=0 (bulge)
    # + a broad Gaussian (disk) + a uniform floor (halo/background)
    n = len(b_fake)
    l_bulge  = rng.normal(0, np.deg2rad(15),  int(n * 0.25))   # central bulge
    l_disk   = rng.normal(0, np.deg2rad(100),  int(n * 0.55))   # inner disk
    l_uniform= rng.uniform(-np.pi, np.pi,     int(n * 0.20))   # outer disk / background
    l_fake   = np.concatenate([l_bulge, l_disk, l_uniform])
    # Wrap to [-π, π]
    l_fake   = (l_fake + np.pi) % (2 * np.pi) - np.pi
    l_fake   = l_fake[:n]

    n_bins_x, n_bins_y = 720, 360

    if coordtype == "eq":
#         l_fake = rng.uniform(0, 2*np.pi, int(n_stars))[:len(b_fake)]

        coords  = SkyCoord(l=l_fake*u.rad, b=b_fake*u.rad, frame='galactic').icrs
        x_fake  = np.clip(coords.ra.rad - np.pi, -np.pi, np.pi)  # clip after shift
        y_fake  = coords.dec.rad

        h, xe, ye = np.histogram2d(x_fake, y_fake, bins=[n_bins_x, n_bins_y],
                                   range=[[-np.pi, np.pi], [-np.pi/2, np.pi/2]])
        xc = 0.5 * (xe[:-1] + xe[1:])
        yc = 0.5 * (ye[:-1] + ye[1:])
        X, Y = np.meshgrid(xc, yc)

    elif coordtype == "gal":
#         l_fake = rng.uniform(-np.pi, np.pi, int(n_stars))[:len(b_fake)]

        h, xe, ye = np.histogram2d(l_fake, b_fake, bins=[n_bins_x, n_bins_y],
                                   range=[[-np.pi, np.pi], [-np.pi/2, np.pi/2]])
        xc = 0.5 * (xe[:-1] + xe[1:])
        yc = 0.5 * (ye[:-1] + ye[1:])
        X, Y = np.meshgrid(xc, yc)

    # Same normalization logic for both branches
    h_masked = np.ma.masked_where(h < angle_mask, h)
    vmax = np.percentile(h[h > 0], 99.5)

    ax.pcolormesh(X, Y, h_masked.T,
                  norm=LogNorm(vmin=1, vmax=vmax),
                  cmap='inferno', shading='auto',
                  zorder=0, rasterized=True, alpha=0.6)

#############
### Plot HD correlations
#############
    
def get_zetas(psrs):
    zetas = []
    for i, psr1 in enumerate(psrs):
        for j, psr2 in enumerate(psrs):
            if j <= i:
                continue
            angle = np.arccos(np.clip(np.dot(psr1.pos, psr2.pos), -1, 1)) * 180 / np.pi
            zetas.append(angle)
    return np.array(zetas)
    
def plot_HD_from_psrlist(psrs, bin_size_deg):
    # Get pulsar pair angles
    zetas = get_zetas(psrs)
    
    plt.figure(figsize=(8,5))
    ax = plt.gca()
    plt.suptitle(f"{len(psrs)} pulsars - {len(zetas)} pairs", fontsize=20)

    HD = get_HD_curve(zetas)
    ax.plot(zetas, HD, '.', c='k', ms=5, alpha=1)
    ax.axhline(0, c='k', zorder=0, alpha=.4)
    ax.axvline(0, c='k', ls='--', alpha=.3)
    ax.axvline(180, c='k', ls='--', alpha=.3)

    ax.set_xlabel("Angle $\zeta_{ij}$ between Earth-pulsar baseline pairs [deg.]", fontsize=18)
    ax.set_ylabel("Arrival time correlation $\Gamma(\zeta_{ij})$", fontsize=18)
    ax.tick_params(axis='both', which='major', labelsize=14)
    ax.grid(alpha=.3, zorder=0)

    ax2 = ax.twinx()
    bin_edges   = np.arange(0, 180 + bin_size_deg, bin_size_deg)  # [0, 5, 10, ..., 180]
    bin_centers = bin_edges[:-1] + bin_size_deg / 2               # [2.5, 7.5, ..., 177.5]
    bin_counts, _ = np.histogram(zetas, bins=bin_edges)  
    ax2.bar(bin_centers, bin_counts, edgecolor='black', width=bin_size_deg, alpha=0.3, color='cornflowerblue', zorder=1, label='Number of pairs per bin')

    ax2.set_ylabel('Histogram: number of pairs', fontsize=16)
    ax2.set_ylim(0, 2*np.max(bin_counts))
    ax2.tick_params(axis='both', which='major', labelsize=14)

    plt.tight_layout()
    plt.show()
    
#############
### Plot timing residuals
#############

def mjd2greg(mjd):
    return 2000 + (np.array(mjd)-51544.5)/365.25

def greg2mjd(greg):
    return (np.array(greg) - 2000) * 365.25 + 51544.5
    
def plot_pulsar_timing(psr, prefit_res=None, plot_histo=False, color=None, colorname=None, addtitle=""):
    
    # --- Layout logic ---
    nrows = 2 if prefit_res is not None else 1
    ncols = 2 if plot_histo else 1

    fig, axes = plt.subplots(
        nrows, ncols,
        figsize=(15 if plot_histo else 13, 5 if nrows == 1 else 4),
        gridspec_kw={'width_ratios': [6, 1] if plot_histo else [1],
                     'hspace': 0., 'wspace': 0.}
    )

    # Normalize axes array shape
    if nrows == 1:
        axes = np.array([axes])
    if ncols == 1:
        axes = axes.reshape(nrows, 1)

    if prefit_res is None:
        y = 1.05
    else:
        y = 1.09
        
        
    fontsize = 14
    
    fig.suptitle(f"{psr.name}{addtitle}", fontsize=fontsize+3, y=y)

    # --- Color handling ---
    if color is None:
        color = psr.freqs
        colorname = "Frequency [$MHz$]"

    cmap = mcolors.LinearSegmentedColormap.from_list(
        'coral_blue', ['coral', 'cornflowerblue']
    )
    norm = mcolors.Normalize(vmin=np.min(color), vmax=np.max(color))
    
    

    def plot_residual_panel(ax, residuals):
        for toa, res, err, cval in zip(psr.toas, residuals, psr.toaerrs, color):
            ax.errorbar(
                toa/86400,
                res*1e6,
                yerr=err*1e6,
                fmt='.',
                ms=6,
                capsize=2,
                color=cmap(norm(cval))
            )

        ax.axhline(0., ls=':', c='k', lw=2)
        ax.grid(alpha=0.2)

    # Prefit (top)
    if prefit_res is not None:
        ax_prefit = axes[0, 0]
        plot_residual_panel(ax_prefit, prefit_res)
        ax_prefit.set_ylabel('Prefit residuals\n[$\mu s$]', fontsize=fontsize)
        ax_prefit.tick_params(axis='both', which='major', labelsize=fontsize)
        plt.setp(ax_prefit.get_xticklabels(), visible=False)
        
        # Upper x-axis with Gregorian dates
        ax_top = ax_prefit.twiny()
        ax_top.set_xlim(ax_prefit.get_xlim())
        # Pick ~6 evenly spaced tick positions in MJD
        display_cadence = 2 # years
        greg_labels = [int(i) for i in np.arange(int(mjd2greg(psr.toas.min()/86400)), int(mjd2greg(psr.toas.max()/86400))+1, display_cadence)]
        upperticks_mjd = greg2mjd(greg_labels)
        ax_top.set_xticks(upperticks_mjd)
        ax_top.set_xticklabels(greg_labels, fontsize=fontsize)
        ax_top.set_xlabel('Date', fontsize=fontsize)
        ymin, ymax = ax_top.get_ylim()
        ylim = max(abs(ymin), abs(ymax))
        ax_top.set_ylim(-ylim, ylim)
        ax_prefit.set_xticks(upperticks_mjd, minor=True)
        ax_prefit.grid(which='minor', axis='x', alpha=1, lw=1, ls='--')

    # Postfit (bottom or only)
    ax_post = axes[-1, 0]
    plot_residual_panel(ax_post, psr.residuals)
    ax_post.set_xlabel('Epoch [MJD]', fontsize=fontsize)
    ax_post.set_ylabel('Timing residuals\n[$\mu s$]', fontsize=fontsize)
    ax_post.tick_params(axis='both', which='major', labelsize=fontsize)
    ymin, ymax = ax_post.get_ylim()
    ylim = max(abs(ymin), abs(ymax))
    ax_post.set_ylim(-ylim, ylim)
    
    if prefit_res is None:
        # Upper x-axis with Gregorian dates
        ax_top = ax_post.twiny()
        ax_top.set_xlim(ax_post.get_xlim())
        # Pick ~6 evenly spaced tick positions in MJD
        display_cadence = 2 # years
        greg_labels = [int(i) for i in np.arange(int(mjd2greg(psr.toas.min()/86400)), int(mjd2greg(psr.toas.max()/86400))+1, display_cadence)]
        upperticks_mjd = greg2mjd(greg_labels)
        ax_top.set_xticks(upperticks_mjd)
        ax_top.set_xticklabels(greg_labels, fontsize=fontsize)
        ax_top.set_xlabel('Date', fontsize=fontsize)
    
    ax_post.set_xticks(upperticks_mjd, minor=True)
    ax_post.grid(which='minor', axis='x', alpha=1, lw=1, ls='--')
    
    # --- Colorbar ---
    cax = ax_post.inset_axes([0.02, 0.01, 0.4, 0.03])
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    sm.set_array([])
    cb = fig.colorbar(sm, cax=cax, orientation='horizontal')
    cax.xaxis.set_label_position('top')
    cax.xaxis.tick_top()
    
    if colorname is None:
        cb.set_label('Color scale', size=fontsize, labelpad=5)
    else:
        cb.set_label('Frequency [MHz]', size=fontsize, labelpad=5)
    cb.ax.tick_params(labelsize=fontsize)
    

    # --- Histogram(s) ---
    if plot_histo:
        def plot_hist(ax, residuals):
            wr = residuals / psr.toaerrs
            counts, bins = np.histogram(wr, bins=30)
            counts = counts / counts.max()
            ax.stairs(counts, bins, orientation="horizontal", lw=2)

            # Gaussian reference
            g = np.random.randn(10000)
            c2, b2 = np.histogram(g, bins=30)
            c2 = c2 / c2.max()
            ax.stairs(c2, b2, orientation="horizontal", lw=2, color='green', label="$N(0,1)$")
            ax.legend(fontsize=fontsize)

            ax.grid(alpha=0.4)
            plt.setp(ax.get_xticklabels(), visible=False)
            plt.setp(ax.get_yticklabels(), visible=False)
            ymin, ymax = ax.get_ylim()
            ylim = max(abs(ymin), abs(ymax))
            ax.set_ylim(-ylim, ylim)
            ax.yaxis.set_label_position("right")
            ax.yaxis.tick_right()
            ax.set_ylabel("Weighted residuals", fontsize=fontsize, rotation=-90, labelpad=20)
            # Remove y ticks AND tick labels
            ax.set_yticks([])
            ax.tick_params(axis='y', which='both', length=0)

        if prefit_res is not None:
            plot_hist(axes[0, 1], prefit_res)

        plot_hist(axes[-1, 1], psr.residuals)

    plt.show()

2.1 Main properties

We work with a simulated IPTA data set, for which par/tim files are located here

./DataSim/ideal/PTAdata/pulsarname/pulsarname.par
./DataSim/ideal/PTAdata/pulsarname/pulsarname.tim

Dataset: "ideal"

Main properties of the data set

• Pulsars: all the same isolated pulsars, but positions placed at real PTA pulsar positions
• Epoch properties: $30$ days cadence from MJD 55000 to MJD 61000, thus spanning $6000$ days (~$16$ years),
• Errorbars randomly sampled from $0.05\,\mu$s to $1\,\mu$s
• 3 observing radio frequencies among 400, 1200 and 1400 MHz, not important here as we do not consider chromatic (i.e., radio-frequency dependent) noise

What the timing residuals contain

• White noise: EFAC only
• GWB: An Isotropic, unpolarized, stationary and Gaussian Common Red Signal with Hellings-Downs spatial correlations, with amplitude $A = 10^{-15}$ and $\gamma = 13/3 \approx 4.33$ (predicted from a GW-driven population of circular SMBHBs).

</div>

Let us finally start to play with the notebook, by reading and vizualizing the chosen data set.

2.2 Reading the parfiles & timfiles with libstempo (a python wrapper for Tempo2)

currentdir = os.getcwd()
print('Current directory: '+currentdir)

# Choose the data set
dataset = "ideal_for_IPTASW"

# Define data dir
datadir = f'{currentdir}/DataSim/{dataset}/'

# Find all par/tim directories
partimdirs = np.sort(glob.glob(f"{datadir}/PTAdata/*"))

# Define pulsar list
psrnames = [os.path.basename(d).split("/")[0] for d in partimdirs]

# Define lists of parfiles and timfiles
parfiles = [f"{d}/{psrnames[i]}.par" for i,d in enumerate(partimdirs)]
timfiles = [f"{d}/{psrnames[i]}.tim" for i,d in enumerate(partimdirs)]

# Read the data with libstempo and fit for the timing model
psrs = []
prefit_res = {}
for parfile, timfile in zip(parfiles, timfiles):
    ltpsr = LT.tempopulsar(parfile, timfile)
    prefit_res.update({ltpsr.name:np.copy(ltpsr.residuals())})
    ltpsr.fit()
    psr = Pulsar(ltpsr)
    psrs.append(psr)
    print(ltpsr.name)

2.3 Sky distribution of the pulsars

The millisecond pulsars (MSPs) are not uniformly distributed on the sky, they are predominantly found along the Galactic plane, reflecting both theit intrinsic distribution in the Milky Way and observational selection effects (sensitivity of radio telescopes, scattering from the ionized interstellar medium).

# Equatorial coordinates: "eq" ; Galactic coordinates: "gal"
coordtype = "eq"

plot_SkyMap(psrs, coordtype, show_mw=True, psrcolor = "w")

2.4 Pulsar pair distribution of the data set

To constrain the Arrival time correlations and search for the HD signature, we need to have a lot of pulsar pairs, so a lot of pulsars. Let us plot our coverage for the HD curve.

FYI, the number of pairs for $n$ pulsars is $n(n-1)/2$

# Play with pulsar numbers ; Max is 133 here (our full data set)
Npsrs = 133

# Reducing the number of pulsars to play with the curve
selpsrs = np.random.choice(psrs, size=Npsrs, replace=False)

# Bin size for the histogram (in deg.)
bin_size_deg=5

# Plot !
plot_HD_from_psrlist(psrs=selpsrs, bin_size_deg=bin_size_deg)

2.5 Prefit vs. Postfit residuals

Let us now plot the timing residuals (before and after fitting for the timing model parameters)

for psr in psrs[-35:-30]:
    plot_pulsar_timing(psr, prefit_res=prefit_res[psr.name], plot_histo=False, addtitle=f" - Data set: {dataset}")
#     plot_pulsar_timing(psr, plot_histo=True, addtitle=f" - Data set: {dataset}")

Some selected references related to Section 2

• Tempo2, a new pulsar timing package. II: The timing model and precision estimates (Edwards et al., 2006), The Tempo2 timing model
• libstempo: Python wrapper for Tempo2 (Vallisneri et al., 2020), reference for Libstempo

3

3. Modeling the noise and GWB

</h1>

Useful functions

def set_pta_enterprise(psrs, gwb_model, noise_model="ideal", Npsrs=None, verbose=False):
    
    # Choose a random subset of Npsrs pulsars
    if Npsrs is None:
        Npsrs = len(psrs)
    selpsrs = np.random.choice(psrs, size=Npsrs, replace=False)
    
    # Sort pulsar objects from names 
    idxs = np.argsort([p.name for p in selpsrs])
    selpsrs = [psr for psr in selpsrs[idxs]]
    
    ###########################
    ## Set up the noise model
    ###########################
    
    ########################### Timing Model marginalization
    # Initiate the enterprise signal_collection object with the TimingModel
    tm = gp_signals.TimingModel()
    
    ########################### White noise
    ## EFAC - Here we don't consider EQUAD/ECORR, just for simplicity
    efac = white_noise_block(vary=False, tnequad=None, select=None)
    
    ########################### Achromatic Red Noise - Not used here
    rn = red_noise_block(psd="powerlaw", components=30, prior="log-uniform", name="red_noise")
    
    ########################### DM variations - Not used here
    dmgp = dm_noise_block(psd="powerlaw", components=30, prior="log-uniform", name="dm_gp")
    
    ########################### Common Red Signal
    # Common Red Signal uncorrelated power-law
    if gwb_model == "curn_pl":
        orf = None
        crn_psd = "powerlaw"
        crn_name = "gw_curn_pl"
    
    # Common Red Signal uncorrelated free-spectrum
    elif gwb_model == "curn_fs":
        orf = None
        crn_psd = "spectrum"
        crn_name = "gw_curn_fs"
    
    # Helling-Downs correlated power-law
    elif gwb_model == "hd_pl":
        orf = "hd"
        crn_psd = "powerlaw"
        crn_name = "gw_hd_pl"
    
    crn = common_red_noise_block(psd=crn_psd, components=30, prior='log-uniform', orf=orf, name=crn_name)
    
    ###########################
    ## Set up the enterprise SignalCollection object
    ###########################
    if noise_model=="ideal":
        signal = tm + efac + crn
    elif noise_model=="realistic":
        signal = tm + efac + rn + dmgp + crn
    
    model = [signal(psr) for psr in selpsrs]
    
    ###########################
    ## Set up the enterprise PTA object
    ###########################
    pta = signal_base.PTA(model)
    
    ########################### Fix EFAC values to 1.
    params = {}
    for psrname in psrnames:
        parname = f"{psrname}_efac"
        parval = 1.
        params.update({parname:parval})
        
    pta.set_default_params(params)
    
    if verbose:
        print(f"PTA object set for {Npsrs} pulsars, using '{noise_model}' noise model and '{gwb_model}' CRS model.")
    
    return selpsrs, pta

3.1 The PTA Likelihood

In order to constrain models in PTA data, we build a likelihood function: the probability of observing the data given a chosen model and its parameter values.

The full PTA likelihood is a multivariate Gaussian over all pulsars jointly: $$\mathcal{L}(\delta\mathbf{t} \mid \boldsymbol{\theta}) = \frac{1}{\sqrt{(2\pi)^{n} \det \mathbf{C}}} \exp\left(-\frac{1}{2} \delta\mathbf{t}^\top \mathbf{C}^{-1} \delta\mathbf{t}\right)$$ where $\delta\mathbf{t} = (\delta t_1, \ldots, \delta t_{N_p})$ is the concatenated residual vector over all pulsars, and $\mathbf{C}$ is the full covariance matrix.

In this tutorial, we consider three components: the timing model solution errors, the white noise (EFAC only) and the Gravitational Wave Background, all included in the covariance matrix as:

$$\mathbf{C} = \mathbf{C}^{\rm TM} + \mathbf{C}^{\rm WN} + \mathbf{C}^{\rm GWB}$$

Components of the covariance matrix

1. White noise $\mathbf{C}^{\rm WN}$

Each TOA has an independent measurement uncertainty $\sigma_{\rm TOA}$, so the white noise covariance is purely diagonal: $$\mathbf{C}^{\rm WN}_{ij} = {\rm EFAC} \ \sigma_{\rm TOA}^2(t_i) \ \delta_{ij} = \sigma_{\rm TOA}^2(t_i) \ \delta_{ij}$$ Here EFAC is fixed to 1 (its injected value), so it does not rescale the uncertainties. No free parameters.

2. Timing model $\mathbf{C}^{\rm TM}$

Rather than fitting for timing model parameters explicitly, we marginalize over them analytically. See references as Haasteren & Vallisneri (2014) for more details. No free parameters.

3. Gravitational Wave Background $\mathbf{C}^{\rm GWB}$

The GWB is modeled as a stationary Gaussian process with Fourier basis functions $\Phi$ at frequencies indexed with $\mu$ and $\nu$. Its covariance between TOAs observed at times $t_{i, a}$ and $t_{j, b}$, respectively for pulsars $a$ and $b$ is:
$$\mathbf{C}^{\rm GWB}_{ab}(t_{i,a}, t_{jb}) = \Gamma_{ab} \cdot \sum_{\mu, \nu} \Phi_{\mu}(t_{i,a}) \Sigma_{\mu \nu} \Phi_{\nu}(t_{j,b})$$
where $\Sigma_{\mu \nu} = S(f_{\mu}) \cdot \delta_{\mu \nu} / T$, with $S(f)$ the GWB power spectral density (PSD) at Fourier frequency $f$, $\delta_{\mu \nu}$ the kronecker delta and $T$ the total observing time span of the data set. $\Gamma_{ab}$ is the overlap reduction function (ORF), which encodes the spatial correlation between pulsar pairs. The ORF depends on the angular separation $\zeta_{ab}$ on the sky; for an isotropic GWB it takes the form of the Hellings-Downs curve (HD). Setting $\Gamma_{ab} = \delta_{ab}$ instead gives a common uncorrelated red noise (CURN), a common red process with no spatial correlations.
For the PSD model, the two choices considered in this tutorials are:

• Power law parametrized by an amplitude $A$ and spectral index $\gamma$: $$S(f) = \frac{A^2}{12\pi^2} \left(\frac{f}{f_{\rm yr}}\right)^{-\gamma} f_{\rm yr}^{-3}$$ For a GW-driven population of circular SMBHBs, $\gamma = 13/3$.
  The free parameters are $\{A, \gamma\}$.
• Free spectrum, parametrized by the PSD amplitude $S(f_k)$ at each Fourier frequency bin $f_k$, allowing a non-parametric reconstruction of the GWB spectrum. The free parameters are $\{S(f_1), \ldots, S(f_{N_f})\}$.

In this tutorial, we consider a CURN ($\Gamma_{ab} = \delta_{ab}$), eith with a power-law or a free spectrum model for the Bayesian analysis. The HD model will be applied for the Optimal Statistics section.

3.2 How to compute it

Let us now compute the PTA likelihood (and priors for the Bayesian analysis) with Enterprise "by hand".

# Initiate the enterprise signal_collection object with the TimingModel
tm = gp_signals.TimingModel(use_svd=True)

########################### White noise
## EFAC - Here we don't consider EQUAD/ECORR, just for simplicity
efac_prior = parameter.Constant()
efac = white_signals.MeasurementNoise(efac=efac_prior)

########################### Common Red Signal
# Common Red Signal uncorrelated power-law
orf = None # "hd"
crn_psd = "powerlaw" # "spectrum"
crn_name = "gw_curn_pl"

crn = common_red_noise_block(psd=crn_psd, components=30, prior='log-uniform', orf=orf, name=crn_name)

###########################
## Set up the enterprise SignalCollection object
###########################
signal = tm + efac + crn
model = [signal(psr) for psr in psrs]

###########################
## Set up the enterprise PTA object
###########################
pta = signal_base.PTA(model)

########################### Fix EFAC values to 1.
params = {}
for psrname in psrnames:
    parname = f"{psrname}_efac"
    parval = 1.
    params.update({parname:parval})

pta.set_default_params(params)


print(f"PTA object set for {len(psrs)} pulsars, using {str(orf)} as ORF and '{crn_psd}' as PSD.")

3.3 Playing with the PTA object

y = np.array([4.33, -15])
cov = np.diag([.1, .01])

N = 3000
    
xs = np.zeros((3, N))
for i in range(N):
    print(f"{(i+1)} on {N}", end="\r")
    g, A = np.random.multivariate_normal(y, cov=cov)

    x = {
        'gw_curn_pl_gamma':g, 
        'gw_curn_pl_log10_A':A
    }

    xs[0,i] = g
    xs[1,i] = A
    xs[2,i] = pta.get_lnlikelihood(x)

fig, ax = plt.subplots(figsize=(8, 6))
sc = ax.scatter(xs[0], xs[1], c=xs[2]-xs[2].max(), vmin=-20, vmax=0, cmap='viridis', s=10)
cb = plt.colorbar(sc, ax=ax)
cb.set_label(label=r'$\ln\mathcal{L} - \ln\mathcal{L}_{\rm max}$ (from -20 to 0)', size=16, labelpad=5)

# Mark injected values
ax.axvline(4.33, ls='--', color='k', lw=1.5, label='Injected values', alpha=.8)
ax.axhline(-15,  ls='--', color='k', lw=1.5, alpha=.8)
ax.set_xlabel(r'$\gamma$', fontsize=16)
ax.set_ylabel(r'$\log_{10} A$', fontsize=18)
ax.set_title(r'$\ln\mathcal{L}$ map around the injected value', fontsize=16)
ax.legend(fontsize=16)
plt.tight_layout()
plt.show()

pta.summary()

Let us now perform a Bayesian analysis to constrain the GWB signal We will then use the function "set_pta_enterprise" to compute the PTA object

Some selected references related to Section 3

• The Nanohertz Gravitational Wave Astronomer (Taylor, 2021), very useful book about PTA, again related to this part.
• New advances in the Gaussian-process approach to pulsar-timing data analysis (van Haasteren & Vallisneri, 2014), reference to see how the PTA likelihood is built and how the stochastic signals are modeled as Gaussian processes.
• ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE (Ellis, Siemens & Hazboun, 2019), reference for Enterprise software

4

4. Bayesian analysis

</h1>

Useful functions

# Module-level globals, populated before forking, inherited by workers
_GLOBAL_PTA = None
_GLOBAL_RUN_KWARGS = None

def _single_PTMCMC_run(outdir):
    """Worker function: reads pta and kwargs from globals, no pickling needed."""
    pta = _GLOBAL_PTA
    nsamples    = _GLOBAL_RUN_KWARGS["nsamples"]
    SCAMweight  = _GLOBAL_RUN_KWARGS["SCAMweight"]
    AMweight    = _GLOBAL_RUN_KWARGS["AMweight"]
    DEweight    = _GLOBAL_RUN_KWARGS["DEweight"]
    overwrite_dir = _GLOBAL_RUN_KWARGS["overwrite_dir"]

    if os.path.exists(outdir) and not overwrite_dir:
        print(f"\nOutdir already exists: {outdir}")
        return

    if os.path.exists(outdir) and overwrite_dir:
        shutil.rmtree(outdir)
        os.makedirs(outdir)
        print(f"\nCleared and recreated: {outdir}")
    elif not os.path.exists(outdir):
        os.makedirs(outdir)
        print(f"\nCreated folder: {outdir}")

    x0 = np.hstack([p.sample() for p in pta.params])
    ndim = len(x0)
    cov = np.diag(np.ones(ndim) * 0.01**2)
    sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, outDir=outdir)
    save_runtime_info(pta, outdir=outdir)

    jp = JumpProposal(pta, None, empirical_distr=None)
    sampler.addProposalToCycle(jp.draw_from_prior, 5)
    sel_sig = {"red_noise": 10, "dm": 10, "gw": 30}
    for sig, val in sel_sig.items():
        if any([sig in p for p in pta.param_names]):
            sampler.addProposalToCycle(jp.draw_from_par_prior(sig), val)

    print(f"\nPTMCMC run started in '{outdir}' with {nsamples} iterations.\n")
    sampler.sample(x0, int(nsamples), SCAMweight=SCAMweight, AMweight=AMweight, DEweight=DEweight)


def run_PTMCMC(outdir, pta, nsamples=1e5, SCAMweight=30, AMweight=15, DEweight=50,
               overwrite_dir=False, NRuns=1):

    global _GLOBAL_PTA, _GLOBAL_RUN_KWARGS  # single declaration at top of function
    _GLOBAL_PTA = pta
    _GLOBAL_RUN_KWARGS = dict(nsamples=nsamples, SCAMweight=SCAMweight,
                              AMweight=AMweight, DEweight=DEweight,
                              overwrite_dir=overwrite_dir)

    if NRuns == 1:
        _single_PTMCMC_run(f"{outdir}")

    else:
        outdirs = [f"{outdir}/ptmcmc_{i}" for i in range(NRuns)]
        print(f"\nLaunching {NRuns} parallel PTMCMC runs...")
        ctx = mp.get_context("fork")
        with ctx.Pool(processes=NRuns) as pool:
            pool.map(_single_PTMCMC_run, outdirs)
        print(f"\nAll {NRuns} runs completed.")
        
def run_PTMCMC_old(outdir, pta, nsamples=1e5, SCAMweight=30, AMweight=15, DEweight=50, overwrite_dir=False):
    
    if os.path.exists(outdir) and not overwrite_dir:
        print(f"\nOutdir already exists: {outdir}")
    else:
        if os.path.exists(outdir) and overwrite_dir:
            shutil.rmtree(outdir)
            os.makedirs(outdir)
            print(f"\nCleared and recreated: {outdir}")
        
        elif not os.path.exists(outdir):
            os.makedirs(outdir)
            print(f"\nCreated folder: {outdir}")
    
        ######################### Generic Set-up part
        # Define the initial point
        x0 = np.hstack([p.sample() for p in pta.params])

        # Initialize the parameter covariance matrix used by SCAM and AM Jump Proposals
        ndim = len(x0)
        cov = np.diag(np.ones(ndim) * 0.01**2)

        # Initialize the PTMCMC objectz
        sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, outDir=outdir)

        # Save set-up information
        save_runtime_info(pta, outdir=outdir)

        ######################### Additionals Jump Proposals
        # Set Jump Proposals
        jp = JumpProposal(pta, None, empirical_distr=None)

        # always add draw from prior
        sampler.addProposalToCycle(jp.draw_from_prior, 5)

        # Jump proposals from priors of selected params
        sel_sig = {
            "red_noise":10, 
            "dm":10, 
            "gw":30
        }
        for sig, val in sel_sig.items():
            if any([sig in p for p in pta.param_names]):
                sampler.addProposalToCycle(jp.draw_from_par_prior(sig), val)

        ######################### Sample !
        print(f"\nPTMCMC run started with {nsamples} iterations.")
        sampler.sample(x0, int(nsamples), SCAMweight=SCAMweight, AMweight=AMweight, DEweight=DEweight)
        
def read_chains(chaindir, burnin=0.3):
    ch = np.loadtxt(chaindir + '/chain_1.txt')
    ch = ch[int(len(ch) * burnin):]
    pars = np.loadtxt(chaindir + '/pars.txt', dtype='str')
    return ch, pars


def read_all_chains(outdir, NRuns=None, burnin=0.3):
    chains, pars = [], None

    if NRuns is None:
        # Single run, no subdirectory
        chains, pars = read_chains(outdir, burnin=burnin)
        return chains, pars

    if NRuns == 'auto':
        # Find all ptmcmc_{i} folders automatically
        import os
        subdirs = sorted([
            d for d in os.listdir(outdir)
            if d.startswith('ptmcmc_') and os.path.isdir(os.path.join(outdir, d))
        ])
        run_indices = [int(d.split('_')[1]) for d in subdirs]
    else:
        run_indices = range(NRuns)

    for i in run_indices:
        chaindir = f"{outdir}/ptmcmc_{i}"
        ch, p = read_chains(chaindir, burnin=burnin)
        chains.append(ch)
        if pars is None:
            pars = p

    chains = np.concatenate(chains, axis=0)
    return chains, pars

def compute_rho(log10_A, gamma, f, T):
    """
    Converts from power to residual RMS.
    """
    fyr = 1 / (365.25*86400)
    return np.sqrt((10**log10_A)**2 / (12.0*np.pi**2)
                   * fyr**(gamma-3) * f**(-gamma) / T)

4.1 Parameter estimation on the CURN power-law PSD model

We infer the GWB parameters $\theta = \{A, \gamma\}$ from the data using Bayes' theorem:

$$p(\theta \mid \delta t) = \frac{p(\delta t \mid \theta)\, p(\theta)}{p(\delta t)}$$

where:

• $p(\delta t \mid \theta)$ is the likelihood: for our Gaussian noise model, it is a multivariate Gaussian in the timing residuals $\delta t$
• $p(\theta)$ is the prior: log-uniform on $A$ (in practice it is uniform on $\log_{10} A$) and uniform on $\gamma$
• $p(\delta t)$ is the Bayesian evidence, used for model comparison via Bayes factors

We sample the posterior using MCMC via the enterprise + PTMCMCSampler packages.

Two GWB models are run here:

• CURN power-law (curn_pl): a Common Uncorrelated Red Noise with a power-law PSD, fast to sample, but ignores inter-pulsar spatial correlations.
• HD power-law (hd_pl): the same power-law PSD but with the Hellings-Downs ORF enforced, this is the true GWB model, but computing the HD covariance matrix at every likelihood evaluation is much more expensive. We won't really use this one.

overwrite_dir = False
gwb_model = "curn_pl"
outdir_curnpl = f"{datadir}/chains/CURN_pl/"

nsamples = 1e5

selpsrs, pta = set_pta_enterprise(psrs, gwb_model, verbose=True)

if os.path.exists(outdir_curnpl) and not overwrite_dir:
    print(f"\nOutdir already exists: {outdir_curnpl}")
else:
    if os.path.exists(outdir_curnpl) and overwrite_dir:
        shutil.rmtree(outdir_curnpl)
        os.makedirs(outdir_curnpl)
        print(f"\nCleared and recreated: {outdir_curnpl}")

    elif not os.path.exists(outdir_curnpl):
        os.makedirs(outdir_curnpl)
        print(f"\nCreated folder: {outdir_curnpl}")

    ######################### Generic Set-up part
    # Define the initial point
    x0 = np.hstack([p.sample() for p in pta.params])

    # Initialize the parameter covariance matrix used by SCAM and AM Jump Proposals
    ndim = len(x0)
    cov = np.diag(np.ones(ndim) * 0.01**2)

    # Initialize the PTMCMC objectz
    sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, outDir=outdir_curnpl)

    # Save set-up information
    save_runtime_info(pta, outdir=outdir_curnpl)

    ######################### Additionals Jump Proposals
    # Set Jump Proposals
    jp = JumpProposal(pta, None, empirical_distr=None)

    # always add draw from prior
    sampler.addProposalToCycle(jp.draw_from_prior, 5)

    # Jump proposals from priors of selected params
    sel_sig = { 
        "gw":30
    }
    for sig, val in sel_sig.items():
        if any([sig in p for p in pta.param_names]):
            sampler.addProposalToCycle(jp.draw_from_par_prior(sig), val)

    ######################### Sample !
    print(f"\nPTMCMC run started with {nsamples} iterations.")
    sampler.sample(x0, int(nsamples), SCAMweight=30, AMweight=15, DEweight=50)

Let it run for a bit...

4.2 Parameter estimation on the HD power-law PSD model

You can also change the gwb_model to "hd_pl", and see how adding "HD" correlations makes it way longer !

overwrite_dir = False
gwb_model = "hd_pl"
outdir_hdpl = f"{datadir}/chains/HD_pl/"

nsamples = 1e5

selpsrs, pta = set_pta_enterprise(psrs, gwb_model, verbose=True)

if os.path.exists(outdir_hdpl) and not overwrite_dir:
    print(f"\nOutdir already exists: {outdir_hdpl}")
else:
    if os.path.exists(outdir_hdpl) and overwrite_dir:
        shutil.rmtree(outdir_hdpl)
        os.makedirs(outdir_hdpl)
        print(f"\nCleared and recreated: {outdir_hdpl}")

    elif not os.path.exists(outdir_hdpl):
        os.makedirs(outdir_hdpl)
        print(f"\nCreated folder: {outdir_hdpl}")

    ######################### Generic Set-up part
    # Define the initial point
    x0 = np.hstack([p.sample() for p in pta.params])

    # Initialize the parameter covariance matrix used by SCAM and AM Jump Proposals
    ndim = len(x0)
    cov = np.diag(np.ones(ndim) * 0.01**2)

    # Initialize the PTMCMC objectz
    sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, outDir=outdir_hdpl)

    # Save set-up information
    save_runtime_info(pta, outdir=outdir_hdpl)

    ######################### Additionals Jump Proposals
    # Set Jump Proposals
    jp = JumpProposal(pta, None, empirical_distr=None)

    # always add draw from prior
    sampler.addProposalToCycle(jp.draw_from_prior, 5)

    # Jump proposals from priors of selected params
    sel_sig = { 
        "gw":30
    }
    for sig, val in sel_sig.items():
        if any([sig in p for p in pta.param_names]):
            sampler.addProposalToCycle(jp.draw_from_par_prior(sig), val)

    ######################### Sample !
    print(f"\nPTMCMC run started with {nsamples} iterations.")
    sampler.sample(x0, int(nsamples), SCAMweight=30, AMweight=15, DEweight=50)

Just stop it, it might finish next month...

The key takeaway is that recovering the full Hellings-Downs signal in a Bayesian framework requires dedicated high-performance computing (e.g. GPU-accelerated likelihoods or HPC clusters), not a single laptop core. In practice, PTA collaborations run these analyses over weeks on clusters and now use Discovery on GPUs.

4.3 Parameter estimation on the HD free-spectrum PSD model.

Rather than constraining a single power-law, a free-spectrum model recovers the GWB power independently at each frequency bin. This gives a model-agnostic view of the spectrum.

gwb_model = "curn_fs"
outdir_curnfs = f"{datadir}/chains/CURN_fs/"

selpsrs, pta = set_pta_enterprise(psrs, gwb_model, verbose=True)
run_PTMCMC(outdir_curnfs, pta, nsamples=1e5, SCAMweight=30, AMweight=15, DEweight=50)

4.4 Visualizing the recovered posterior distributions

Once the MCMC chains have converged, we read them back, discard the first 30% as burn-in, and visualize the posterior in two ways:

1. Trace plots: show the chain value of each parameter vs. iteration. A well-mixed chain should look like white noise around a stable mean, with no long-range drift.
2. Corner plot: shows the 1D and 2D marginal posteriors for all parameter pairs. Tight, unimodal contours centered near the injected values confirm a successful recovery.

Let us focus on the CURN power-law model in this part.

Injected values: $\log_{10} A = -15$, $\gamma = 4.33$ (SMBHB power-law). The recovered posterior should be consistent with these.

ch_curnpl, pars_curnpl = read_all_chains(outdir_curnpl, NRuns='auto', burnin=0.3)

chain_to_plot = ch_curnpl ; pars_to_plot = pars_curnpl ; titleplot = "CURN Power-law PSD"

for i, p in enumerate(pars_to_plot):
    plt.figure(figsize=(8,5))
    plt.title(titleplot, fontsize=20)
    plt.plot(chain_to_plot[:,i])
    plt.xlabel("iterations", fontsize=16)
    plt.ylabel(p, fontsize=16)
    plt.show()

corner.corner(chain_to_plot[:,:-4], labels=pars_to_plot, truths=[4.33, -15], color='cornflowerblue', truth_color='k', show_titles=True)
plt.show()

4.5 Visualizing the recovered GWB spectrum

We now overlay three representations of the recovered GWB spectrum on the same plot:

• Injected power-law (black): the ground truth used to simulate the data
• Recovered CURN power-law (coral): 1000 draws from the posterior of $(\log_{10}A, \gamma)$ converted to $\rho(f) = \sqrt{S_h(f)/T}$, the RMS timing residual contribution per frequency bin
• Recovered CURN free-spectrum (blue violins): the posterior of $\rho$ at each frequency bin, each plotted as a violin. These are model-agnostic and should trace the injected power-law.

What to look for: The free-spectrum violins should align with the injected power-law at low frequencies, and become consistent with the white noise level at high frequencies where the GWB signal is sub-dominant. The CURN power-law posterior should encompass the injected values.

outdir_curnpl = f"{datadir}/chains/CURN_pl_precomputed/"
outdir_curnfs = f"{datadir}/chains/CURN_fs_precomputed/"

ch_curnpl, pars_curnpl = read_all_chains(outdir_curnpl, NRuns='auto', burnin=0.3)
ch_curnfs, pars_curnfs = read_all_chains(outdir_curnfs, NRuns='auto', burnin=0.3)

Tspan = np.max([psr.toas.max()-psr.toas.min() for psr in psrs])
freqs = np.linspace(1/Tspan, 30/Tspan, 30)

plt.figure(figsize=(12,7))

#################
### Plot injected powerlaw
rho_rn_inj = compute_rho(log10_A=np.log10(1e-15), gamma=4.33, f=freqs, T=Tspan)
plt.plot(freqs, rho_rn_inj, color='k', label="Injected GWB", lw=1)

#################
### Get WN level, as S(f)=2 * median(sigma)^2 * median(cadence) (single-pulsar) and S(f)=2 * median(sigma)^2 * median(cadence) / Npsr

Npsrs = len(psrs)
all_sigmas = np.array([np.median(psr.toaerrs) for psr in psrs])
all_cadences = np.array([np.median(psr.toas[1::2] - psr.toas[::2]) for psr in psrs])
sigma_eff = np.median(all_sigmas)
cadence_eff = np.median(all_cadences)
# Single-pulsar WN floor
rho_wn_single = np.sqrt(2 * sigma_eff**2 * cadence_eff / Tspan)
plt.plot(freqs, np.repeat(rho_wn_single, len(freqs)), '--', color='r', alpha=1, label=f"Single-pulsar WN level")
# CURN: averages over Npsrs auto-correlations
rho_wn_curn = rho_wn_single / np.sqrt(Npsrs)
plt.plot(freqs, np.repeat(rho_wn_curn, len(freqs)), ':', color='k', alpha=1, label=f"CURN WN level ({Npsrs} pulsars)")

#################
### Plot powerlaw

# Define the indices to use from the chain
N = 1000
idxs = np.random.choice(range(len(ch_curnpl)), size=N, replace=False)
As = ch_curnpl[:, list(pars_curnpl).index('gw_curn_pl_log10_A')]
Gs = ch_curnpl[:, list(pars_curnpl).index('gw_curn_pl_gamma')]

for i in range(N):
    rho_rn = compute_rho(log10_A=As[i], gamma=Gs[i], f=freqs, T=Tspan)
    if i==0:
        plt.plot(freqs, rho_rn, color='coral', alpha=.02, label="recovered CURN Power-law", zorder=1)
    else:
        plt.plot(freqs, rho_rn, color='coral', alpha=.02, zorder=1)

#################
### Plot free-spectrum

rhos = ch_curnfs[:,:-4]
parts = plt.violinplot(10**rhos, positions=freqs, widths=freqs*0.07, showextrema=False)
for pc in parts['bodies']:
    pc.set_facecolor('cornflowerblue')
    # pc.set_edgecolor('black')
    pc.set_alpha(0.6)
    pc.set_zorder(0)
    
handles, labels = plt.gca().get_legend_handles_labels()
handles.append(Patch(facecolor='cornflowerblue', alpha=0.6, label="Recovered CURN Free spectrum"))
labels.append("Recovered CURN Free spectrum")

#################
### Others
plt.legend(handles=handles, labels=labels, fontsize=14)
plt.axvline(1/86400/365.25, c='k', ls='--', lw=1, alpha=.2)
plt.xscale('log')
plt.yscale('log')
plt.xlabel("Frequency [Hz]", fontsize=20)
plt.ylabel(r"$\rho$ [s]", fontsize=20)
plt.tick_params(labelsize=18)
plt.text(1/86400/365.25 - 3.9e-9, 1.7e-7, "1 $\mathrm{yr}^{-1}$", rotation=90, fontsize=20)
plt.grid(which='both', alpha=.2)
plt.show()

Some selected references related to Section 4

• The Nanohertz Gravitational Wave Astronomer (Taylor, 2021), very useful book about PTA... again related to this part !!
• The International Pulsar Timing Array second data release: Search for an isotropic gravitational wave background (Antoniadis et al. 2022), GWB search with the IPTA DR2
• Comparing Recent Pulsar Timing Array Results on the Nanohertz Stochastic Gravitational-wave Background (Agazie et al. 2025), Comparison of results from multiple PTAs that were published in 2023
• The second data release from the European Pulsar Timing Array — III. Search for gravitational wave signals (Antoniadis et al., EPTA (2023), GWB search with the EPTA DR2 and InPTA
• The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background (Agazie et al., 2023), GWB search with the NANOGrav 15yr
• Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array (Reardon et al., 2023), GWB search with PPTA DR3
• Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array (Xu et al.,, 2023), GWB search with the CPTA DR1
• The MeerKAT Pulsar Timing Array: the first search for gravitational waves with the MeerKAT radio telescope (Miles et al., 2025), GW search with MeerKAT
• PTMCMCSampler: Parallel tempering MCMC sampler (Ellis & van Haasteren, 2017), reference for PTMCMCSampler

5

5. Optimal Statistics (OS)

</h1>

Useful functions

# This function computes the weighted average and uncertainty
# for a set of values (rho) with associated uncertainties (sig).
# It returns the weighted average and the corresponding standard error.
def weightedavg(rho, sig):
    weights, avg = 0., 0.
    for r,s in zip(rho,sig):
        weights += 1./(s*s)
        avg += r/(s*s)
        
    return avg/weights, np.sqrt(1./weights)


# This function bins the cross-correlation data (rho) and associated uncertainties (sig)
# according to angular separation values (xi), using bins defined by zeta.
# Each bin spans 10 degrees in angle. It returns the binned weighted average
# cross-correlation and corresponding uncertainties.
def bin_crosscorr_npairs_tmp(zeta, xi, rho, sig):
    """
    Bin cross-correlations into equal-number of pulsar pairs
    """
    
    rho_avg, sig_avg = np.zeros(len(zeta)), np.zeros(len(zeta))
    
    for i,z in enumerate(zeta[:-1]):
        myrhos, mysigs = [], []
        for x,r,s in zip(xi,rho,sig):
            if x >= z and x < (z+10.):
                myrhos.append(r)
                mysigs.append(s)
        rho_avg[i], sig_avg[i] = weightedavg(myrhos, mysigs)
        
    return rho_avg, sig_avg

def bin_crosscorr_npairs(xi_rad, rho, sig, NPairPerBin=10):
    """
    Bin cross-correlations into equal number of pulsar pairs
    """
    # sort
    idx       = np.argsort(xi)
    xi_sorted   = xi_rad[idx]
    rho_sorted  = rho[idx]
    sig_sorted  = sig[idx]

    xi_mean, xi_err, rho_avg, sig_avg = [], [], [], []
    i = 0
    while i < len(xi_sorted):
        xi_mean.append(np.mean(xi_sorted[i:i+NPairPerBin]))
        xi_err.append(np.std(xi_sorted[i:i+NPairPerBin]))
        r, s = weightedavg(rho_sorted[i:i+NPairPerBin], sig_sorted[i:i+NPairPerBin])
        rho_avg.append(r)
        sig_avg.append(s)
        i += NPairPerBin
    
    return (np.array(xi_mean), np.array(xi_err), np.array(rho_avg), np.array(sig_avg))

def bin_crosscorr_binwidth(xi_rad, rho, sig, bin_size_deg=10.0):
    """
    Bin cross-correlations into equal width angular bins.
    
    Parameters
    ----------
    xi_rad      : array, angular separations in radians
    rho         : array, cross-correlations (normalized)
    sig         : array, uncertainties
    bin_size_deg: float, bin width in degrees
    
    Returns
    -------
    xi_mean, xi_err, rho_avg, sig_avg : arrays of shape (nbins,)
    """
    # sort
    idx       = np.argsort(xi)
    xi_sorted   = xi_rad[idx]
    rho_sorted  = rho[idx]
    sig_sorted  = sig[idx]
    
    bin_size_rad = bin_size_deg * np.pi / 180
    
    xi_deg = xi_sorted

    bin_edges  = np.arange(0, np.pi + bin_size_rad, bin_size_rad)
    bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])

    xi_mean, xi_err, rho_avg, sig_avg = [], [], [], []
    for lo, hi in zip(bin_edges[:-1], bin_edges[1:]):
        mask = (xi_deg >= lo) & (xi_deg < hi)

        if mask.sum() == 0:
            continue  # skip empty bins

        xi_mean.append(np.mean(xi_deg[mask]))
        xi_err.append(np.std(xi_deg[mask]))
        r, s = weightedavg(rho_sorted[mask], sig_sorted[mask])
        rho_avg.append(r)
        sig_avg.append(s)

    return (np.array(xi_mean), np.array(xi_err), np.array(rho_avg), np.array(sig_avg))

def get_histo_equal_Npairs(xi_mean, NPairPerBin):
    idx       = np.argsort(xi_mean)
    xi_sorted   = xi_mean[idx]
    bin_counts  = [NPairPerBin] * len(xi_sorted)
    
    # Estimate bin edges as midpoints between consecutive xi_mean values
    bin_edges        = np.zeros(len(xi_sorted) + 1)
    bin_edges[1:-1] = (xi_sorted[:-1] + xi_sorted[1:]) / 2   # midpoints between consecutive means
    bin_edges[0]     = xi_sorted[0] - 0.5 * (xi_sorted[1]  - xi_sorted[0])   # extrapolate left
    bin_edges[-1]    = xi_sorted[-1] + 0.5 * (xi_sorted[-1] - xi_sorted[-2])  # extrapolate right
    bin_widths       = np.diff(bin_edges)
    bin_centers      = bin_edges[:-1] + bin_widths / 2  # centers from edges, not xi_mean    
    
    return bin_centers, bin_widths, bin_counts

def get_histo_equal_BinWidth(xi_rad, bin_size_deg):
    xi_deg = xi_rad * 180 / np.pi
    
    bin_edges   = np.arange(0, 180 + bin_size_deg, bin_size_deg)  # edges from 0 to 180
    bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])          # centers at 5, 15, 25, ...
    bin_counts, _ = np.histogram(xi_deg, bins=bin_edges)           # use actual edges
    
    return bin_centers * np.pi / 180, bin_counts

def GPfit(x, y, yerr, x_output):
    kernel = RBF(length_scale=30) + WhiteKernel()
    gp     = GaussianProcessRegressor(kernel=kernel, alpha=yerr**2)
    gp.fit(x.reshape(-1,1), y)

    y_fit, y_std_gp = gp.predict(x_smooth.reshape(-1, 1), return_std=True)

    # Interpolate measurement noise onto smooth grid
    yerr_smooth = np.interp(x_smooth, x, yerr)

    # Total uncertainty = GP epistemic + measurement noise in quadrature
    y_std = np.sqrt(y_std_gp**2 + yerr_smooth**2)

    return y_fit, y_std

def build_Q_matrix_defiant(OS_obj, params):
    """
    Build the normalized Q matrix such that OS = d^T Q d for d ~ N(0,I).
    Follows the Hazboun et al. construction (QtildeIJ block matrix).
    """
    npsr = OS_obj.npsr
    nfreq = OS_obj.nfreq
    ngw = 2 * nfreq

    # GW spectrum
    phi = OS_obj._get_phi(params)   # shape (ngw,)
    sPhi = np.sqrt(phi)

    # ORF values for HD
    orf_full = OS_obj.orf_matrix[0]  # shape (npsr, npsr)
    pairs = [(i, j) for i in range(npsr) for j in range(i+1, npsr)]
    orfs = np.array([orf_full[i, j] for (i, j) in pairs])

    # noise-marginalized GW overlap matrices
    # FCF[i] = F^T C_i^{-1} F (marginalized over timing model)
    FCFs = []
    for i in range(npsr):
        FNF    = OS_obj._get_FNF(i, params)   # (ngw, ngw)
        FNT    = OS_obj._get_FNT(i, params)   # (ngw, ntim)
        TNT    = OS_obj._get_TNT(i, params)   # (ntim, ntim)
        phiinv = OS_obj._get_phiinv(i, params) # (ntim, ntim)
        inner  = phiinv + TNT
        FCF    = FNF - FNT @ np.linalg.solve(inner, FNT.T)
        FCFs.append(FCF)

    # P_i = sPhi @ FCF[i] @ sPhi  (= D_i in Hazboun notation)
    Ps = [sPhi[:, None] * FCF * sPhi[None, :] for FCF in FCFs]

    # Cholesky factors: P_i = L_i L_i^T
    Ls = []
    for P in Ps:
        P = 0.5 * (P + P.T)
        eps = 1e-10 * np.trace(P) / P.shape[0]
        L = np.linalg.cholesky(P + eps * np.eye(ngw))
        Ls.append(L)

    # normalization: bottom = sum_{i<j} tr(Pinv_i S_ij Pinv_j S_ij^T)
    # in our notation: bottom = sum_{i<j} orf^2 * tr(P_i P_j)
    # = sum_{i<j} orf^2 * tr(L_i^T L_j L_j^T L_i)
    bottom = 0.0
    for w, (i, j) in zip(orfs, pairs):
        LiTLj = Ls[i].T @ Ls[j]          # (ngw, ngw)
        bottom += w**2 * np.trace(LiTLj.T @ LiTLj)
    
#     norm = 1.0 / np.sqrt(bottom)
    norm = 1.0 / (2.0 * np.sqrt(bottom))
#     print(f'norm: {norm:.6f}, bottom: {bottom:.6f}')

    # build block Q matrix
    cnt = npsr * ngw
    inds = [slice(i * ngw, (i + 1) * ngw) for i in range(npsr)]
    Q = np.zeros((cnt, cnt))

    for w, (i, j) in zip(orfs, pairs):
        # Q_ij = norm * orf_ij * L_i^T L_j  (off-diagonal block)
        Bij = norm * w * (Ls[i].T @ Ls[j])
        Q[inds[i], inds[j]] += Bij
        Q[inds[j], inds[i]] += Bij.T

    # sanity checks
#     print(f'Q symmetric: {np.max(np.abs(Q - Q.T)):.2e}')
#     print(f'tr(Q): {np.trace(Q):.6f}             # should be 0')
#     print(f'sqrt(2*tr(Q^2)): {np.sqrt(2*np.trace(Q@Q)):.6f}  # should be 1')

    return Q

def gx2_survival_hybrid(snr_values, eigs, crossover_snr=0.0):
    """
    Use saddlepoint (Lugannani-Rice) for all SNR values.
    Imhof is unreliable with the current eigenvalue set.
    The saddlepoint is inaccurate near SNR~0 but exact in the tail,
    which is the only region that matters for p-value reporting.
    """
    from scipy import optimize
    
    eigs = np.array(eigs)
    t_max = 0.999 / (2 * np.max(eigs))
    t_min = 0.999 / (2 * np.min(eigs)) if np.min(eigs) < 0 else -10.0

    def K(t):
        return -0.5 * np.sum(np.log(1 - 2 * eigs * t))
    def K1(t):
        return np.sum(eigs / (1 - 2 * eigs * t))
    def K2(t):
        return np.sum(2 * eigs**2 / (1 - 2 * eigs * t)**2)

    results = []
    methods = []

    for s in snr_values:
        try:
            t_hat = optimize.brentq(lambda t: K1(t) - s, t_min, t_max,
                                    xtol=1e-12, maxiter=1000)
            w = np.sign(t_hat) * np.sqrt(2 * (t_hat * s - K(t_hat)))
            u_sp = t_hat * np.sqrt(K2(t_hat))
            if np.abs(w) > 1e-6:
                p = float(stats.norm.sf(w) + stats.norm.pdf(w) * (1/w - 1/u_sp))
                p = np.clip(p, 0, 1)
            else:
                p = 0.5
        except Exception as e:
            p = np.nan
        results.append(p)
        methods.append('saddlepoint')

    return np.array(results), methods

inj_params = {
    f'gw_curn_pl_gamma':4.33,
    f'gw_curn_pl_log10_A': -15
}

5.1 Motivation and description

The Bayesian approach is powerful but computationally expensive, especially when spatial correlations (HD ORF) are included. The Optimal Statistic (OS) is a frequentist cross-correlation estimator that:

• Is fast to compute (no MCMC required for the HD step itself)
• Provides an amplitude estimate $\hat{A}^2$ and its uncertainty $\sigma(\hat{A}^2)$
• Enables a frequentist significance estimate via the signal-to-noise ratio $\hat{\rho} = \hat{A}^2 / \sigma(\hat{A}^2)$
• Allows per-pair cross-correlation recovery, directly tracing the Hellings-Downs curve

The OS estimator for the GWB amplitude squared $\hat{A}^2$ is:

$$\hat{A}^2 = \frac{\displaystyle\sum_{a<b} \delta t_a^T C_a^{-1} \tilde{S}_{ab} C_b^{-1} \delta t_b}{\displaystyle\sum_{a<b} \mathrm{tr}\!\left(C_a^{-1} \tilde{S}_{ab} C_b^{-1} \tilde{S}_{ba}\right)}$$

where:

• $\delta t_a$ are the timing residuals of pulsar $a$
• $C_a$ is the noise covariance matrix of pulsar $a$ (including white noise and pulsar-intrinsic red noise)
• $\tilde{S}_{ab}(f) = \Gamma_{ab} \cdot S_h(f)$ is the expected cross-power, with $\Gamma_{ab}$ the HD coefficient for pair $(a,b)$

In the noise-dominated (weak signal) regime, $\langle \hat{A}^2 \rangle = 0$, while the variance of the estimator is:

$$\sigma(\hat{A}^2) = \left[\sum_{a<b} \mathrm{tr}\!\left(C_a^{-1} \tilde{S}_{ab} C_b^{-1} \tilde{S}_{ba}\right)\right]^{-1/2}$$

Key subtlety: the noise marginalization: The single-parameter OS above uses fixed noise parameters. In practice, the noise covariance $C_a$ depends on estimated white-noise and red-noise parameters, which themselves carry uncertainty. The Noise-Marginalized OS (NMOS) marginalizes over this uncertainty by drawing noise parameters from the posterior of a noise-only MCMC chain and averaging the OS over those draws, this is what compute_OS(N=100, ...) does below.

5.2 Amplitude recovery and signal-to-noise ratio

We now build the OS object using defiant.OptimalStatistic, feeding it the CURN power-law MCMC chain so that the noise marginalization can draw from the posterior.

The cells below compute the Noise-Marginalized OS (NMOS) with $N=100$ noise draws. Each draw gives a different $\hat{A}^2$ estimate; the resulting distribution captures both the statistical uncertainty of the OS and the propagated noise-parameter uncertainty.

Output of compute_OS(N=100):

• A2: array of $\hat{A}^2$ values, one per noise draw
• A2s: corresponding $\sigma(\hat{A}^2)$ values
• idx: indices of the noise draws used

# Since we created our PTA object with 'gw' as the name, make sure to set that!
gwb_model = "curn_pl"
outdir_curnpl = f"{datadir}/chains/CURN_pl_precomputed/"

selpsrs, pta = set_pta_enterprise(psrs, gwb_model="curn_pl", verbose=True)
# run_PTMCMC(outdir_curnpl, pta, nsamples=1e5, SCAMweight=30, AMweight=15, DEweight=50, overwrite_dir=True, NRuns=10)

ch, pars = read_all_chains(outdir_curnpl, NRuns='auto')
os_obj = OptimalStatistic(psrs, pta=pta, chain=ch, param_names=list(pars), gwb_name='gw_curn_pl', orfs=['hd'])

# os_obj.set_orf(['hd'])

# When N>1, params is ignored. The default value of params is None, so we are good!
output_NMOS = os_obj.compute_OS(N=100, return_pair_vals=False)
A2, A2s, idx = (output_NMOS[k] for k in ['A2','A2s','idx'])

# A function to implement uncertainty sampling to account for underlying uncertainty in the optimal statistic A^2
full_A2 = utils.uncertainty_sample(A2,A2s,pfos=False,mcos=False)

plt.figure(figsize=(8,5))

plt.hist(A2,bins='auto',histtype='step', density=True, label='A^2 distribution')
plt.hist(full_A2,bins='auto',histtype='step', density=True, label='Full A^2 distribution')

plt.axvline(10**(2*(-15)),linestyle='dashed',color='k',label='Injected')

plt.title('Noise Marginalized Optimal Statistics (NMOS)', fontsize=18)
plt.xlabel('$A^2$', fontsize=16)
plt.ylabel('$p(A^2)$', fontsize=16)
plt.grid()
plt.legend(fontsize=14)
plt.show()

5.3 Per-pair cross-correlations

Rather than a single amplitude, we can compute the OS per angular separation bin, recovering the shape of the inter-pulsar correlation as a function of $\zeta$.

For each angular bin $[\zeta_i, \zeta_{i+1}]$, we define a binned ORF $\hat{\Gamma}_i(\zeta_{ab})$ and compute the amplitude in that bin:

$$\hat{A}^2_i = \frac{\sum_{a<b,\, \zeta_{ab}\in\text{bin }i} \delta t_a^T C_a^{-1} \tilde{S}_{ab} C_b^{-1} \delta t_b} {\sum_{a<b,\, \zeta_{ab}\in\text{bin }i} \mathrm{tr}\!\left(C_a^{-1} \tilde{S}_{ab} C_b^{-1} \tilde{S}_{ba}\right)}$$

Plotting $\hat{A}^2_i$ vs. $\zeta_i$ and comparing to the theoretical HD curve $\Gamma(\zeta)$ is the most direct visual evidence for the spatial correlations expected from a GWB.

Let us start with the enterprise_extensions' method

## Here we compute the Optimal Statistic and the per-pair cross-correlations
OS_ee = ostat.OptimalStatistic(psrs, pta=pta, orf='hd')

# xi: angular separation [rad] for each pulsar pair
# rho: correlation coefficient for each pulsar pair
# sig: 1-sigma uncertainty on correlation coefficient for each pulsar pair.
# OS: Optimal statistic value (units of A_gw^2)
# OS_sig: 1-sigma uncertainty on OS
xi, rho, sig, OS, OS_sig = OS_ee.compute_os(params=inj_params)

# Normalize by OS to get correlation in [-1, 1]
rho_norm = rho / OS
sig_norm  = sig / OS

# adjust as wanted

plotype = 1 # 1 = Equal Npairs, 2 = Equal Bin Width

NPairPerBin = 200 # adjust as wanted - Used if plotype==1
bin_widths = 5 # degree - adjust as wanted - Used if plotype==2

##############

if plotype == 1:
    xi_mean, xi_err, rho_avg, sig_avg = bin_crosscorr_npairs(xi, rho_norm, sig_norm, NPairPerBin=NPairPerBin)
    bin_centers, bin_widths, bin_counts = get_histo_equal_Npairs(xi_mean, NPairPerBin)
elif plotype == 2:
    xi_mean, xi_err, rho_avg, sig_avg = bin_crosscorr_binwidth(xi, rho_norm, sig_norm, bin_size_deg=bin_widths)
    bin_centers, bin_counts = get_histo_equal_BinWidth(xi, bin_widths)
    bin_widths *= np.pi / 180

fig, ax1 = plt.subplots(figsize=(10, 6))

if plotype==1:
    add_title = f"Equal N pairs per bin - {NPairPerBin} pairs per bin"
if plotype==2:
    add_title = f"Equal bin widths - {bin_widths * 180/np.pi:.1f} degrees per bin"
    
fig.suptitle(f"OS correlation - {add_title}", fontsize=18)

# ---- Left y-axis: correlations ----
##############################
ax1.errorbar(xi_mean*180/np.pi, rho_avg, xerr=xi_err*180/np.pi/2, yerr=sig_avg,
             ls='', color='cornflowerblue', fmt='o',
             capsize=4, elinewidth=1.2, zorder=2, label='OS data')

##############################
x_smooth = np.linspace(0.01, 180, 500)
y_fit, y_std = GPfit(x=xi_mean * 180/np.pi, y=rho_avg, yerr=sig_avg, x_output=x_smooth)
ax1.plot(x_smooth, y_fit, color='green', lw=3, label='GP fit', zorder=4, alpha=.5)
ax1.fill_between(x_smooth,
                 y_fit - y_std,
                 y_fit + y_std,
                 color='green', alpha=0.1, zorder=0, label='GP uncertainty')

##############################
zeta = np.linspace(0.01, 180, 100)
HD = get_HD_curve(zeta + 1)
ax1.plot(zeta, HD, ls='--', label='Hellings-Downs', color='gray', lw=3)


##############################
ax1.axhline(0, c='k', zorder=0, alpha=.6)

##############################
ax1.set_xlabel('Angle of separation (degrees)', fontsize=16)
ax1.set_ylabel('Correlation', fontsize=16)
ax1.grid(alpha=.2)


# ---- Right y-axis: number of pairs per bin ----
###################################
ax2 = ax1.twinx()
ax2.bar(bin_centers * 180/np.pi, 
        bin_counts, 
        edgecolor='black', 
        width=bin_widths * 180/np.pi, 
        alpha=0.1, color='gray', zorder=1, label='Number of pairs per bin')

ax2.set_ylabel('Number of pairs', fontsize=16)

if plotype==1:
    ax2.set_ylim(0, max(bin_counts) * 5)  # push histogram to bottom so it doesn't crowd correlations
elif plotype==2:
    ax2.set_ylim(0, max(bin_counts) * 1.5)  # push histogram to bottom so it doesn't crowd correlations
    
# ---- Legends ----
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, fontsize=13)

plt.tight_layout()

plt.show()

Let us now plot it with Defiant tools

ch, pars = read_all_chains(outdir_curnpl, NRuns='auto')
os_obj = OptimalStatistic(psrs, pta=pta, chain=ch, param_names=list(pars), gwb_name='gw_curn_pl', orfs=['hd'])

Nbins = 15

# Set our ORF to HD
os_obj.set_orf(['hd'])
hd_orf = orf_functions.get_orf_function('hd') 
xi_range = np.linspace(0,np.pi,1000)[1:] # Define our range of pulsar separations
hd_mod = hd_orf(xi_range)

# If params=None, then DEFIANT will use maximum likelihood values from OS_obj.lfcore
# xi,rho,sig,C,A2,A2s,idx = OS_obj.compute_OS(params=None)
output_OS = os_obj.compute_OS(params=None)
A2_OS = output_OS['A2']
A2s_OS = output_OS['A2s']
idx_OS = output_OS['idx']
xi_OS = output_OS['xi']
rho_OS = output_OS['rho']
sig_OS = output_OS['sig']
C_OS = output_OS['C']

# Plot the binned pair correlation plot. 
f, ax = defplot.create_correlation_plot(xi_OS, rho_OS, sig_OS, C_OS, A2_OS, A2s_OS, bins=Nbins, figsize=(10,6))
ax.plot(xi_range,10**(2*inj_params[f'gw_curn_pl_log10_A'])*hd_mod,'--k',label='HD prediction')

ax.xaxis.label.set_size(18)
ax.yaxis.label.set_size(18)
plt.title('OS correlation vs. HD curve', fontsize=20)

plt.legend()
# plt.savefig("os.png")
# plt.close()
plt.show()

5.5 The Pair-Covariance OS (PCOS)

The standard OS treats pulsar pairs as statistically independent, which is an approximation: in reality, pairs sharing a pulsar are correlated through the common noise of that pulsar. The Pair-Covariance OS (PC+OS) accounts for these inter-pair correlations by computing the full covariance matrix $\mathbf{C}$ of the pair cross-correlations.

This improves the amplitude estimate when the GWB is strong relative to the noise, but at the cost of computing an $N_{\rm pairs} \times N_{\rm pairs}$ covariance matrix.

⏳ This cell takes ~15 minutes. It is included for completeness, feel free to run it after the tutorial. The standard OS and NMOS results are sufficient for the main analysis.

"""
PCOS - Takes ~15 min !
"""
print("DOING PC+OS...")

output_PCOS = os_obj.compute_OS(inj_params, pair_covariance=True)

Nbins = 15

A2 = output_PCOS['A2']
A2s = output_PCOS['A2s']
idx = output_PCOS['idx']
xi = output_PCOS['xi']
rho = output_PCOS['rho']
sig = output_PCOS['sig']
C = output_PCOS['C']

# Plot the binned pair correlation plot. 
f, ax = defplot.create_correlation_plot(xi_OS, rho_OS, sig_OS, C_OS, A2_OS, A2s_OS, bins=Nbins, figsize=(10,6))
ax.plot(xi_range,10**(2*inj_params[f'gw_curn_pl_log10_A'])*hd_mod,'--k',label='HD prediction')

ax.xaxis.label.set_size(18)
ax.yaxis.label.set_size(18)
plt.title('OS correlation vs. HD curve', fontsize=20)

plt.title('PC+OS')
plt.legend()
plt.show()

5.6 Measuring the significance

We can also use the OS approach to estimate the significance of the GWB signal. Let us start with the signal-to-noise ratio.

5.6.1 The signal-to-noise ratio (S/N)

The OS signal-to-noise ratio is simply: $$\hat{\rho} = \frac{\hat{A}^2}{\sigma(\hat{A}^2)}$$ A value $\hat{\rho} \gg 1$ indicates that the measured cross-correlations are inconsistent with noise alone.

Here we use the NMOS output to obtain a distribution of $\hat{\rho}$ values, one per noise draw, whose median summarizes our detection significance.

A2 = output_NMOS['A2']
A2s = output_NMOS['A2s']

# S/N per NM iteration, no uncertainty sampling needed
snr_nm = A2 / A2s
snr_nm_mean = np.median(snr_nm)

plt.figure(figsize=(8, 5))
plt.hist(snr_nm, bins='auto', histtype='step', density=True)
plt.axvline(snr_nm_mean, ls='--', color='k', label=f'Median S/N = {np.median(snr_nm):.2f}')
plt.xlabel(r'S/N $= A^2 / \sigma_{A^2}$', fontsize=16)
plt.ylabel(r'$p(\mathrm{S/N})$', fontsize=16)
plt.title('Noise Marginalized Optimal Statistics (NMOS)', fontsize=18)
plt.grid()
plt.legend(fontsize=14)
plt.show()

5.6.2 The p-value

From a S/N, we can compute a p-value, the probability of obtaining a S/N at least as large as observed, assuming the null hypothesis (no GWB) is true:

$$p = P(\hat{\rho} \geq \hat{\rho}_{\rm obs} \mid H_0)$$

The smaller the p-value, the higher the significance. The tricky part for PTAs is that we need to build a null hypothesis distribution empirically, since the GWB signal is always present in the data (i.e., it is present in the timing residuals of all pulsars). Two standard approaches are used to construct this null distribution:

Phase shifts

The timing residuals of each pulsar are phase-shifted in the Fourier domain by a random offset $\phi_a \in [0, 2\pi)$, independently for each pulsar. This destroys the inter-pulsar correlations (and thus the HD signature) while preserving the spectral properties of the noise and any common red process. The OS is then recomputed on each shifted dataset, building a distribution of $\hat{\rho}$ under $H_0$.

Sky scrambles

Instead of shifting phases, the sky positions of pulsars are randomly reassigned (or pulsars are swapped between each other), which changes the ORF values $\Gamma_{ab}$ to effectively random values inconsistent with HD. This destroys the spatial correlation pattern while keeping the noise properties intact. The OS is recomputed for each scrambled configuration, again building the null $\hat{\rho}$ distribution.

Both methods can be combined simultaneously, applying phase shifts and sky scrambles together, to more efficiently decorrelate the data from any HD-like pattern.

The p-value is then estimated as: $$p = \frac{\#\{\hat{\rho}_{\rm null} \geq \hat{\rho}_{\rm obs}\}}{N_{\rm null}}$$ where $\hat{\rho}_{\rm null}$ denotes S/N values from either phase-shifted, sky-scrambled (or both applied) datasets.

### Building the null distribution with "phase shifts" To convert the measured S/N into a **p-value**, we need to know the distribution of $\hat{\rho}$ under the null hypothesis $H_0$ (no GWB). Since the GWB signal is always present in simulated data, we cannot simply remove it, instead, we **destroy the HD spatial correlations** while preserving all spectral properties. In a **phase shift**, the timing residuals of each pulsar are independently shifted by a random phase $\phi_a \in [0, 2\pi)$ in the Fourier domain. This scrambles the relative phase between pulsars and erases any coherent inter-pulsar correlation, while keeping the single-pulsar noise properties unchanged.

Note: we use fixed (injected) noise parameters here rather than marginalizing. Noise marginalization could be included by sampling params from the MCMC chain for each shift, but this increases the computational cost significantly.

n_shifts = 500

# p_Phase (float): The p-value of the OS
# snr_Phase (float): The measured SNR of the OS
# n_dist_Phase (np.ndarray): The null distribution of the SNR
p_Phase, snr_Phase, n_dist_Phase = phase_shift_OS(os_obj, params=inj_params, n_shifts=n_shifts)
print('Measured P-value:',p_Phase,'Upper limit: P<=',1/(len(n_dist_Phase)+2))

plt.title('OS [HD] S/N and null distribution: phase shift')

g = np.random.randn(int(1e5))
plt.hist(g, bins=100, histtype='step', alpha=0.7, density =True, label="$\mathcal{N}(0,1)$")

plt.hist(n_dist_Phase,bins='auto',density=True,histtype='step',label='Phase-shifted Null');
plt.axvline(snr_Phase, color='r', linestyle='--',label='Measured SNR')
plt.xlabel('SNR')
plt.legend()
plt.show()

### Alternative null distribution: The sky scrambles An alternative to phase shifts is the **sky scramble**: instead of shifting phases, the sky positions of pulsars are randomly reassigned. This changes the ORF values $\Gamma_{ab}$ to effectively random values inconsistent with HD, while all timing properties remain intact. Sky scrambles are particularly powerful because they test specifically the **angular pattern** of the correlations, any signal whose spatial signature is not HD-like should not survive scrambling.

Phase shifts vs. sky scrambles: Both methods are valid and complementary. Phase shifts operate in the time/frequency domain; sky scrambles operate in the angular domain. Using both (or combining them) gives a more robust null distribution. The 2023 PTA papers used both approaches to cross-check significance estimates. These methods were highly discussed since 2023, but we won't go into the details in this tutorial.

n_scrambles = 500

# p_Sky (float): The p-value of the OS
# snr_Sky (float): The measured SNR of the OS
# n_dist_Sky (np.ndarray): The null distribution of the SNR
p_Sky, snr_Sky, n_dist_Sky = sky_scramble_OS(os_obj, inj_params, n_scrambles=n_scrambles, swap_pos=False)
print('Measured P-value:',p_Sky,'Upper limit: P<=',1/(len(n_dist_Sky)+2))

plt.title('OS (ORF: HD) S/N and null distribution: sky scramble')

g = np.random.randn(int(1e5))
plt.hist(g, bins=100, histtype='step', alpha=0.7, density =True, label="$\mathcal{N}(0,1)$")

plt.hist(n_dist_Sky,bins='auto', density=True, histtype='step',label='Sky-scrambled Null');
plt.axvline(snr_Sky, color='r', linestyle='--',label='Measured SNR')
plt.xlabel('SNR')
plt.legend()
plt.show()

5.6.3 The equivalent Gaussian significance

A common way to communicate a p-value is to convert it to an equivalent number of Gaussian standard deviations, the so-called sigma ($\sigma$) statistics. This is the value $\mathcal{N}_\sigma$ such that a one-sided Gaussian tail gives the same probability as the measured p-value. The higher the number of sigma, the more incompatible the data are with the null hypothesis.

A key subtlety from PTAs is that the null distribution of $\hat{\rho}$ is not Gaussian. Previous discussions of the OS incorrectly assumed that the analytic null distribution of $\hat{\rho}$ is well-approximated by a zero-mean unit-variance Gaussian. In reality, the null distribution has tails that differ significantly from a Gaussian, but which follows a generalized chi-squared (GX2) distribution, i.e. a linear combination of chi-squared distributions. Assuming Gaussianity therefore gives the wrong p-value, and hence a misleading sigma. A correct assessment of the statistical significance requires fitting the null distribution with a GX2 model and deriving the p-value from this fit. The sigma value is then still reported as the Gaussian equivalent of that p-value, a convenient way to communicate significance in familiar units, even though the underlying distribution is not Gaussian.

The standard thresholds in GW astronomy are:

• $\mathcal{N}_\sigma \gtrsim 3$: evidence for a signal ($p \lesssim 1.3 \times 10^{-3}$)
• $\mathcal{N}_\sigma \gtrsim 5$: detection ($p \lesssim 2.9 \times 10^{-7}$)

The 2023 PTA papers (CPTA DR1, EPTA DR2, NANOGrav 15yr, PPTA DR3) reported HD correlations at the $2$–$4\sigma$ level depending on the dataset and method, motivating the use of "compelling evidence for" rather than "detection of" a GWB.

### The GX2 null distribution We now run both phase shifts and sky scrambles with $N=500$ realizations each, and compare them to the analytic **generalized chi-squared (GX2)** null distribution. The GX2 distribution arises because the OS statistic is a quadratic form in Gaussian random variables: its null distribution is a weighted sum of independent $\chi^2_1$ variables, whose weights are the eigenvalues of the OS normalization matrix $Q$. This is computed analytically below using the **saddlepoint approximation** (Lugannani-Rice method), which is accurate in the tails where p-values are small.

Why the null distribution matters for sigma: Converting a p-value to "$n\sigma$" implicitly assumes a Gaussian null distribution. But for PTAs, the OS null distribution has non-Gaussian tails, assuming Gaussianity gives the wrong sigma. The correct procedure is to compute the p-value from the actual (GX2 or empirical) null distribution, then convert that p-value to a Gaussian-equivalent sigma $\mathcal{N}_\sigma = \Phi^{-1}(1-p)$.

# Compute with new Q matrix construction for the GX2. We won't go into the details here, we just place it here as it's long to compute
Q = build_Q_matrix_defiant(os_obj, inj_params)
eigvals_raw = np.linalg.eigvalsh(Q)

Let us start with phase shifts

n_shifts = 500
p_Phase, snr_Phase, n_dist_Phase = phase_shift_OS(os_obj, params=inj_params, n_shifts=n_shifts)

print('Measured P-value:', p_Phase,'Upper limit: P<=',1/(len(n_dist_Phase)+2))

# --- p-values ---
n_phase = len(n_dist_Phase)
p_obs_phase, _ = gx2_survival_hybrid([snr_Phase], eigvals_raw)
p_obs_phase = p_obs_phase[0]
p_empirical_phase = np.mean(n_dist_Phase >= snr_Phase)

p_obs_phase_str = f"{p_obs_phase:.2e} ({stats.norm.isf(p_obs_phase):.1f}σ)"
p_phase_str     = f"{p_empirical_phase:.2e} ({stats.norm.isf(p_empirical_phase):.1f}σ)" if p_empirical_phase > 0 else f"<{1/(n_phase+2):.4f} (>{stats.norm.isf(1/(n_phase+2)):.1f}σ)"

plot_measurements = False
Nsig_max_plot = 5 # try 5 or 23
GX2_SNRmax_plot = 5 # try 5 or 60
Gaussian_SNRmax_plot = 5 # try 5 or 60

########################
fig, ax = plt.subplots(figsize=(9, 6))

## Phase shift null distribution
n_dist_sorted_phase = np.sort(n_dist_Phase)
survival_phase = 1 - np.arange(1, n_phase + 1) / (n_phase + 1)
ax.step(np.append(n_dist_sorted_phase[0], n_dist_sorted_phase),
        np.append(1.0, survival_phase),
        where='post', label=f'Phase shifts  p={p_phase_str}', color='darkorange')

## GX2
x_gx2 = np.linspace(-2, GX2_SNRmax_plot, 100)
survival_gx2, _ = gx2_survival_hybrid(x_gx2, eigvals_raw)
ax.plot(x_gx2, survival_gx2, color='green', label=f'GX2 (saddlepoint)  p={p_obs_phase_str}')

## Gaussian
x = np.linspace(-2, Gaussian_SNRmax_plot, 1000)
ax.plot(x, stats.norm.sf(x), 'k--', label=r'$\mathcal{N}(0,1)$')

## Observed SNR
if plot_measurements:
    ax.axvline(snr_Phase, color='r', linestyle='--', label=f'S/N Phase = {snr_Phase:.1f}')

## Sigma levels
for ns in range(1, Nsig_max_plot + 1):
    p_level = stats.norm.sf(ns)
    ax.axhline(p_level, color='gray', linestyle=':', lw=1)
    ax.text(ax.get_xticks()[-2]-1, p_level*1.2, f' ${ns}\\sigma$', ha='left', va='bottom', fontsize=16)

ax.set_yscale('log')
ax.set_xlabel('SNR', fontsize=16)
ax.set_ylabel('p-value', fontsize=16)
ax.set_title('OS (ORF: HD) — Phase shift null distribution')
ax.tick_params(axis='both', which='major', labelsize=14)
ax.legend(fontsize=10)
ax.grid(which='both', alpha=0.2)
plt.tight_layout()
plt.show()

Now the sky scrambles

n_scrambles = 500

os_obj.set_orf(['hd'])
p_Sky, snr_Sky, n_dist_Sky = sky_scramble_OS(os_obj,inj_params, n_scrambles=n_scrambles, swap_pos=False)

print('Measured P-value:', p_Sky,'Upper limit: P<=',1/(len(n_dist_Sky)+2))

# Calculate p-values
n_sky = len(n_dist_Sky)
p_obs_sky, _ = gx2_survival_hybrid([snr_Sky], eigvals_raw)
p_obs_sky = p_obs_sky[0]
p_empirical_sky = np.mean(n_dist_Sky >= snr_Sky)

p_obs_sky_str = f"{p_obs_sky:.2e} ({stats.norm.isf(p_obs_sky):.1f}σ)"
p_sky_str     = f"={p_empirical_sky:.2e} ({stats.norm.isf(p_empirical_sky):.1f}σ)" if p_empirical_sky > 0 else f"<{1/(n_sky+2):.4f} (>{stats.norm.isf(1/(n_sky+2)):.1f}σ)"

plot_measurements = False
Nsig_max_plot = 5 # try 5 or 23
GX2_SNRmax_plot = 5 # try 5 or 60
Gaussian_SNRmax_plot = 5 # try 5 or 60

########################
fig, ax = plt.subplots(figsize=(9, 6))

## Sky scramble null distribution
n_dist_sorted_sky = np.sort(n_dist_Sky)
survival_sky = 1 - np.arange(1, n_sky + 1) / (n_sky + 1)
ax.step(np.append(n_dist_sorted_sky[0], n_dist_sorted_sky),
        np.append(1.0, survival_sky),
        where='post', label=f'Sky scrambles  p{p_sky_str}', color='steelblue')

## GX2
x_gx2 = np.linspace(-2, GX2_SNRmax_plot, 100)
survival_gx2, _ = gx2_survival_hybrid(x_gx2, eigvals_raw)
ax.plot(x_gx2, survival_gx2, color='green', label=f'GX2 (saddlepoint)  p={p_obs_sky_str}')

## Gaussian
x = np.linspace(-2, Gaussian_SNRmax_plot, 1000)
ax.plot(x, stats.norm.sf(x), 'k--', label=r'$\mathcal{N}(0,1)$')

## Observed SNR
if plot_measurements:
    ax.axvline(snr_Sky, color='darkred', linestyle='--', label=f'S/N Sky = {snr_Sky:.1f}')

## Sigma levels
for ns in range(1, Nsig_max_plot + 1):
    p_level = stats.norm.sf(ns)
    ax.axhline(p_level, color='gray', linestyle=':', lw=1)
    ax.text(ax.get_xticks()[-2]-1, p_level*1.2, f' ${ns}\\sigma$', ha='left', va='bottom', fontsize=16)

ax.set_yscale('log')
ax.set_xlabel('SNR', fontsize=16)
ax.set_ylabel('p-value', fontsize=16)
ax.set_title('OS (ORF: HD) — Sky scramble null distribution')
ax.tick_params(axis='both', which='major', labelsize=14)
ax.legend(fontsize=10)
ax.grid(which='both', alpha=0.2)
plt.tight_layout()
plt.show()

Some selected references related to Section 5

• The Nanohertz Gravitational Wave Astronomer (Taylor, 2021), very useful book about PTA. Interesting, again here.
• Optimal strategies for gravitational-wave stochastic background searches in pulsar timing arrays (Anholm et al., 2022), basis Optimal Statistic
• Noise-marginalized optimal statistic: A robust hybrid frequentist-Bayesian statistic for the stochastic gravitational-wave background in pulsar timing arrays (Vigeland et al., 2018), noise marginalized Optimal Statistics
• Analytic distribution of the optimal cross-correlation statistic for stochastic gravitational-wave-background searches using pulsar timing arrays (Hazboun et al., 2023), reference for the GX2 null distribution.
• Spatial and spectral characterization of the gravitational-wave background with the PTA optimal statistic (Gersbach et al., 2025), Per-frequency Optimal Statistics, relevant publication for defiant, good review on the optimal statistic method

6

6. Bonus: Sensitivity as a Function of Pulsar Number

</h1>

How does the significance grow with the array?

One of the most powerful diagnostics of a PTA detection is to ask: does the significance grow as expected when we add more pulsars?

For an ideal array with $N_{\rm psr}$ pulsars and $N_{\rm pairs} = N_{\rm psr}(N_{\rm psr}-1)/2$ pairs, the expected S/N scales roughly as (Siemens et al. 2013):

$$\hat{\rho} \propto \sqrt{N_{\rm pairs}} \propto N_{\rm psr}$$

Deviations from this scaling can reveal:

• Pulsars with anomalously high noise
• A signal that is not perfectly isotropic
• Systematics in the dataset

Our task

Using the precomputed OS results, let us rank the pulsars by their contribution to the S/N and plot:

• $\hat{\rho}$ vs. $N_{\rm psr}$ (adding pulsars one by one, ranked by sensitivity) and compare with predicted scaling from Siemens et al., 2013

iterperpsr = 50
Npsrlist = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130]

SNRs = []

# Read the chain for CURN PL
ch, pars = read_all_chains(outdir_curnpl, NRuns='auto')

for Npsrs in Npsrlist:
    print(f"Run for {Npsrs} pulsars.")
    
    # The same process is done "iterperpsr" times for each pulsar number
    SNR_per_psr = []
    for i in range(iterperpsr):
        print(f"{(i+1)*100/iterperpsr:.1f}%", end="\r")
        
        # Create pta object for Npsrs randomly chosen pulsars
        selpsrs, pta = set_pta_enterprise(psrs, gwb_model, Npsrs=Npsrs, verbose=False) 
        
        # Compute the Noise-marginalized Optimal Statistics, and save the corresponding S/N
        os_obj = OptimalStatistic(selpsrs, pta=pta, chain=ch, param_names=list(pars), gwb_name='gw_curn_pl', orfs=['hd'])
        output_NMOS = os_obj.compute_OS(N=1, return_pair_vals=False)
        A2, A2s = (output_NMOS[k] for k in ['A2','A2s'])
        SNR_per_psr.append(A2 / A2s)
        
    print("\nOK.\n")
    SNRs.append(SNR_per_psr)

# Compute median and spread per Npsr
SNRs_median = np.array([np.median(snr_list) for snr_list in SNRs])
SNRs_std    = np.array([np.std(snr_list)    for snr_list in SNRs])
SNRs_p16    = np.array([np.percentile(snr_list, 16) for snr_list in SNRs])
SNRs_p84    = np.array([np.percentile(snr_list, 84) for snr_list in SNRs])

Npsrlist_arr = np.array(Npsrlist)

#############################
### Fit for a power-law model
#############################

# Define the power-law model
def power_law(N, alpha, C):
    return C * N**alpha

# Fit in linear space (robust with median)
popt, pcov = curve_fit(power_law, Npsrlist_arr, SNRs_median, p0=[0.5, 1.0])
alpha_fit, C_fit = popt
alpha_err, C_err = np.sqrt(np.diag(pcov))
N_fine = np.linspace(Npsrlist_arr.min(), Npsrlist_arr.max(), 300)
f_fit      = power_law(N_fine, *popt)

#############################
### Same but with predicted scaling relation from Siemens et al. 2013, alpha=1. Here C is anchored to the data
#############################

power_law_fixed = lambda N, C: C * N**1.0
popt_th, _ = curve_fit(power_law_fixed, Npsrlist_arr, SNRs_median)
C_theory = popt_th[0]

### Plot
plt.figure(figsize=(14, 5))

plt.errorbar(Npsrlist_arr, SNRs_median, yerr=[SNRs_median - SNRs_p16, SNRs_p84 - SNRs_median], fmt='o', color='steelblue', capsize=4, label='Median ± 1σ (16–84%)')
plt.plot(N_fine, f_fit,   'r--', linewidth=1.5, label=fr'Fit: SNR $\propto N^{{{alpha_fit:.2f} \pm {alpha_err:.2f}}}$')
plt.plot(N_fine, C_theory * N_fine**1.0, 'g-', linewidth=1.5, label=r'Siemens et al. 2013: SNR $\propto N^{1}$')

plt.xlabel('Number of pulsars $N$', fontsize=22)
plt.ylabel('SNR', fontsize=22)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.legend(fontsize=14)
plt.grid(True, which='both', alpha=0.3)

plt.suptitle(fr'PTA SNR scaling law: SNR $\propto N^{{{alpha_fit:.2f} \pm {alpha_err:.2f}}}$', fontsize=22)
plt.tight_layout()
plt.show()

Some selected references related to Section 6

• The Nanohertz Gravitational Wave Astronomer (Taylor, 2021), very useful book about PTA, why not adding it ?
• The stochastic background: scaling laws and time to detection for pulsar timing arrays (Siemens et al., 2013), Scaling laws for the GWB significance in PTAs