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 !")

Loading modules ...

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/tmp/ipykernel_3216/2801485294.py in <cell line: 0>()
     29 import matplotlib.colors as mcolors
     30 from matplotlib.colors import LogNorm
---> 31 import corner
     32 
     33 # Astro general

ModuleNotFoundError: No module named 'corner'

---------------------------------------------------------------------------
NOTE: If your import is failing due to a missing package, you can
manually install dependencies using either !pip or !apt.

To view examples of installing some common dependencies, click the
"Open Examples" button below.
---------------------------------------------------------------------------

pip install corner

Collecting corner
  Downloading corner-2.3.0-py3-none-any.whl.metadata (2.3 kB)
Requirement already satisfied: matplotlib>=2.1 in /usr/local/lib/python3.12/dist-packages (from corner) (3.10.0)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (1.3.3)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (4.63.0)
Requirement already satisfied: kiwisolver>=1.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (1.5.0)
Requirement already satisfied: numpy>=1.23 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (2.0.2)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (26.2)
Requirement already satisfied: pillow>=8 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (11.3.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (3.3.2)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.12/dist-packages (from matplotlib>=2.1->corner) (2.9.0.post0)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.12/dist-packages (from python-dateutil>=2.7->matplotlib>=2.1->corner) (1.17.0)
Downloading corner-2.3.0-py3-none-any.whl (16 kB)
Installing collected packages: corner
Successfully installed corner-2.3.0

# 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()

<>:7: SyntaxWarning: invalid escape sequence '\z'
<>:8: SyntaxWarning: invalid escape sequence '\G'
<>:7: SyntaxWarning: invalid escape sequence '\z'
<>:8: SyntaxWarning: invalid escape sequence '\G'
/tmp/ipykernel_1151/4174761798.py:7: SyntaxWarning: invalid escape sequence '\z'
  plt.xlabel("Angle $\zeta_{ij}$ between Earth-pulsar baseline pairs [deg.]", fontsize=18)
/tmp/ipykernel_1151/4174761798.py:8: SyntaxWarning: invalid escape sequence '\G'
  plt.ylabel("Arrival time correlation $\Gamma(\zeta_{ij})$", fontsize=18)

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

<>:138: SyntaxWarning: invalid escape sequence '\z'
<>:139: SyntaxWarning: invalid escape sequence '\G'
<>:226: SyntaxWarning: invalid escape sequence '\m'
<>:250: SyntaxWarning: invalid escape sequence '\m'
<>:138: SyntaxWarning: invalid escape sequence '\z'
<>:139: SyntaxWarning: invalid escape sequence '\G'
<>:226: SyntaxWarning: invalid escape sequence '\m'
<>:250: SyntaxWarning: invalid escape sequence '\m'
/tmp/ipykernel_1151/910448248.py:138: SyntaxWarning: invalid escape sequence '\z'
  ax.set_xlabel("Angle $\zeta_{ij}$ between Earth-pulsar baseline pairs [deg.]", fontsize=18)
/tmp/ipykernel_1151/910448248.py:139: SyntaxWarning: invalid escape sequence '\G'
  ax.set_ylabel("Arrival time correlation $\Gamma(\zeta_{ij})$", fontsize=18)
/tmp/ipykernel_1151/910448248.py:226: SyntaxWarning: invalid escape sequence '\m'
  ax_prefit.set_ylabel('Prefit residuals\n[$\mu s$]', fontsize=fontsize)
/tmp/ipykernel_1151/910448248.py:250: SyntaxWarning: invalid escape sequence '\m'
  ax_post.set_ylabel('Timing residuals\n[$\mu s$]', fontsize=fontsize)

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)

Current directory: /content

# 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")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_1151/67580338.py in <cell line: 0>()
      2 coordtype = "eq"
      3 
----> 4 plot_SkyMap(psrs, coordtype, show_mw=True, psrcolor = "w")

/tmp/ipykernel_1151/910448248.py in plot_SkyMap(psrs, coordtype, psrcolor, show_mw)
     12     # Milky Way disk
     13     if show_mw:
---> 14         Plot_MW(ax, coordtype=coordtype, n_stars=5e5, angle_mask=10)
     15 
     16     # Pulsars

/tmp/ipykernel_1151/910448248.py in Plot_MW(ax, coordtype, n_stars, angle_mask)
     80 #         l_fake = rng.uniform(0, 2*np.pi, int(n_stars))[:len(b_fake)]
     81 
---> 82         coords  = SkyCoord(l=l_fake*u.rad, b=b_fake*u.rad, frame='galactic').icrs
     83         x_fake  = np.clip(coords.ra.rad - np.pi, -np.pi, np.pi)  # clip after shift
     84         y_fake  = coords.dec.rad

NameError: name 'SkyCoord' is not defined

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)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/tmp/ipykernel_1151/501929753.py in <cell line: 0>()
      3 
      4 # Reducing the number of pulsars to play with the curve
----> 5 selpsrs = np.random.choice(psrs, size=Npsrs, replace=False)
      6 
      7 # Bin size for the histogram (in deg.)

numpy/random/mtrand.pyx in numpy.random.mtrand.RandomState.choice()

ValueError: 'a' cannot be empty unless no samples are taken

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_1151/3243691878.py in <cell line: 0>()
      1 # Initiate the enterprise signal_collection object with the TimingModel
----> 2 tm = gp_signals.TimingModel(use_svd=True)
      3 
      4 ########################### White noise
      5 ## EFAC - Here we don't consider EQUAD/ECORR, just for simplicity

NameError: name 'gp_signals' is not defined

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)

1 on 3000

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_1151/1978923334.py in <cell line: 0>()
     16     xs[0,i] = g
     17     xs[1,i] = A
---> 18     xs[2,i] = pta.get_lnlikelihood(x)

NameError: name 'pta' is not defined

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

'enterprise v3.4.4,  Python v3.11.15\n==========================================================================================\n\nSignal Name                              Signal Class                   no. Parameters      \n==========================================================================================\nJ0023+0923_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0023+0923_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0023+0923_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0023+0923_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0030+0451_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0030+0451_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0030+0451_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0030+0451_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0034-0534_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0034-0534_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0034-0534_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0034-0534_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0101-6422_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0101-6422_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0101-6422_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0101-6422_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0125-2327_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0125-2327_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0125-2327_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0125-2327_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0154+1833_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0154+1833_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0154+1833_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0154+1833_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0218+4232_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0218+4232_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0218+4232_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0218+4232_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0340+4130_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0340+4130_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0340+4130_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0340+4130_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0406+3039_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0406+3039_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0406+3039_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0406+3039_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0437-4715_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0437-4715_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0437-4715_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0437-4715_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0509+0856_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0509+0856_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0509+0856_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0509+0856_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0557+1551_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0557+1551_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0557+1551_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0557+1551_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0605+3757_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0605+3757_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0605+3757_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0605+3757_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0610-2100_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0610-2100_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0610-2100_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0610-2100_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0613-0200_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0613-0200_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0613-0200_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0613-0200_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0614-3329_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0614-3329_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0614-3329_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0614-3329_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0621+1002_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0621+1002_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0621+1002_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0621+1002_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0636+5128_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0636+5128_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0636+5128_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0636+5128_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0636-3044_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0636-3044_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0636-3044_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0636-3044_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0645+5158_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0645+5158_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0645+5158_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0645+5158_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0709+0458_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0709+0458_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0709+0458_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0709+0458_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0711-6830_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0711-6830_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0711-6830_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0711-6830_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0732+2314_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0732+2314_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0732+2314_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0732+2314_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0740+6620_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0740+6620_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0740+6620_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0740+6620_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0751+1807_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0751+1807_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0751+1807_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0751+1807_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0824+0028_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0824+0028_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0824+0028_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0824+0028_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0900-3144_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0900-3144_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0900-3144_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0900-3144_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0931-1902_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0931-1902_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0931-1902_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0931-1902_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ0955-6150_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ0955-6150_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ0955-6150_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ0955-6150_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1012+5307_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1012+5307_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1012+5307_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1012+5307_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1012-4235_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1012-4235_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1012-4235_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1012-4235_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1017-7156_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1017-7156_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1017-7156_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1017-7156_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1022+1001_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1022+1001_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1022+1001_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1022+1001_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1024-0719_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1024-0719_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1024-0719_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1024-0719_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1036-8317_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1036-8317_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1036-8317_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1036-8317_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1045-4509_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1045-4509_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1045-4509_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1045-4509_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1101-6424_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1101-6424_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1101-6424_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1101-6424_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1125+7819_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1125+7819_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1125+7819_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1125+7819_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1125-5825_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1125-5825_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1125-5825_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1125-5825_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1125-6014_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1125-6014_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1125-6014_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1125-6014_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1216-6410_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1216-6410_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1216-6410_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1216-6410_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1231-1411_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1231-1411_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1231-1411_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1231-1411_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1312+0051_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1312+0051_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1312+0051_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1312+0051_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1327-0755_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1327-0755_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1327-0755_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1327-0755_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1421-4409_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1421-4409_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1421-4409_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1421-4409_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1431-5740_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1431-5740_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1431-5740_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1431-5740_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1435-6100_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1435-6100_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1435-6100_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1435-6100_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1446-4701_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1446-4701_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1446-4701_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1446-4701_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1453+1902_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1453+1902_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1453+1902_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1453+1902_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1455-3330_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1455-3330_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1455-3330_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1455-3330_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1513-2550_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1513-2550_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1513-2550_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1513-2550_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1514-4946_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1514-4946_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1514-4946_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1514-4946_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1518+4904_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1518+4904_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1518+4904_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1518+4904_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1525-5545_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1525-5545_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1525-5545_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1525-5545_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1543-5149_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1543-5149_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1543-5149_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1543-5149_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1545-4550_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1545-4550_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1545-4550_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1545-4550_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1547-5709_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1547-5709_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1547-5709_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1547-5709_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1600-3053_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1600-3053_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1600-3053_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1600-3053_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1603-7202_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1603-7202_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1603-7202_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1603-7202_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1614-2230_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1614-2230_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1614-2230_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1614-2230_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1629-6902_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1629-6902_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1629-6902_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1629-6902_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1630+3734_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1630+3734_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1630+3734_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1630+3734_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1640+2224_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1640+2224_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1640+2224_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1640+2224_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1643-1224_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1643-1224_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1643-1224_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1643-1224_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1652-4838_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1652-4838_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1652-4838_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1652-4838_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1653-2054_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1653-2054_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1653-2054_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1653-2054_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1658-5324_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1658-5324_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1658-5324_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1658-5324_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1705-1903_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1705-1903_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1705-1903_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1705-1903_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1708-3506_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1708-3506_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1708-3506_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1708-3506_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1710+4923_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1710+4923_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1710+4923_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1710+4923_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1713+0747_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1713+0747_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1713+0747_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1713+0747_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1719-1438_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1719-1438_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1719-1438_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1719-1438_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1721-2457_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1721-2457_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1721-2457_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1721-2457_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1730-2304_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1730-2304_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1730-2304_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1730-2304_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1732-5049_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1732-5049_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1732-5049_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1732-5049_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1737-0811_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1737-0811_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1737-0811_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1737-0811_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1738+0333_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1738+0333_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1738+0333_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1738+0333_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1741+1351_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1741+1351_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1741+1351_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1741+1351_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1744-1134_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1744-1134_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1744-1134_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1744-1134_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1745+1017_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1745+1017_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1745+1017_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1745+1017_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1745-0952_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1745-0952_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1745-0952_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1745-0952_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1747-4036_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1747-4036_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1747-4036_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1747-4036_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1751-2857_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1751-2857_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1751-2857_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1751-2857_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1757-5322_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1757-5322_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1757-5322_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1757-5322_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1801-1417_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1801-1417_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1801-1417_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1801-1417_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1802-2124_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1802-2124_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1802-2124_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1802-2124_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1804-2717_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1804-2717_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1804-2717_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1804-2717_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1804-2858_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1804-2858_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1804-2858_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1804-2858_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1811-2405_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1811-2405_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1811-2405_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1811-2405_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1825-0319_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1825-0319_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1825-0319_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1825-0319_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1832-0836_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1832-0836_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1832-0836_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1832-0836_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1843-1113_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1843-1113_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1843-1113_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1843-1113_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1843-1448_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1843-1448_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1843-1448_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1843-1448_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1853+1303_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1853+1303_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1853+1303_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1853+1303_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1857+0943_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1857+0943_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1857+0943_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1857+0943_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1902-5105_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1902-5105_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1902-5105_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1902-5105_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1903+0327_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1903+0327_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1903+0327_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1903+0327_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1903-7051_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1903-7051_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1903-7051_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1903-7051_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1909-3744_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1909-3744_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1909-3744_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1909-3744_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1910+1256_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1910+1256_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1910+1256_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1910+1256_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1911+1347_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1911+1347_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1911+1347_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1911+1347_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1911-1114_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1911-1114_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1911-1114_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1911-1114_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1918-0642_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1918-0642_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1918-0642_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1918-0642_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1923+2515_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1923+2515_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1923+2515_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1923+2515_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1933-6211_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1933-6211_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1933-6211_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1933-6211_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1939+2134_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1939+2134_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1939+2134_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1939+2134_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1944+0907_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1944+0907_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1944+0907_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1944+0907_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1946+3417_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1946+3417_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1946+3417_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1946+3417_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1946-5403_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1946-5403_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1946-5403_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1946-5403_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ1955+2908_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ1955+2908_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ1955+2908_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ1955+2908_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2010-1323_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2010-1323_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2010-1323_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2010-1323_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2017+0603_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2017+0603_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2017+0603_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2017+0603_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2019+2425_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2019+2425_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2019+2425_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2019+2425_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2033+1734_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2033+1734_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2033+1734_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2033+1734_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2039-3616_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2039-3616_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2039-3616_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2039-3616_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2043+1711_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2043+1711_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2043+1711_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2043+1711_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2055+3829_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2055+3829_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2055+3829_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2055+3829_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2124-3358_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2124-3358_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2124-3358_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2124-3358_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2129-5721_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2129-5721_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2129-5721_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2129-5721_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2145-0750_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2145-0750_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2145-0750_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2145-0750_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2150-0326_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2150-0326_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2150-0326_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2150-0326_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2205+6012_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2205+6012_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2205+6012_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2205+6012_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2214+3000_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2214+3000_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2214+3000_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2214+3000_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2222-0137_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2222-0137_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2222-0137_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2222-0137_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2229+2643_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2229+2643_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2229+2643_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2229+2643_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2234+0611_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2234+0611_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2234+0611_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2234+0611_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2234+0944_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2234+0944_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2234+0944_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2234+0944_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2236-5527_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2236-5527_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2236-5527_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2236-5527_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2241-5236_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2241-5236_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2241-5236_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2241-5236_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2302+4442_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2302+4442_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2302+4442_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2302+4442_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2317+1439_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2317+1439_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2317+1439_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2317+1439_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2322+2057_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2322+2057_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2322+2057_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2322+2057_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\nJ2322-2650_linear_timing_model_svd       TimingModel                    0                   \n\nparams:\n__________________________________________________________________________________________\nJ2322-2650_measurement_noise             MeasurementNoise               0                   \n\nparams:\nJ2322-2650_efac:Constant=1.0                                                              \n__________________________________________________________________________________________\nJ2322-2650_gw_curn_pl                    BasisGP                        2                   \n\nBasis shape (Ntoas x N basis functions): (200, 60)\nN selected toas: 200\n\nparams:\ngw_curn_pl_log10_A:Uniform(pmin=-18, pmax=-11)                                            \ngw_curn_pl_gamma:Uniform(pmin=0, pmax=7)                                                  \n__________________________________________________________________________________________\n==========================================================================================\nTotal params: 135\nVarying params: 2\nCommon params: 266\nFixed params: 133\nNumber of pulsars: 133\n'

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)

PTA object set for 133 pulsars, using 'ideal' noise model and 'curn_pl' CRS model.

Outdir already exists: /home/pulsar/Downloads/GWB-search/DataSim/ideal_for_IPTASW//chains/CURN_pl/

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)

PTA object set for 133 pulsars, using 'ideal' noise model and 'hd_pl' CRS model.

Outdir already exists: /home/pulsar/Downloads/GWB-search/DataSim/ideal_for_IPTASW//chains/HD_pl/

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)

PTA object set for 133 pulsars, using 'ideal' noise model and 'curn_fs' CRS model.

Outdir already exists: /home/pulsar/Downloads/GWB-search/DataSim/ideal_for_IPTASW//chains/CURN_fs/

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.

outdir = '/home/pulsar/Downloads/GWB-search/DataSim/ideal_for_IPTASW//chains/CURN_pl_precomputed'
ch_curnpl, pars_curnpl = read_all_chains(outdir, 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)

PTA object set for 133 pulsars, using 'ideal' noise model and 'curn_pl' CRS model.

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'])

NMOS Iters: 100%|██████████████████████████████████████████| 100/100 [00:15<00:00,  6.50it/s]

# 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()

/opt/build_psrsoft/workspace/micromamba/envs/IPTA_Env/lib/python3.11/site-packages/sklearn/gaussian_process/kernels.py:445: ConvergenceWarning: The optimal value found for dimension 0 of parameter k2__noise_level is close to the specified lower bound 1e-05. Decreasing the bound and calling fit again may find a better value.

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)

'\nPCOS - Takes ~15 min !\n'

"""
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()
"""

"\nNbins = 15\n\nA2 = output_PCOS['A2']\nA2s = output_PCOS['A2s']\nidx = output_PCOS['idx']\nxi = output_PCOS['xi']\nrho = output_PCOS['rho']\nsig = output_PCOS['sig']\nC = output_PCOS['C']\n\n# Plot the binned pair correlation plot. \nf, ax = defplot.create_correlation_plot(xi_OS, rho_OS, sig_OS, C_OS, A2_OS, A2s_OS, bins=Nbins, figsize=(10,6))\nax.plot(xi_range,10**(2*inj_params[f'gw_curn_pl_log10_A'])*hd_mod,'--k',label='HD prediction')\n\nax.xaxis.label.set_size(18)\nax.yaxis.label.set_size(18)\nplt.title('OS correlation vs. HD curve', fontsize=20)\n\nplt.title('PC+OS')\nplt.legend()\nplt.show()\n"

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

shifts: 100%|██████████████████████████████████████████████| 500/500 [01:16<00:00,  6.52it/s]

Measured P-value: 0.0 Upper limit: P<= 0.00199203187250996

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

scrambles: 100%|███████████████████████████████████████████| 500/500 [00:30<00:00, 16.41it/s]

Measured P-value: 0.0 Upper limit: P<= 0.00199203187250996

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

shifts: 100%|██████████████████████████████████████████████| 500/500 [01:16<00:00,  6.54it/s]

Measured P-value: 0.0 Upper limit: P<= 0.00199203187250996

# --- 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))

scrambles: 100%|███████████████████████████████████████████| 500/500 [00:30<00:00, 16.14it/s]

Measured P-value: 0.0 Upper limit: P<= 0.00199203187250996

# 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)

Run for 10 pulsars.
100.0%
OK.

Run for 20 pulsars.
100.0%
OK.

Run for 30 pulsars.
100.0%
OK.

Run for 40 pulsars.
100.0%
OK.

Run for 50 pulsars.
100.0%
OK.

Run for 60 pulsars.
100.0%
OK.

Run for 70 pulsars.
100.0%
OK.

Run for 80 pulsars.
100.0%
OK.

Run for 90 pulsars.
100.0%
OK.

Run for 100 pulsars.
100.0%
OK.

Run for 110 pulsars.
100.0%
OK.

Run for 120 pulsars.
100.0%
OK.

Run for 130 pulsars.
100.0%
OK.

# 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