# Spatial and Spectral Characterization of the Gravitational-wave Background
with the PTA Optimal Statistic

Kyle A. Gersbach<EMAIL_ADDRESS>Department of Physics &
Astronomy, Vanderbilt University,
2301 Vanderbilt Place, Nashville, TN 37235, USA  Stephen R. Taylor  Department
of Physics & Astronomy, Vanderbilt University,
2301 Vanderbilt Place, Nashville, TN 37235, USA  Patrick M. Meyers  Division
of Physics, Mathematics, and Astronomy, California Institute of Technology,
Pasadena, CA 91125, USA  Joseph D. Romano  Department of Physics and
Astronomy, University of Texas Rio Grande Valley,
One West University Boulevard, Brownsville, TX 78520, USA

###### Abstract

Pulsar timing arrays (PTAs) have made tremendous progress and are now showing
strong evidence for the gravitational-wave background (GWB). Further probing
the origin and characteristics of the GWB will require more generalized
analysis techniques. Bayesian methods are most often used but can be
computationally expensive. On the other hand, frequentist methods, like the
PTA Optimal Statistic (OS), are more computationally efficient and can produce
results that are complementary to Bayesian methods, allowing for stronger
statistical cases to be built from a confluence of different approaches. In
this work we expand the capabilities of the OS through a technique we call the
Per-Frequency Optimal Statistic (PFOS). The PFOS removes the underlying power-
law assumption inherent in previous implementations of the OS, and allows one
to estimate the GWB spectrum in a frequency-by-frequency manner. We have also
adapted a recent generalization from the OS pipeline into the PFOS, making it
capable of accurately characterizing the spectrum in the intermediate and
strong GW signal regimes using only a small fraction of the necessary
computational resources when compared with fully-correlated Bayesian methods,
while also empowering many new types of analyses not possible before. We find
that even in the strong GW signal regime, where the GWB dominates over noise
in all frequencies, the injected value of the signal lies within the
$50^{\mathrm{th}}$-percentile of the PFOS uncertainty distribution in
$41-45\%$ of simulations, remaining $3\sigma$-consistent with unbiased
estimation.

##  I Introduction

Pulsar timing arrays (PTAs) [1] exploit the precision timing of millisecond-
period pulsars to search for interpulsar-correlated timing residuals over
campaigns of years to decades. This timespan sets the period (and frequencies,
i.e., nanohertz) of GWs to which PTAs are sensitive. At such low frequencies
and with initially limited resolution, PTAs are primarily sensitive to a GW
background (GWB) [2, 3]. If this background is isotropic then the timing
correlation signature is called the Hellings & Downs (HD) curve [4], a
quadrupolar-like function of only the angular separation of pulsars. In June
2023, several major PTA collaborations announced evidence for this GWB,
including the North American Nanohertz Observatory for Gravitational waves
(NANOGrav) with $\sim 4\sigma$ significance [5], the European and Indian
collaborations (EPTA and InPTA) with joint $\sim 3\sigma$ significance [6],
and the Parkes Pulsar Timing Array (PPTA) with $\sim 2\sigma$ significance
[7]. The International Pulsar Timing Array (IPTA) conducted a comparison of
all results, which validated their consistency in terms of relative detection
sensitivity and their measurements of the GWB spectrum [8]. 111The more
recently established Chinese PTA also announced evidence, albeit only at one
frequency and using statistical techniques that are not comparable with those
used by the other PTAs [9].

As datasets grow in sensitivity to the GWB, its spatial and spectral
properties will become ever more of interest in the goal of understanding its
origin [10, 11, 12]. The primary signal hypothesis is of a large population of
supermassive black-hole binary (SMBHB) systems ($\sim
10^{8}-10^{10}M_{\odot}$) at sub-centiparsec separations, all emitting GWs
that sum incoherently to produce a stochastic timing signal [13, 5, 6]. Other
interpretations posit that part or all of the GWB is from inflationary or
early-Universe phase-transition processes, topological defects, or dark matter
[14, 15]. While many current inference techniques assume power-law models for
the GWB characteristic strain spectrum, the effects of discreteness [e.g., 16,
17] or non-GW–driven orbital evolution in a SMBHB population [18, 19], or the
details of cosmological processes [20, 21], can all induce departures from
such simple behavior. There is a growing need for detection statistics and
parameter estimation techniques that can generalize to arbitrary spectral and
spatial power distributions, and do so rapidly for large-scale programs of
simulations and explorations.

The most commonly used framework in PTA searches is that of Bayesian
inference, whereupon a generative model of pulsar timing residuals across an
array is constructed. This model is then assessed against data in a likelihood
function, and with parameter prior assumptions, the joint posterior
probability distribution of model parameters can be sampled using Markov Chain
Monte Carlo (MCMC) techniques [22, 23]. One can then construct marginals of
this joint distribution to derive parameter estimates and uncertainties. Using
a variety of approaches, ratios of the fully marginalized likelihood (i.e.,
Bayes factors) can also be calculated to act as detection statistics; in PTAs
the relevant ratio is of a steep-spectrum red process that is either HD-
correlated or uncorrelated between pulsars. The flexibility of Bayesian
modeling makes it the primary statistical approach, yet its major downside is
computational expense, with the bottleneck of the likelihood function being
covariance-matrix inversions that must be performed sequentially in MCMC
analyses that are run for millions of steps. Despite this, there are new
methods that are accelerating these steps [24, 11, 25, 26].

An alternative approach, and one which is most often used for determining
detection significance, is the PTA optimal statistic (OS) [27]. This is a
frequentist statistic originally developed to look for GWBs in ground-based
gravitational-wave detectors [28]. In a PTA context, the OS was conceived to
compute the amplitude and signal-to-noise ratio (S/N) of a GWB with a power-
law strain spectrum, usually fixed to the expected shape given by a SMBHB
population [27, 29]. Unlike the Bayesian approach, the OS focuses exclusively
on the cross-correlated signal between pulsars, leveraging the distinctiveness
of the HD curve and treating the auto-correlated GWB signal as a noise process
that adds to intrinsic noise in the pulsars. It is also a much faster
statistic to evaluate, which means that it is often the first approach carried
out to provide detection statistics, GWB amplitude estimates, reconstructed
correlation signatures, and cross-validation tests. Its speed also makes
calibration of the OS S/N with large numbers of data augmentation operations
(e.g., phase shifting [30] and sky scrambling [31]) more tractable.

Much work has recently been done to extend the OS’s applicability and
usefulness. It was originally developed under the weak-signal hypothesis,
where the GWB signal would always be buried in noise. However, this is no
longer true of the PTA GWB signal [32], meaning that the assumptions of the OS
required revision to account for pairwise correlation measurements of the GWB
no longer being independent [33, 34]. Furthermore, while the OS was originally
designed to search for just an isotropic GWB, new work has extended it to
search for anisotropy [10], as well as multiple processes among which may be
systematics that induce spatial correlations (e.g., clock errors or Solar
System ephemeris offsets) [35]. Some of these recent advances were used in
reporting NANOGrav’s evidence for the GWB [5, 36, 37]. These new methods and
extensions, along with the ones we propose in this paper, are briefly
explained in Table 1. Importantly, all of these methods can be combined
together to form even more powerful methods.

Method Name | Abbreviation | Description | Section | Citation
---|---|---|---|---
Optimal Statistic | OS |  | The first PTA optimal statistic. Tests the significance of a
---
single spatial correlation pattern. Biased in intermediate
and strong signal regimes.
II | [27, 29]
Noise-Marginalized OS | NMOS |  | An additional procedure to account for the spread in
---
parameter estimates from a Bayesian CURN search by
iteratively running the OS with many parameter estimates.
More accurate in cases with significant intrinsic red noise.
II.2 | [38]
Pair-Covariant OS | PCOS |  | An additional procedure to incorporate covariance between
---
pulsar pairs in cases with a significant GWB. Fixes OS bias
in intermediate and strong signal regimes.
II.3 | [22, 32, 33]
Multiple-Component OS | MCOS |  | A generalization of the OS which allows for simultaneous
---
estimates of multiple spatial correlation pattern amplitudes.
– | [35]
Per-Frequency OS* | PFOS |  | A generalization to the OS which allows for estimating a
---
single correlation pattern at individual frequencies.
Accounts for GWB amplitude at other frequencies.
III | –
| Narrowband-normalized
---
Per-Frequency OS*
| Narrowband
---
PFOS
| A simplification of the PFOS which estimates the GWB
---
amplitude at each frequency assuming other frequencies
don’t contribute. Biased in strong signal regimes.
III | –
Uncertainty Sampling* | UC |  | A method that more accurately describes the total GWB
---
amplitude from a NMOS analysis by accounting for the
uncertainties from each OS iteration.
II.2 | –
Table 1: A table listing the many methods and extensions of the PTA optimal
statistic. These methods can also be combined together to form a compound
method such as the Pair-Covariant Noise-Marginalized Optimal Statistic
(PC+NMOS). Names marked with a ‘*’ indicate methods which are new in this
work.

However, a core assumption of the OS still remains that the GWB has a power-
law strain spectrum, often kept at the fiducial value for a SMBHB population.
In this work, we generalize the OS beyond this assumption, allowing for our
new per-frequency optimal statistic (PFOS) to characterize the GWB spectrum
and its spatial-correlation properties using frequentist statistics. Our
method not only provides a more general frequentist framework for PTA GWB
searches, but also allows for comparisons with Bayesian measurements to
strengthen future GWB evidence, and potentially even show early signs of an
emerging individual GW source.

This paper is laid out as follows. In Section II we review the PTA optimal
statistic as it is currently implemented, as well as its recent generalization
to the arbitrary signal-strength regime. Section III introduces our new PTA
per-frequency statistic, allowing for generalized spectral estimation that can
nevertheless be subsequently fit to parameterized spectral models (e.g., power
laws, models with low-frequency turnovers, or models with single-frequency
excess due to discrete binary population behavior). We discuss the design of
simulations with which to assess the efficacy of this new statistic in Section
IV followed by the results of these assessments in Section V. Our conclusions
and discussion of future work are presented in Section VI, followed by several
appendices of supporting calculations.

##  II Traditional Optimal statistic

The optimal statistic (OS) is a cross-correlation statistic employed in GWB
searches with PTA experiments [27, 29, 39]; its name comes from the use of an
optimal filter in the cross-correlation calculation that maximizes the signal-
to-noise ratio [28, 27]. It is primarily used as a frequentist detection
statistic since it can be quickly evaluated, and is often employed before or
in parallel with more expensive Bayesian analyses. To explain these statements
in more detail, this section will focus on the OS as originally derived
(henceforth called the traditional OS for clarity) and discuss its primary use
cases. A reader familiar with the traditional OS may wish to skip to Section
III where we introduce our new techniques.

###  II.1 Creating the Traditional OS

The measurements made by PTAs come in the form of vectors of pulse times of
arrival (TOAs) for all pulsars in our PTA, $\bm{t}_{\mathrm{TOA}}$. While many
possible processes can change the expected arrival time of these pulses, each
process can be broadly grouped into either deterministic or stochastic
processes [23]. The leading-order deterministic contribution is the timing
model, whose parameters are estimated through individual pulsar analyses and
are used to predict the expected TOAs, $\bm{t}_{\mathrm{det}}$ [23]. By
subtracting the best-fit timing model from the measured TOAs, we can create a
set of timing residuals,

$\bm{\delta t}=\bm{t}_{\mathrm{TOA}}-\bm{t}_{\mathrm{det}}.$ (1)

When modeling these residuals, we must also consider stochastic processes, as
well as allowing for linearized offsets in the deterministic model fits. Hence
we model these residuals as

$\bm{\delta t}=\bm{M}\vec{\epsilon}+\bm{F}\vec{a}+\vec{n}_{w},$ (2)

where the linearized timing model design matrix, $\bm{M}$, contains partial
derivatives of TOAs with respect to each timing model parameter, allowing for
deviations from the best-fit deterministic model found in individual pulsar
analyses; $\vec{\epsilon}$ are linear timing model parameter deviations;
$\bm{F}$ is the total red-noise Fourier design matrix which contains
alternating TOA-length columns of sine and cosine elements evaluated at each
of the lowest frequencies, modeling both intrinsic red noise in a single
pulsar as well as the common red noise from the GWB; $\vec{a}$ are the
corresponding Fourier coefficients; and $\vec{n}_{\mathrm{w}}$ is an
additional white noise term intrinsic to each pulsar and mostly accounted for
by TOA measurement uncertainties.

We can write a Gaussian likelihood function over all residuals in the PTA,

$p(\bm{\delta
t}|\vec{\epsilon},\vec{a})=\frac{1}{\sqrt{\mathrm{det}(2\pi\bm{\Sigma})}}\exp{\left(-\frac{1}{2}\bm{\delta
t}^{T}\bm{\Sigma}^{-1}\bm{\delta t}\right)}.$ (3)

The residual covariance matrix, $\bm{\Sigma}$, is modeled through blocks of
individual pulsar auto- and cross-covariances,

$\bm{\Sigma}=\begin{bmatrix}\bm{P}_{1}&\bm{S}_{12}&\bm{S}_{13}&\\\
\bm{S}_{21}&\bm{P}_{2}&\bm{S}_{23}&...\\\
\bm{S}_{31}&\bm{S}_{32}&\bm{P}_{3}&\\\ &\vdots&\end{bmatrix},$ (4)

where we define sub-matrices for auto-covariance, $\bm{P}_{a}=\left<\bm{\delta
t}_{a}\bm{\delta t}_{a}^{T}\right>$, and cross-covariances,
$\bm{S}_{ab}=\left<\bm{\delta t}_{a}\bm{\delta t}_{b}^{T}\right>$ for $a\neq
b$. In these equations, the subscripts $a$ and $b$ label pulsars, such that,
e.g., $\bm{\delta t}_{a}$ is a column vector of timing residuals for pulsar
$a$. The GWB induces cross-covariances in the timing residuals of pulsars,
such that

$\bm{S}_{ab}=\Gamma_{ab}\;\bm{F}_{a}\phi\bm{F}_{b}^{T},$ (5)

where $\Gamma_{ab}$ represents the constant factor of the overlap reduction
function (ORF) between pulsars $a$ and $b$, $\bm{F}_{a}$ is pulsar $a$’s
Fourier design matrix, and $\phi$ is the diagonal Fourier-domain covariance
matrix of timing residuals induced by the GWB, represented here as the power
spectral density (PSD) of the GWB in the sine and cosine elements of each
Fourier frequency. These frequencies are modeled as multiples of a base
Fourier frequency, $f_{k}=k/T_{\mathrm{span}}$, where $T_{\mathrm{span}}$ is
the total time-span of the PTA, $k\in[1,2,\ldots,k_{\mathrm{max}}]$, and
$k_{\mathrm{max}}\ll n_{\mathrm{TOA}}$.

The traditional OS assumes that the GWB’s PSD, $S(f_{k})$, is described by a
power-law, such that for the different frequencies, $f_{k}$, of our PTA we can
write the diagonal elements of the Fourier domain covariance matrix as

$\phi_{kk^{\prime}}=\delta_{kk^{\prime}}S(f_{k})\Delta
f=\frac{A_{\mathrm{gw}}^{2}}{12\pi^{2}T_{\mathrm{span}}}\left(\frac{f_{k}}{f_{\mathrm{yr}}}\right)^{-\gamma}f_{\mathrm{yr}}^{-3},$
(6)

where $\Delta f=T_{\mathrm{span}}^{-1}$ is the width of a frequency bin,
$A_{\mathrm{gw}}$ is the amplitude of the GWB characteristic strain at a
reference frequency of $f_{\mathrm{yr}}=1/\mathrm{year}$, and $\gamma$ is the
spectral index of the power-law, which is expected to be $13/3$ for a GWB
composed of SMBHBs evolving entirely due to gravitational waves [13]. Note
that the off-diagonals of this covariance matrix are zero since we model the
GWB and all red processes as stationary. The shape of the $\phi$ matrix is
($2k_{\mathrm{max}}\times 2k_{\mathrm{max}}$) since $\phi$ includes both the
sine and cosine variances of each frequency. Also note that for ease of
notation, we suppress the $\Delta f$ factor from the PSD for the remainder of
the paper, treating the PSD to be synonymous with variance of timing residuals
at each frequency.

The ORF can be generalized to be any number of other correlation patterns like
monopole, dipole, or even a combination of those correlation patterns [35].
For simplicity, this paper uses only Hellings and Downs (HD) correlations, as
that is what is expected for an isotropic GWB [4]. The HD ORF can be written
as,

$\Gamma_{ab}=\frac{1}{2}-\frac{1-\cos{\xi_{ab}}}{4}\left[\frac{1}{2}-3\ln\left(\frac{1-\cos{\xi_{ab}}}{2}\right)\right]+\frac{1}{2}\delta_{ab},$
(7)

where $\xi_{ab}$ is the separation angle between pulsars $a$ and $b$, and
$\delta_{ab}$ is a Kronecker delta function that accounts for the pulsar term
being fully correlated only when $a=b$.

With these now defined, we can introduce the traditional OS. Chamberlin  _et
al._ [29] show that by making a first-order Taylor expansion of the PTA log-
likelihood, the traditional OS can be derived as a maximum-likelihood
estimator of $A^{2}_{\mathrm{gw}}$. Importantly, this derivation assumes that
the PTA is in the weak-signal regime, meaning that the GWB is weaker than the
sum of all other noise sources at each frequency analyzed. This assumption
will be touched on later in subsection II.3. For this paper, we opt to use the
version presented in Appendix A of Chamberlin  _et al._ [29], which presents
the traditional OS in terms of pair-wise estimators,

$\rho_{ab}=\mathcal{N}_{ab}\,\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\bm{\delta t}_{b}.$
(8)

In this equation, $\bm{\delta t}_{a}$ and $\bm{P}_{a}$ are as previously
defined, and $\hat{\bm{S}}_{ab}$ is the cross-covariance matrix between $a$
and $b$ divided by the ORF factor and amplitude of the GWB power-law, such
that $\bm{S}_{ab}=A^{2}_{\mathrm{gw}}\Gamma_{ab}\hat{\bm{S}}_{ab}$. The
normalization, $\mathcal{N}_{ab}$, is chosen such that the expectation is the
ORF-modulated power-law GWB amplitude,

$\left<\rho_{ab}\right>=\Gamma_{ab}A^{2}_{\mathrm{gw}}.$ (9)

This calculation is given in full in the Appendix A1, and leads us to our
final pair-wise estimator,

$\rho_{ab}=\frac{\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\>\hat{\bm{S}}_{ab}\>\bm{P}_{b}^{-1}\bm{\delta
t}_{b}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\>\hat{\bm{S}}_{ba}\>\bm{P}_{a}^{-1}\>\hat{\bm{S}}_{ab}\right]}.$
(10)

The uncertainty in this estimator is derived by calculating its variance. As
originally constructed and presented elsewhere, this quantity is derived under
the null hypothesis of zero inter-pulsar correlations, i.e.,
$\left<\rho_{ab}\right>_{0}=0$. Using this, and as derived in Appendix B1, the
uncertainty on our pair-wise estimator is

$\sigma_{ab,0}=\mathrm{tr}\left[\bm{P}_{b}^{-1}\>\hat{\bm{S}}_{ba}\>\bm{P}_{a}^{-1}\>\hat{\bm{S}}_{ab}\right]^{-1/2},$
(11)

where the subscript 0 indicates that this uncertainty has been derived under
the null hypothesis.

Using the set of all pair-wise estimators we then can create the full
traditional PTA optimal statistic through an inverse noise-weighted sum of the
pair-wise estimators [29],

$\hat{A}^{2}_{\mathrm{gw}}=\frac{\sum_{a<b}\>\rho_{ab}\Gamma_{ab}/\sigma^{2}_{ab,0}}{\sum_{a<b}\>\Gamma^{2}_{ab}/\sigma^{2}_{ab,0}}.$
(12)

The associated uncertainty in this estimator is then calculated as the square
root of the variance

$\sigma_{\hat{A}^{2},0}=\left(\sum_{a<b}\frac{\Gamma_{ab}^{2}}{\sigma_{ab,0}^{2}}\right)^{-1/2}.$
(13)

Most of the matrices involved in these calculations are large and require
correspondingly computationally expensive operations, such as the inverse
matrix $\bm{P}^{-1}_{a}$ which has dimension equal to the squared number of
TOAs for pulsar $a$. To help reduce the computational cost, we use rank-
reduction strategies as detailed in Appendix A of Pol  _et al._ [10]. They
define two matrix products,

$\bm{X}_{a}=\bm{F}_{a}^{T}\bm{P}_{a}^{-1}\bm{\delta t}_{a},$ (14)
$\bm{Z}_{a}=\bm{F}_{a}^{T}\bm{P}_{a}^{-1}\bm{F}_{a}.$ (15)

Since these matrices now represent quantities in the Fourier domain, we will
also need to decompose our residual cross-covariance matrices with the
substitution, $\hat{\bm{S}}_{ab}=\bm{F}_{a}\hat{\phi}\bm{F}_{b}^{T}$, where
$\phi=A^{2}_{\mathrm{gw}}\hat{\phi}$. Using these definitions, we can rewrite
Equation 10 and Equation 11 in a more computationally efficient form,

$\rho_{ab}=\frac{\bm{X}_{a}^{T}\>\hat{\phi}\>\bm{X}_{b}}{\mathrm{tr}\left[\bm{Z}_{b}\>\hat{\phi}\>\bm{Z}_{a}\>\hat{\phi}\right]},$
(16)
$\sigma_{ab,0}=\mathrm{tr}\left[\bm{Z}_{b}\>\hat{\phi}\>\bm{Z}_{a}\>\hat{\phi}\right]^{-1/2}.$
(17)

It is important to understand the many versions of $\phi$ that will appear in
this paper. On its own $\phi$ represents the covariance matrix of the Fourier
coefficients of the timing residuals induced by the GWB, which we take to be
diagonal. $\hat{\phi}$ factors out the GWB amplitude making it a
representation of the unit-amplitude power-law spectral template of the GWB.

To further generalize the traditional OS for future expansions, we frame
Equation 12 as the solution to a chi-squared minimization scheme [29], with

$\chi^{2}=\left(\vec{\rho}-A^{2}_{\mathrm{gw}}\vec{\Gamma}\right)^{T}\bm{C}_{0}^{-1}\left(\vec{\rho}-A^{2}_{\mathrm{gw}}\vec{\Gamma}\right),$
(18)

where we now represent all pulsar-pair estimators, $\rho_{ab}$, and
corresponding ORF factors, $\Gamma_{ab}$, as column vectors of length equal to
the number of distinct pulsar pairs in the PTA. $\bm{C}_{0}$ is the pulsar-
pair covariance matrix assuming the null-hypothesis (indicated by the 0
subscript), whose elements are defined as

$\bm{C}_{0}=\delta_{ac}\,\delta_{bd}\;\bm{C}_{ab,cd,0}=\delta_{ac}\,\delta_{bd}\;\sigma^{2}_{ab,0}.$
(19)

In the traditional OS derivation, this matrix is diagonal and trivial to
invert. A more detailed discussion about this covariance matrix is given in
subsection II.3.

Minimizing this chi-square results in a square amplitude estimator and a more
generalizable form of Equation 12,

$\hat{A}^{2}_{\mathrm{gw}}=\left(\vec{\Gamma}^{T}\bm{C}_{0}^{-1}\vec{\Gamma}\right)^{-1}\vec{\Gamma}^{T}\bm{C}_{0}^{-1}\vec{\rho}.$
(20)

We then compute the uncertainty in this estimator under the null hypothesis by
taking the square root of its variance,

$\sigma_{\hat{A}^{2},0}=\left(\vec{\Gamma}^{T}\bm{C}_{0}^{-1}\vec{\Gamma}\right)^{-1/2}.$
(21)

Since our method modifies the traditional OS, we will frequently refer back to
this section to explain what we are changing.

###  II.2 Uncertainty Sampling

As constructed, the traditional OS requires estimates of every pulsar’s noise
auto-covariance, $\bm{P}_{a}$. For these, all deterministic signals, intrinsic
pulsar noise, and common pulsar noise for each pulsar must be estimated in
some way [27, 29, 38]. These pulsar noise estimates often come from a quick
PTA Bayesian analysis. Specifically, we get pulsar noise and common process
estimates from a PTA common uncorrelated red-noise (CURN) search, which
samples from a posterior with a similar likelihood to the one detailed in
Equation 3, except that we set the cross-correlated terms $\bm{S}_{ab}$ of the
covariance matrix in Equation 4 to zero. We opt for a CURN search as it is
roughly two orders-of-magnitude faster to analyze than a fully inter-pulsar
correlated search, and provides reasonably accurate estimates for pulsar noise
parameters in the presence of such a GWB-like process.

Naively, we could use the maximum likelihood noise and CURN parameters from
this MCMC analysis. This however does not account for uncertainty in these
parameters as represented by the posterior spread, and may result in poor
performance of the estimator. Vigeland  _et al._ [38] present a more robust
methodology dubbed the Noise Marginalized Optimal Statistic (NMOS), in which
they marginalize over the spread in parameter posteriors by computing the
traditional OS over many draws from the MCMC chain. This results in a
distribution of both $\hat{A}^{2}_{\mathrm{gw}}$ and $\sigma_{\hat{A}^{2},0}$

The distribution of OS quantities given by the NMOS can be expressed as a
Monte Carlo integral over the CURN posterior samples,

$\displaystyle p(X|d)$ $\displaystyle=\int\delta(X-X_{\eta})p(\eta|d)\,d\eta$
(22) $\displaystyle\approx\frac{1}{N}\sum_{i}^{N}\delta(X-X_{\eta_{i}}),$

where $d$ represents the PTA data, $\eta$ is a vector of noise parameters for
the PTA, $p(\eta|d)$ is the joint posterior distribution of these PTA noise
parameters from the CURN search, $\\{\eta_{i}\\}_{N}$ are a set of $N$ random
draws of these PTA noise parameters from $p(\eta|d)$, and
$X\in\\{\hat{A}^{2}_{\mathrm{gw}},\sigma_{\hat{A}^{2},0}\\}$ is either the
estimated $\hat{A}^{2}_{\mathrm{gw}}$ or its uncertainty
$\sigma_{\hat{A}^{2},0}$ from the OS calculation, and $X_{\eta}$ is notation
to indicate when these quantities have been computed for a specific choice of
$\eta$.

While the NMOS solves some problems with the traditional OS, it also
introduces a new issue. Each vector of parameter estimates, $\eta$, creates
its own $\hat{A}^{2}_{\mathrm{gw}}$ with Gaussian standard-deviation
uncertainty $\sigma_{\hat{A}^{2},0}$. In previous analyses the noise-
marginalized distribution of $\hat{A}^{2}_{\mathrm{gw}}$ is shown while
neglecting the $\sigma_{\hat{A}^{2},0}$ values [40, 38]. More recent
publications have omitted OS parameter estimation entirely [5], instead
focusing on it as a detection statistic in which $\sigma_{\hat{A}^{2},0}$ is
only used to compute the SNR value for a particular $\eta$. As shown in Figure
1, plotting only the distribution of $\hat{A}^{2}_{\mathrm{gw}}$ values and
ignoring the underlying uncertainty in the estimator can cause narrow peaks in
the estimated distribution, making the NMOS amplitude estimate seem more
constrained than it really is.

Figure 1: A normalized histogram of both the NMOS $\hat{A}^{2}_{\mathrm{gw}}$
estimator and the uncertainty-sampled estimator distributions for a realistic
simulation with an injected GWB $A_{\mathrm{gw}}=10^{-14}$. Ignoring the
contributions of the NMOS’s uncertainty results in a peaked distribution which
excludes the injected GWB amplitude.

To combat this problem, we propose an additional step in all OS pipelines,
which we dub uncertainty sampling. This step replaces each
$\hat{A}^{2}_{\mathrm{gw}}$ estimate with many random draws from a Gaussian
with mean, $\hat{A}^{2}_{\mathrm{gw}}$, and variance,
$\sigma^{2}_{\hat{A}^{2},0}$, implied by a given vector of pulsar parameter
estimates, $\eta$. 222This is akin to using a Gaussian kernel density
estimator where the width of each $\hat{A}^{2}$ point from NMOS sampling is
determined by the point’s corresponding uncertainty, $\sigma_{\hat{A}^{2},0}$.
Importantly, while a Gaussian distribution is an appropriate approximation in
usual circumstances, the true distribution is a generalized-$\chi^{2}$ [41].
We can frame this as an additional Monte-Carlo integral to the one already
performed during the NMOS calculation, such that we can now approximate the
full distribution on the GWB amplitude from the OS as

$\displaystyle p(A^{2}_{\mathrm{gw}}|d)$ $\displaystyle=\int
p(A^{2}_{\mathrm{gw}}|\\{\hat{A}^{2}_{\mathrm{gw}},\sigma_{\hat{A}^{2},0}\\}_{\eta})\,p(\eta|d)\,d\eta$
(23)
$\displaystyle\approx\frac{1}{N}\sum_{i}^{N}p(A^{2}_{\mathrm{gw}}|\\{\hat{A}^{2}_{\mathrm{gw}},\sigma_{\hat{A}^{2},0}\\}_{\eta_{i}})$
$\displaystyle\approx\frac{1}{N\times
M}\sum_{i}^{N}\sum_{j}^{M}\delta(A^{2}_{\mathrm{gw}}-Y_{\eta_{i},j}),$

where $\\{Y_{\eta_{i},j}\\}_{M}\sim
N(\hat{A}^{2}_{\mathrm{gw}},\sigma_{\hat{A}^{2},0})_{\eta_{i}}$ are $M$
samples from the uncertainty distribution of the OS (assumed to be Gaussian)
with mean and standard deviation given by $\hat{A}^{2}_{\mathrm{gw}}$ and
$\sigma_{\hat{A}^{2},0}$, which depends on a given set of PTA noise parameters
$\eta_{i}$.

Using our realistic simulations, we find that $10^{3}$ NMOS samples, and using
uncertainty sampling of 100 draws per NMOS sample, was sufficient in
stabilizing the final estimator distribution. However, the exact numbers for
these vary depending on the data set, and the stability of the distribution
should always be checked. In this case, we did so by comparing the full
uncertainty sampled $A^{2}_{\mathrm{gw}}$ distribution with one made with a
smaller portion of those of the same distribution, e.g., 90% of the
$A^{2}_{\mathrm{gw}}$ distribution. If these two distributions closely match,
the solution is likely stable.

With uncertainty sampling, we now have a distribution of possible
$A^{2}_{\mathrm{gw}}$ values that accounts for both the spread from the PTA
noise parameter estimation and the spread in the GWB amplitude estimation.
Shown in orange in Figure 1, we see that the widening of the distribution now
covers the injected value of the simulation. All subsequent analysis in this
paper will be performed using uncertainty sampling, unless otherwise noted.

###  II.3 Parameter estimation

As so far discussed, the traditional OS is designed to work for PTA data in
the weak-signal regime, defined such that the contribution of the GWB to the
total noise is negligible in all measured frequencies [39]. This assumption
leads to using a null-hypothesis diagonal pulsar-pair covariance matrix
$\bm{C}_{0}$ in Equation 20. When evaluating the traditional OS in the weak-
signal regime, the amplitude estimator is unbiased [38].

The problem, however, is that some PTA datasets are no longer in the weak-
signal regime, as the GWB signal is becoming the dominant contribution to the
total red noise in the lowest frequencies of many pulsars [5]. Revisiting our
covariance matrix calculation,

$\bm{C}=\bm{C}_{ab,cd}=\left<\rho_{ab}\,\rho_{cd}\right>-\left<\rho_{ab}\right>\left<\rho_{cd}\right>,$
(24)

in larger signal regimes, the correlated GWB signals in the $\bm{S}_{ab}$
terms we originally assumed to be zero in Appendix B1 become more significant.
Ignoring these GWB contributions results in an estimator with a variance that
is too small [33]. This was also influential in the recent NANOGrav 15-year
GWB analysis, in which it was found that neglecting these terms results in
$20-40\%$ smaller uncertainties [5]. The solution to this is a new technique
of including the pulsar-pair covariance [32, 33, 22], which involves a
replacement of the approximate diagonal $\bm{C}_{0}$ with a its exact dense
form $\bm{C}$ in Equation 20 and Equation 21. While mathematically
straightforward, calculating all elements of this covariance matrix can be
computationally expensive without simplifying strategies. We use the
covariance matrices described in Johnson  _et al._ [22] and apply the same
rank-reducing strategies found in Pol  _et al._ [10] to drastically reduce the
computational resource needs. We detail the calculation of the pulsar-pair
covariance matrix elements in Appendix C.

A key step in calculating the pair covariance matrix is the inclusion of an
assumed GWB amplitude (see equations 81, 82, and 83). This presents a somewhat
circular problem, as we must assume a GWB amplitude to then measure the GWB
amplitude. Allen and Romano [33] proposes a treatment with the assumed
amplitude and estimated amplitude as separate quantities, and use
$\chi^{2}$-fitting to find regions in which the assumed and measured
amplitudes agree. Attempting this solution however, we find that the inversion
of the pulsar-pair covariance matrix can become numerically unstable in some
cases, making it impossible to find regions in which the assumed and measured
amplitudes agreed. The scenarios and regimes when this occurs have not yet
been fully explored and are needed to implement this proposed solution.

We instead decide to sidestep this issue by implementing the same strategies
that the traditional OS does. Rather than varying the assumed GWB amplitude,
we set this value to the value found in the preliminary Bayesian CURN search.
This scheme is consistent with the way the current implementation assumes the
values of the auto-correlations in the covariance matrix. In our testing with
simulated data, this scheme resulted in estimated GWB amplitudes that are
consistent with the injected amplitudes, even in strong signal simulations.

We compare the performance of the traditional OS with and without pulsar-pair
covariance using Percentile-Percentile (P-P) plots. We created 100 realistic
simulated datasets with a GWB power-law amplitude of
$A^{2}_{\mathrm{gw}}=10^{-14}$ (additional details about these simulations are
given in Section IV). This amplitude is intentionally chosen to be far
stronger than is seen in current datasets (NANOGrav has a median amplitude of
$2.4\times 10^{-15}$ [5]). This ensures that the GWB dominates in all 10
frequencies used in the analysis in order to specifically strain the abilities
of the traditional OS with and without pulsar-pair covariance.

For each simulation we used the NMOS by directly providing it with the
injected parameter values for the simulation, effectively assuming that our
CURN analysis parameter recovery is perfect. We then used uncertainty sampling
to create a distribution of $A^{2}_{\mathrm{gw}}$ values, after which we found
the percentile at which the injected amplitude lies. In Figure 2, the P-P plot
shows the cumulative distribution function of the injected-amplitude
percentiles from each simulation for the traditional OS with and without
pulsar-pair covariance modeling. We clearly see the impact of incorporating
pulsar-pair covariance, as this new estimator remains within the $95\%$
confidence region throughout, while without it the estimator shows a sideways
‘S’ shape, which is characteristic of an underestimated uncertainty [11].

Figure 2: P-P plot comparing the traditional OS with and without pulsar-pair
covariance incorporated. The pair covariance calculation assumes injected
parameter values. Each line represents the same 100 realistic simulations with
an injected GWB amplitude of $10^{-14}$. The shaded regions indicate the
1$\sigma$, 2$\sigma$, and 3$\sigma$ confidence contours for an unbiased
estimator.

###  II.4 Detection statistic

The second major use case for the optimal statistic is to gauge the
significance of inter-pulsar correlations induced by a GWB signal. We can
compute a signal-to-noise ratio (SNR) by taking the ratio of our estimator
$\hat{A}_{\mathrm{gw}}^{2}$ and its uncertainty calculated assuming the null
hypothesis:

$\mathrm{SNR}=\frac{\hat{A}^{2}_{\mathrm{gw}}}{\sigma_{\hat{A}^{2},0}}.$ (25)

The traditional OS SNR is often used as an informal proxy for the significance
of cross-correlations in terms of standard deviations away from zero; however
the null hypothesis distribution of the traditional OS is formally non-
Gaussian, and therefore the significance must be calibrated empirically or
through specification of the true generalized-$\chi^{2}$ null hypothesis
distribution [41].

Expanding this scheme with the NMOS is trivial, as each vector of PTA
parameter estimates results in an SNR. Over many draws of PTA parameter
estimates, we can construct a distribution of SNRs where the median more
closely follows the injected-value SNR [38]. Unlike in the case of parameter
estimation, the SNR calculation works well without the inclusion of pulsar-
pair covariance, since it will be calibrated anyway through data augmentation
and bootstrapping techniques. The two current methods of calibration are
numerical estimation through schemes like sky scrambles and phase shifts [42,
30], or analytic derivation by calculating the generalized-$\chi^{2}$
distribution that describes the true shape of the null distribution [41].

##  III A New Per-Frequency Optimal Statistic

The primary goal of this work is to further generalize the traditional OS to
work for any GWB spectral shape, not just power-law behavior. We call the
method presented here the Per-Frequency Optimal Statistic (PFOS). The
traditional OS was constructed to measure the correlated amplitude and SNR of
a power-law spectrum GWB signal, using a unit-amplitude power-law spectral
template for the GWB, $\hat{\phi}=\phi/A^{2}_{\mathrm{gw}}$, where $\phi$ is
represents the power-law PSD of the GWB in Equation 6. This broadband
estimation condenses the correlated information from each frequency into a
single, more sensitive measurement [27]. Recalling the rank-reduced forms in
Equation 16 and Equation 17 for the cross-correlation estimator, we can change
the spectrum which that estimator assumes by simply swapping $\hat{\phi}$ with
another spectral template.

Here we change the spectral template to be entirely agnostic to the spectral
shape. Since the traditional OS already deals in discrete Fourier frequencies
in both the $\bm{X}_{a}$ and $\bm{Z}_{a}$ terms, we replace the unit-amplitude
Fourier domain covariance matrix, $\hat{\phi}$, with a frequency-selector,
$\tilde{\phi}(f_{k})$, which multiplies and adds the adjacent sine and cosine
components of the two pulsars at the specified frequency. This frequency-
selector will still have the same shape as $\phi$, ($2k_{\mathrm{max}}\times
2k_{\mathrm{max}}$) elements and still be diagonal, however, all of the
elements of this diagonal matrix will be zero except for those which
correspond to frequency $f_{k}$. This is best shown through examples,

$\displaystyle\tilde{\phi}(f_{1})$
$\displaystyle=\mathrm{diag}(1,1,0,0,...,0,0),$ (26)
$\displaystyle\tilde{\phi}(f_{2})$
$\displaystyle=\mathrm{diag}(0,0,1,1,...,0,0),$
$\displaystyle\tilde{\phi}(f_{k_{\mathrm{max}}})$
$\displaystyle=\mathrm{diag}(0,0,0,0,...,1,1),$

where $f_{k}$ represents the frequency to be analyzed (which must be a
multiple of the lowest PTA frequency). However, swapping the spectral template
from $\hat{\phi}$ to $\tilde{\phi}(f_{k})$ does introduce a new assumption
that measurements at each frequency are independent of one another; we test
this assumption in subsection V.3.

We now construct a new pair-wise cross-correlation estimator at the specified
frequency,

$\rho_{ab}(f_{k})=\mathcal{N}_{ab}(f_{k})\,\bm{X}^{T}_{a}\tilde{\phi}(f_{k})\bm{X}_{b},$
(27)

where $\mathcal{N}_{ab}(f_{k})$ is now a frequency-dependent normalization,
defined such that that the expectation is the ORF-independent GWB PSD,
$S(f_{k})$, at frequency $f_{k}$,
$\left<\rho_{ab}(f_{k})\right>=\Gamma_{ab}\,S(f_{k})$. This calculation is
presented in full in Appendix A2. The resulting estimator is then

$\rho_{ab}(f_{k})=\frac{\bm{X}^{T}_{a}\tilde{\phi}(f_{k})\bm{X}_{b}}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\Phi(f_{k})\right]},$
(28)

and the variance in this estimator—calculated in Appendix B2—is

$\sigma_{ab,0}(f_{k})^{2}=\frac{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\tilde{\phi}(f_{k})\right]}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\;\bm{Z}_{b}\Phi(f_{k})\right]^{2}}.$
(29)

In both of these equations, we introduce a new matrix quantity
$\Phi(f_{k})=\phi/S(f_{k})$. This represents the shape of the Fourier-domain
covariance matrix with unit-value at frequency $f_{k}$. This template will be
unique for each frequency measured and allows us to measure each frequency one
at a time. Similar to the traditional OS, the uncertainties in Equation 29 are
calculated under the null hypothesis where $\left<\rho_{ab}(f_{k})\right>=0$.
Since the normalization term introduces some broadband dependence to the PFOS,
we refer to this as the broadband-normalized PFOS (or just the PFOS).

Including $\Phi(f_{k})$ may seem counterintuitive, as it requires information
on the spectrum we are currently trying to estimate. However, this is similar
to the problem of incorporating $A^{2}_{\mathrm{gw}}$ as an assumed amplitude
when using pulsar-pair covariance with the traditional OS. We employ the same
strategy that was discussed in subsection II.3 where we use the CURN estimates
of the GWB as the assumed spectral shape. As we will soon show, this scheme
produces estimates that are consistent with injection values in our
simulations while also allowing us to fit frequencies one at a time, and
keeping our covariance matrices smaller and more stable. Details about effects
resulting from the choice of a CURN model can be found in subsection V.2 and
Section VI.

We can now estimate the GWB PSD, $S(f_{k})$, using pair-wise correlation
information from the full PTA, such that

$\hat{S}(f_{k})=\left[\vec{\Gamma}^{T}\bm{C}(f_{k})^{-1}\vec{\Gamma}\right]^{-1}\vec{\Gamma}^{T}\bm{C}(f_{k})^{-1}\vec{\rho(f_{k})},$
(30)

and the corresponding uncertainty at each frequency can be written as,

$\sigma_{\hat{S}(f_{k}),{0}}=\left[\vec{\Gamma}^{T}\bm{C}(f_{k})^{-1}\vec{\Gamma}\right]^{-1/2}.$
(31)

In these equations, we again express our pair-wise estimators as a vector of
all pairs for a particular frequency $f_{k}$ with a corresponding vector of
ORF values. The pulsar-pair covariance matrix for a particular frequency is
then written as,

$\bm{C}(f_{k})=\left<\rho_{ab}(f_{k})\rho_{cd}(f_{k})\right>-\left<\rho_{ab}(f_{k})\right>\left<\rho_{cd}(f_{k})\right>.$
(32)

Similar to the traditional OS, we can choose to incorporate pulsar-pair
covariance or not. If we do not then the covariance matrix is again diagonal,
leading to
$\bm{C}_{ab,ab}(f_{k})=\delta_{ac}\delta_{bd}\,\sigma_{ab,0}(f_{k})^{2}$.
Including pair covariance requires calculating all of the elements from the
$\bm{C}_{ab,cd}(f_{k})$ matrix, which can be found in Appendix C.

We can make an additional simplification to the PFOS method if we assume that
all the GW power comes from only the single frequency being analyzed, i.e.,
there is no measurable GWB signal in frequencies other than $f_{k}$. In this
scheme, our “true” unit-amplitude spectral shape changes to be a Kronecker
delta function at the desired frequency, such that
$\Phi(f_{k})=\tilde{\phi}(f_{k})$ in Equation 28 and Equation 29. 333This
somewhat analogous to a GWB “radiometer” anisotropy search where we use a sky-
position delta function for gravitational wave power. The narrowband-
normalized PFOS instead uses a frequency specific delta function [43, 44]. We
will refer to this form as the narrowband-normalized PFOS.

The narrowband-normalized PFOS has similar cross-correlation estimators and
associated variances given by

$\rho_{ab,\mathrm{narrow}}(f_{k})=\frac{\bm{X}^{T}_{a}\tilde{\phi}(f_{k})\bm{X}_{b}}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\tilde{\phi}(f_{k})\right]},$
(33)
$\sigma_{ab,0,\mathrm{narrow}}(f_{k})^{2}=\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\tilde{\phi}(f_{k})\right]^{-1}.$
(34)

We can then still use Equation 30 and Equation 31 to estimate $S(f_{k})$.

While the narrowband-normalized PFOS removes the need for an initial estimate
of the spectral shape and simplifies the equations, it leads to slightly
biased results in cases where there is a detectable signal in other
frequencies. Based on our testing, this simplification does work well in cases
where individual pulsars have similar observing cadence and length. However,
as we will show in subsection V.2, using more realistic data shows worse
performance than the broadband-normalized PFOS. Because of this, we recommend
using the broadband-normalized PFOS unless the data set is within the weak-
signal regime.

Finally, it is also possible to use the different flavors of PFOS to calculate
SNRs at each frequency. Similar to the traditional OS, the PFOS’s null
hypothesis is the absence inter-pulsar correlations. However, the null
hypothesis should be interpreted differently depending on whether the
broadband-normalized or narrowband-normalized PFOS is used. The broadband-
normalized PFOS implies a null hypothesis in which one assumes that no
measurable correlations exist in a particular frequency bin, but permits them
to be present in other bins. This is less biased in cases where there are
significant correlations in other frequency bins, making it possible to
isolate the significance of a single frequency bin in a detectable spectrum.
By contrast, the narrowband-normalized PFOS implies a null hypothesis in which
one instead assumes that no frequency bins whatsoever have measurable
correlations. Hence the narrowband-normalized PFOS null hypothesis is most
useful in cases where there are not significant correlations present at any
frequency, and as such we are trying to estimate the significance in a single
frequency bin. These interpretations and the potential for other
interpretations of the null hypothesis should be deliberated upon on a case-
by-case basis, and should always be calibrated through empirical techniques
[30, 42].

##  IV Simulation Design

We design sets of simulated datasets in order to test the efficacy of the
PFOS. These simulated data sets are all created using the toasim methods found
within the libstempo python package [45]. Each simulation is based on the
pulsars of NANOGrav’s 12.5 year data set [40], including their observational
schedules and noise characteristics. We use many of the same simulation
methods found in Pol  _et al._ [46]. Briefly, we use all $45$ pulsars which
have timing baselines longer than 3 years, where the full full PTA baseline is
$T_{\mathrm{span}}=12.9$ years. We then condense the many near-simultaneous
TOA measurements resulting from narrow-band timing into single epoch-averaged
TOAs. This reduces the original $410,064$ TOA measurements into just $6,244$
with a median TOA uncertainty shrinking from $2.0\times 10^{-6}\,\mathrm{s}$
to $2.9\times 10^{-7}\,\mathrm{s}$. For each simulation, we opted to search in
the 10 lowest Fourier bins of our PTA, corresponding to
$f_{k}=n/T_{\mathrm{span}}$ where $n\in[1,10]$. These changes simplify our
simulated datasets to make injections and analyses faster, while still
retaining realistic data quality features such as multi-year gaps in TOAs for
some pulsars, observations from multiple telescopes, and differences in pulsar
timing baselines.

For each pulsar, we simulate timing data that includes white noise, intrinsic
red noise, and a GWB signal that is correlated among pulsars. The white noise
in each pulsar is drawn according to the TOA uncertainties. The intrinsic red
noise is generated as a time-domain realization of a process with a power-law
power spectrum, with power-law parameters matching the maximum likelihood
values found in Arzoumanian  _et al._ [40]. For our analyses, we opt to fix
these intrinsic red noise power-law parameters in our analyses. This
eliminates the need to simultaneously search for these parameters, vastly
improving the speed and sampling performance for the initial Bayesian CURN
search upon which the OS and its variants rely. We verified in a small number
of simulation studies that the behavior when we fixed versus varied the
intrinsic red-noise parameters was similar. The major difference was the
larger uncertainty in the final PFOS-derived quantities compared to when
intrinsic red noise was held fixed.

Finally, to each simulated dataset we add a GWB signal that is correlated
across all pulsars according to the Hellings & Downs curve, with a power-law
power spectrum that can be parameterized by the amplitude and spectral index
found in Equation 6. All of our simulations have a spectral index of
$\gamma=13/3$. We split our $300$ total simulations into three groups of
$100$, where each group has a different injected power-law GWB amplitude:
$A_{\mathrm{gw}}\in[10^{-16},10^{-15},10^{-14}]$. These correspond to a median
NMOS SNR of $0.02$, $1.91$, and $7.43$, respectively, creating the conditions
of the weak, intermediate, and strong signal regimes. These regimes are
described by Siemens  _et al._ [39] in terms of the ratio of the GWB PSD to
those of other noise sources. In the weak signal regime, all GW frequencies
analyzed are noise dominated; the intermediate regime has some fraction of
frequencies dominated by the GWB; and the strong regime has all analyzed
frequencies dominated by the GWB.

We perform initial Bayesian analyses of each simulation assuming a power-law
CURN approximate model for the GWB signal (i.e., no inter-pulsar correlations
modeled) with two parameters corresponding to power-law spectral index,
$\gamma$, and the amplitude at $f_{\mathrm{yr}}$. We find that this model
produces PTA parameter estimates that lead to better and more consistent
performance when passed to the PFOS framework than when we use a more agnostic
free-spectrum CURN search, as variations in spectral shape between samples is
far more consistent. See Section VI for further discussion on the choice of
model for the initial pilot Bayesian characterization of the PTA noise.

##  V Results

In this section we demonstrate the accuracy of the broadband-normalized
PFOS—henceforth referred to as the PFOS—in our different simulated dataset
scenarios. This section is split into four parts. First, we analyze the
capabilities of the PFOS on a single case-study data set, both for measuring
the PSD of the GWB and for reconstructing the spatial correlations at each
frequency. Second, using the suite of simulations discussed in Section IV, we
assess the accuracy of PFOS-derived PSD measurements at each frequency. Third,
we test broadband power-law spectral fitting with the PFOS products in order
to compare the performance to that of the traditional OS, in which a broadband
power-law model is assumed from the beginning. Finally, we show results from a
simulation in which an isotropic stochastic GWB signal as well as a prominent
supermassive black-hole binary GW signal have been injected, demonstrating the
prospects for excess-power–based searches for individual sources with the
PFOS.

###  V.1 Case study on a single data set

For our first test we analyze a single simulation with an injected GWB
amplitude of $A_{\mathrm{gw}}=10^{-14}$ at $f={\rm yr}^{-1}$, where the PFOS
calculation is performed with the inclusion of pulsar-pair covariance, and in
which we supply the injected noise and signal hyper-parameters. This gives the
reconstructed spectrum of timing variance induced by the GWB shown in Figure
3. All but one frequency contains the injected GWB spectrum within their
$1\sigma$, where the eighth frequency contains the injected value just outside
this at the $1.03\sigma$ level.

Figure 3: The recovered pulsar-pair correlated GWB spectrum from the
(broadband-normalized) PFOS given the injected simulation parameters with
point estimates on $S(f_{k})$ and associated 1-sigma error bars. The black
dashed line indicates the injected GWB PSD.

As with the traditional OS, the PFOS allows one to construct a binned
estimator to visualize the cross-correlated power as a function of pulsar
angular separations. However, the PFOS’s ability to separate correlations
between frequencies allows us to create these binned estimators at each
frequency, enabling more precise examinations of the structure of cross-
correlations in our data. We use methods detailed in Allen and Romano [33] and
Johnson  _et al._ [22] to create a binned estimator of cross-correlations with
corresponding uncertainty,

$\displaystyle\rho_{\mathrm{opt},i}(f_{k})$
$\displaystyle=\Gamma_{\xi_{i}}\,\left(\vec{\Gamma_{i}}^{T}[\bm{C}_{i}(f_{k})]^{-1}\vec{\Gamma_{i}}\right)^{-1}$
(35) $\displaystyle\hskip
36.135pt\times\vec{\Gamma_{i}}^{T}[\bm{C}_{i}(f_{k})]^{-1}\vec{\rho_{i}(f_{k})},$
$\sigma_{\mathrm{opt},i}(f_{k})=\left(\vec{\Gamma_{i}}^{T}[\bm{C}_{i}(f_{k})]^{-1}\vec{\Gamma}_{i}\right)^{-1/2}.$
(36)

In these equations, $i$ is an index selector which represents the subset of
pulsar pairs for a particular angular separation bin, $\vec{\Gamma}_{i}$ is a
vector of ORF values for each pulsar pair within bin $i$,
$\vec{\rho}_{i}(f_{k})$ is a vector of pulsar-pair correlated PSD estimators
within bin $i$, $\bm{C}_{i}(f_{k})$ is the pulsar-pair covariance matrix for
frequency $f_{k}$ for pulsar-pairs within angular separation bin $i$, and
$\Gamma_{\xi_{i}}$ is the ORF of the average angular separation of the pulsar
pairs in bin $i$.

The resulting set of estimators for the first GW frequency ($f_{1}=2.45$ nHz)
for our simulation using the PFOS is shown in Figure 4. Note that—despite
being suppressed in this visualization—the bin-wise estimators are covariant.
Additionally, the correlated PSD estimate and corresponding uncertainty
represented in red are computed using the full covariance between each pulsar
pair. It would be straightforward to incorporate the newly developed Multiple
Component Optimal Statistic (MCOS) to search for many correlated processes
simultaneously at specific frequencies [35].

Figure 4: The recovered binned pulsar pair correlated power for the first GW
frequency ($f_{1}=2.45$ nHz) for the (broadband-normalized) PFOS. The dashed
black line represents the injected GWB signal. The solid red represents the
best fit signal. The shaded red region represents the 1-sigma uncertainty on
the best fit signal. The blue points represent the binned cross-correlation
measurements, where the bins were chosen so that there are the same number of
pulsar pairs in each bin.

Once we incorporate noise marginalization and uncertainty sampling, we create
a distribution of $\hat{S}(f_{k})$ that will more accurately represent the
total uncertainty in the estimator. Figure 5 shows this in a violin plot with
the medians represented as horizontal lines and extrema represented by caps.
Since we are looking at a full distribution of estimators rather than
$1\sigma$ error bars, it is typical to see each violin include both the
injected value and zero. While negative power is non-physical, the nature of
linear fitting with Equation 30 allows for any real-valued ORF amplitude. In
this case, the negative values represent an inverted HD curve, which, if
small, we interpret as noise fluctuations.

Figure 5: The recovered GWB spectrum from the (broadband-normalized) PFOS
using CURN sampling for parameter estimation. The distributions on each
frequency represent the uncertainty sampled noise marginalized distribution on
the PFOS’ recovery.

###  V.2 Per-frequency PSD recovery

We now more systematically test the PSD recovery frequency-by-frequency using
the suite of simulations described in Section Section IV. Each simulation’s
CURN posterior distribution is sampled sufficiently with MCMC techniques
before $10^{3}$ parameter vectors are randomly selected to pass through the
noise-marginalized PFOS. From the distributions of $\hat{S}(f_{k})$ and
$\sigma_{\hat{S}(f_{k})}$, we use uncertainty sampling with $10^{2}$ draws
from each Gaussian to approximate the full distribution of $S(f_{k})$, thereby
incorporating both the Bayesian CURN parameter uncertainty and the PFOS
estimator uncertainty. We then find the percentiles at which the injected GWB
amplitude lies within these full distributions, and use these to construct P-P
plots for each frequency and each tested GWB amplitude. The behavior of these
P-P plots can be characterized with respect to a diagonal, unbiased estimator
[see, e.g., 11, and references therein].

Figure 6: Individual frequency P-P plots at a power-law amplitude of
$A_{\mathrm{gw}}=10^{-14}$ for the (broadband-normalized) PFOS using pair-
covariant, uncertainty sampling, and noise marginalization from a Bayesian
CURN search. The shaded regions indicate the 1, 2, and 3 sigma contours from a
perfectly unbiased estimator for these 100 simulations. Each frequency is
denoted by its multiple of $T_{\rm span}^{-1}$.

Shown in Figure 6 are P-P plots representing the first five frequencies of the
$100$ simulations with the largest GWB amplitude of
$A_{\mathrm{gw}}=10^{-14}$. These strong-signal regime simulations show the
largest bias, and are representative of a worst-case scenario for the PFOS. We
see the first four frequencies show mostly unbiased behaviors with an
estimator mean that only slightly underestimates the true distribution while
remaining well within the $3\sigma$ contours. The fifth frequency briefly
departs this $3\sigma$ contour, indicating a slight underestimation of
$S(f_{k})$ at that frequency. Higher frequencies (not shown) exhibit the same
unbiased behavior as the lowest four frequencies. The slight bias in the fifth
frequency is also present with weaker GWB amplitudes, but less pronounced.
This isolation of a particular problematic frequency hints at a potential
issue with the PTA configuration used in the simulations. This may result from
the irregular timing cadences and gaps leaving poor sensitivity in this
particular frequency, as more ideal simulations with even timing cadence show
no bias at any frequency. More work must be done to determine exactly the
cause of this, but since it does not appear to be an issue with the PFOS
itself this remains beyond the scope of the current work.

We also test the narrowband-normalized PFOS which is more agnostic to the
assumed spectral shape. The narrowband-normalized PFOS is slightly biased in
situations in which there are detectable correlations at other frequencies of
the simulated data. Figure 7 shows a P-P plot with the narrowband-normalized
PFOS for the first five GW frequencies. Each frequency remains broadly
consistent within the $3\sigma$ contours with only a minor deviation in the
fourth frequency. We note however that the mean at every frequency is
marginally overestimated. The uncertainty in the estimators for each frequency
are also inflated due to the mis-modeled spectral shape. For these reasons,
the more accurate and robust solution provided by the broadband-normalized
PFOS are preferred whenever possible.

Figure 7:  Individual frequency P-P plots for a power-law amplitude of
$A_{\mathrm{gw}}=10^{-14}$ for the narrowband-normalized PFOS using pair-
covariance, uncertainty sampling, and noise marginalization from a Bayesian
CURN search. The shaded regions indicate the $1$, $2$, and $3\sigma$ contours
from a perfectly unbiased estimator for these $100$ simulations. Each
frequency is denoted by its multiple of $T_{\mathrm{span}}^{-1}$.

###  V.3 Broadband spectrum recovery

In addition to agnostic spectral recovery, we also test the performance of
fitting the PFOS’ PSD estimates with a parameterized spectral model. We opt
for a power-law so that we can directly compare the PFOS model fitting to the
traditional OS. As with the individual frequency tests, we source our CURN
power-law estimates from Bayesian CURN analyses. Both the PFOS and traditional
OS methods incorporate pulsar-pair covariance and, to provide a fair
comparison, are provided with the same vectors of parameter estimates. Given
the individual PSD estimates resulting from the PFOS, we carry out a least-
squares fit for a power-law amplitude (as usual referenced to a frequency of
$f={\rm yr}^{-1}$) with the spectral index $\gamma$ set to the value of the
Bayesian CURN power-law spectral index parameter to exactly mimic the process
of the traditional OS’ fitting. With noise marginalization and uncertainty
sampling for both the PFOS and traditional OS, we find the percentiles at
which the injected GWB amplitude lies in each of our $100$ simulations.

Figure 8: Broadband P-P plots for a power-law amplitude of
$A_{\mathrm{gw}}=10^{-14}$ for a linear power-law fit on the pair-covariant
(broadband-normalized) PFOS and pair-covariant OS using uncertainty sampling
and noise marginalization from a Bayesian CURN search. The shaded regions
indicate the 1, 2, and 3 sigma contours from a perfectly unbiased estimator
for these 100 simulations.

The results for $A_{\mathrm{gw}}=10^{-14}$ are shown in Figure 8. As with the
individual frequency analysis, the largest GWB amplitude produces the greatest
deviations from an unbiased estimator, such that Figure 8 again represents a
worst-case scenario. We see that the traditional OS with pulsar-pair
covariance performs well and remains within the $3\sigma$ confidence contours
at all percentiles, even in this strong signal regime. By contrast, the power-
law fit on the PFOS analysis products produce minor deviations. There is a
slight sag, especially with higher percentiles, which indicate that the PFOS’
spectral fit and uncertainty are both slightly too large. The fifth frequency
problem shown in Figure 6 may be causing some of these discrepancies, as more
idealized simulations with even cadences show even better agreement between
the two methods.

This test demonstrates that the PFOS’ ability to accurately measure the GWB
PSD is sufficient to then use in a full spectral model fit. Despite the fact
that we remove the power-law assumption of the traditional OS, this method can
produce nearly equivalent results when both methods correctly account for
pulsar-pair covariance.

###  V.4 Stochastic background plus single GW signal

As a final test, we explore the potential for the PFOS to provide early
indicators of a single supermassive black-hole binary system’s GW signal
resounding above the stochastic GWB. For this test, we opt for a pair of PTA
simulations: one with only a GWB signal, and another with exactly the same GWB
signal realization plus an additional single source of continuous
gravitational waves (CWs).

The parameters of this simulation are chosen to be an interesting test, rather
than a realistic possibility. For the GWB, we use an amplitude of
$A_{\mathrm{gw}}=3\times 10^{-15}$ and a spectral index of
$\gamma_{\mathrm{gw}}=13/3$. The CW was injected as a single face-on
supermassive black-hole binary with a chirp mass of $\mathcal{M}=1.1\times
10^{9}M_{\odot}$, a distance of $D_{\mathrm{CW}}=100$ Mpc, and an Earth-term
GW frequency corresponding exactly with the sixth frequency bin in the PTA,
$f_{\rm CW}=1.46\times 10^{-8}\,\mathrm{Hz}$. We also opt to exclude any
binary evolution and pulsar-term contributions to ensure the the effects of
the binary are constrained to this Earth-term frequency.

The PSD of timing residuals induced by a face-on supermassive black-hole
binary—averaged over phase, polarization, and sky position—can be written as
[47]

$S_{\mathrm{CW}}(f_{\mathrm{CW}})=\frac{8\times
10^{-15}}{12\pi^{2}\,f^{2}_{\mathrm{CW}}}\left(\frac{\mathcal{M}}{10^{9}M_{\odot}}\right)^{10/3}\left(\frac{100\,\mathrm{Mpc}}{D_{\mathrm{CW}}}\right)^{2}\left(\frac{f_{\mathrm{CW}}}{50\,\mathrm{nHz}}\right)^{4/3}.$
(37)

With the aforementioned parameters, the combined PSD (i.e., GWB$+$CW) at the
sixth frequency bin will be roughly five times larger than without the CW,
without affecting any other frequencies. However, unlike the GWB, a CW’s
inter-pulsar correlation signature with a finite number of pulsars will only
look approximately like the Hellings & Downs curve when binned by angular
separations, since some pulsars will be more affected than others based on
their proximity to the GW source [48]. As such, the recovered PSD, which we
model as Hellings & Downs correlated, will only be an approximation, and may
not match well in all cases.

Figure 9: The recovered GWB spectrum from the PFOS using CURN sampling for
parameter estimation with and without the CW present. The distributions on
each frequency represent the uncertainty-sampled noise-marginalized
distribution on the (broadband-normalized) PFOS’ recovery. The dashed black
line represents the injected GWB signal at an amplitude of $3\times 10^{-15}$,
and the green star represents the expected PSD value at the 6th frequency due
to the single source.

For this test, we carry out an initial Bayesian power-law CURN search to
provide our PTA parameter estimates, even when the CW was injected. This then
mimics a situation in which the existence of the CW is unknown, thereby
testing the ability for the PFOS to pick up on the excess PSD as a first-look
discovery pipeline for single GW sources. Figure 9 shows the uncertainty-
sampled and noise-marginalized PFOS (with pair covariance) PSD recovery for
identical GWB simulations, with the top panel representing the GWB-only
injection and the bottom panel containing both the GWB and CW injected at the
sixth frequency. We see that even when the CW is unmodeled in the initial
Bayesian PTA parameter estimation stage, the PFOS can pick up on the excess
pulsar-pair correlated PSD caused by a CW. There are other differences to note
however: most frequencies remain mostly unaffected by the CW (like those above
the sixth frequency), some are modified slightly (like the first two
frequencies), and the fifth frequency changes significantly along with the
sixth frequency. We suspect that the reason for this is the now mis-estimated
spectral shape in the pair-wise correlated measurements in Equation 28.

Despite the differences, nearly all PSD recoveries are performing well. Every
frequency in the GWB-only case includes the injected value within the 95%
credible region of the uncertainty sampled distributions. This remains the
case for the GWB$+$CW injection for every frequency except the CW frequency.
This frequency has the injected GWB+CW signal at just below the $3\sigma$
level; this isn’t too surprising, as neglecting the CW’s influence in the GWB
leads to slightly poorer estimates of the total noise and the fact that the
CW’s correlations will not match those of the HD.

While this test was mostly a proof of concept, we plan to further explore the
ability of the PFOS to find the first CW sources in a future study.

##  VI Discussion & Conclusions

In this paper we presented a method to generalize the PTA optimal statistic to
allow for a more general spectral model. We changed the Fourier-domain
covariance matrix of the PTA from a model which assumes a power-law GWB to one
which measures the GWB power spectrum at individual frequencies. The PFOS is
capable of estimating the spectral shape of a correlated GWB in a fraction of
the computational time of that of fully-correlated Bayesian PTA analyses; in
fact the PFOS’ computational time is dominated by the initial Bayesian common
uncorrelated red noise search needed for noise and signal parameter estimates.

In building the PFOS, we had to introduce an additional assumption of
independent PSD estimates at each frequency across the PTA. Despite this
slight mis-modeling, we showed in subsection V.2 that this assumption is valid
for individual frequency measurements, and in subsection V.3 we showed that we
can use these independent measurements to fit broadband spectral models, even
in our realistic data sets. It is expected that with more pulsars, or more
evenly-spaced timing measurements that are encountered in newer data sets
[49], frequency covariance will become even less important.

Additionally, we found that in order to properly normalize the PFOS, it was
important to introduce an estimate of the GWB spectral shape. This manifested
as the unit-amplitude spectral shape at frequency $f$, $\Phi(f_{k})$, that is
analogous to the unit-power-law spectral shape, $\hat{\phi}$, in the
traditional OS. We opted to treat this spectral shape estimate in a similar
manner to the pulsar-pair covariant OS where this estimate is dictated by the
CURN model parameters used to construct the (broadband-normalized) PFOS.
subsection V.2 and subsection V.3 shows that using this strategy results in
consistent measurements with injected values. We also introduced a separate
form of the PFOS called the narrowband-normalized PFOS, which is entirely
agnostic to the spectral shape. This benefit comes at the cost of introducing
slight bias when GWB-induced correlations are detectable in other frequencies.
subsection V.2 showed that despite this purposeful mis-modeling, the
narrowband-normalized PFOS still delivers valid individual-frequency
measurements. However, we recommend that unless entirely independent frequency
measurements are necessary, the more accurate broadband-normalized PFOS should
be used.

We also touched on the choice of the Bayesian CURN search in Section III and
Section IV, where we found that using a variable spectral index power-law CURN
model gave the best recovery, even in cases where there was an unmodeled
single-source in the spectrum, as shown in subsection V.4. We also attempted
to use models that are even more agnostic to the GWB spectral shape, such as
modeling the CURN GWB PSD estimates independently; however, this resulted in
significantly worse performance when these estimates were passed through to
the broadband-normalized PFOS. This was likely due to the PSD estimates of the
highest frequencies spanning multiple orders of magnitude, resulting in
drastic changes to the spectral shape between the posterior draws passed to
the PFOS during the noise marginalization stage.

The PFOS represents the most general implementation of the PTA optimal
statistic, relaxing the assumption of an assumed-shape GWB spectrum. In
conjunction with other recently-developed methods like the multiple component
optimal statistic [35], we now have the ability to characterize the structure
of the spatial correlations and spectrum induced by a GWB at each GW
frequency. There are several directions to continue development of the PFOS
and to employ its flexibility. One may be able to avoid making any assumptions
about the nature of the GWB spectrum (even for the purpose of providing an
initial estimate) by restructuring the PFOS’ $\chi^{2}$-fit to simultaneously
model the contributions from all frequency bins. Initial work on this scheme
shows some promise, though it is being plagued by numerical instability from
the far larger number of elements in the pair covariance matrix (now
$N_{\mathrm{pairs}}^{2}\times N_{\mathrm{freq}}^{2}$), and the correspondingly
larger condition numbers. We plan to continue development of this approach in
hopes of further strengthening the PFOS’ accuracy. We also plan to continue
exploring the effects of single GW sources on the correlated spectrum
estimated by the PFOS. With the inclusion of anisotropic search techniques
there may even be the potential for source localization [44, 10]. A key goal
is to use the PFOS’ capabilities to create a fast, all-frequency power search
for continuous gravitational wave signals. This would act as a first-look
statistic that may motivate, or act as proposal distributions for, subsequent
template-based searches.

##  VII Acknowledgements

The original coding implementation of the PTA optimal statistic’s pair
covariance was developed by Patrick Meyers and Michele Vallisneri, based on
Allen and Romano [33]. We thank Holly Krynicki for many efficiency
improvements in this code, which led the way for many other improvements made
as part of this work. We thank William Lamb, Polina Petrov, and Levi Schult
for cross-checks on both the mathematical and coding aspects of the techniques
developed in this paper. We also thank Shashwat Sardesai for conversations
relating to pulsar pair-covariance implementation. SRT and KAG acknowledge
support from an NSF CAREER #2146016. SRT also acknowledges support from NSF
AST-2007993 and NSF AST-2307719. JDR acknowledges suppport from start-up funds
from the University of Texas Rio Grande Valley. The authors are members of the
NANOGrav collaboration, which receives support from NSF Physics Frontiers
Center award number 1430284 and 2020265\. This work was conducted in part
using the resources of the Advanced Computing Center for Research and
Education (ACCRE) at Vanderbilt University, Nashville, TN.

## References

  * Foster III [1990] R. S. Foster III,  _Constructing a pulsar timing array_ (University of California, Berkeley, 1990).
  * Rosado  _et al._ [2015] P. A. Rosado, A. Sesana, and J. Gair, Expected properties of the first gravitational wave signal detected with pulsar timing arrays, Monthly Notices of the Royal Astronomical Society 451, 2417 (2015).
  * Kelley  _et al._ [2018] L. Z. Kelley, L. Blecha, L. Hernquist, A. Sesana, and S. R. Taylor, Single sources in the low-frequency gravitational wave sky: properties and time to detection by pulsar timing arrays, Monthly Notices of the Royal Astronomical Society 477, 964 (2018).
  * Hellings and Downs [1983] R. W. Hellings and G. S. Downs, Upper limits on the isotropic gravitational radiation background from pulsar timing analysis, The Astrophysical Journal 265, L39 (1983).
  * Agazie  _et al._ [2023a] G. Agazie, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, P. T. Baker, B. Bécsy, L. Blecha, A. Brazier, P. R. Brook, S. Burke-Spolaor, R. Burnette, R. Case, M. Charisi, S. Chatterjee, K. Chatziioannou, B. D. Cheeseboro, S. Chen, T. Cohen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, K. Crowter, C. J. Cutler, M. E. Decesar, D. Degan, P. B. Demorest, H. Deng, T. Dolch, B. Drachler, J. A. Ellis, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, N. Garver-Daniels, P. A. Gentile, K. A. Gersbach, J. Glaser, D. C. Good, K. Gültekin, J. S. Hazboun, S. Hourihane, K. Islo, R. J. Jennings, A. D. Johnson, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, T. C. Klein, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, N. Lewandowska, T. B. Littenberg, T. Liu, A. Lommen, D. R. Lorimer, J. Luo, R. S. Lynch, C.-P. Ma, D. R. Madison, M. A. Mattson, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, P. M. Meyers, C. M. F. Mingarelli, A. Mitridate, P. Natarajan, C. Ng, D. J. Nice, S. K. Ocker, K. D. Olum, T. T. Pennucci, B. B. P. Perera, P. Petrov, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, J. D. Romano, S. C. Sardesai, A. Schmiedekamp, C. Schmiedekamp, K. Schmitz, L. Schult, B. J. Shapiro-Albert, X. Siemens, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, J. P. Sun, A. Susobhanan, J. K. Swiggum, J. Taylor, S. R. Taylor, J. E. Turner, C. Unal, M. Vallisneri, R. van Haasteren, S. J. Vigeland, H. M. Wahl, Q. Wang, C. A. Witt, O. Young, and Nanograv Collaboration, The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background, ApJ 951, L8 (2023a), arXiv:2306.16213 [astro-ph.HE] .
  * Antoniadis  _et al._ [2023] J. Antoniadis, P. Arumugam, S. Arumugam, S. Babak, M. Bagchi, A.-S. B. Nielsen, C. Bassa, A. Bathula, A. Berthereau, M. Bonetti, _et al._ , The second data release from the european pulsar timing array-iii. search for gravitational wave signals, Astronomy & Astrophysics 678, A50 (2023).
  * Reardon  _et al._ [2023] D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane, J. Askew, _et al._ , Search for an isotropic gravitational-wave background with the parkes pulsar timing array, ApJ 951, L6 (2023).
  * Zic  _et al._ [2023] A. Zic, D. J. Reardon, A. Kapur, G. Hobbs, R. Mandow, M. Curyło, R. M. Shannon, J. Askew, M. Bailes, N. D. R. Bhat, A. Cameron, Z.-C. Chen, S. Dai, V. D. Marco, Y. Feng, M. Kerr, A. Kulkarni, M. E. Lower, R. Luo, R. N. Manchester, M. T. Miles, R. S. Nathan, S. Osłowski, A. F. Rogers, C. J. Russell, R. Spiewak, N. Thyagarajan, L. Toomey, S. Wang, L. Zhang, S. Zhang, and X.-J. Zhu, The parkes pulsar timing array third data release, Publications of the Astronomical Society of Australia , 1 (2023).
  * Xu  _et al._ [2023] H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. Nicolas Caballero, J. Yuan, Y. Xu, J. Wang, L. Hao, J. Luo, K. Lee, J. Han, P. Jiang, Z. Shen, M. Wang, N. Wang, R. Xu, X. Wu, R. Manchester, L. Qian, X. Guan, M. Huang, C. Sun, and Y. Zhu, Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I, Research in Astronomy and Astrophysics 23, 075024 (2023), arXiv:2306.16216 [astro-ph.HE] .
  * Pol  _et al._ [2022] N. Pol, S. R. Taylor, and J. D. Romano, Forecasting pulsar timing array sensitivity to anisotropy in the stochastic gravitational wave background, The Astrophysical Journal 940, 173 (2022).
  * Lamb  _et al._ [2023] W. G. Lamb, S. R. Taylor, and R. van Haasteren, Rapid refitting techniques for bayesian spectral characterization of the gravitational wave background using pulsar timing arrays, Phys. Rev. D 108, 103019 (2023).
  * Taylor  _et al._ [2020] S. R. Taylor, R. van Haasteren, and A. Sesana, From bright binaries to bumpy backgrounds: Mapping realistic gravitational wave skies with pulsar-timing arrays, Physical Review D 102, 084039 (2020).
  * Phinney [2001] E. S. Phinney, A practical theorem on gravitational wave backgrounds, Mon. Not. R. Astron. Soc 000, 1 (2001).
  * Agazie  _et al._ [2023b] G. Agazie, A. Anumarlapudi, A. M. Archibald, P. T. Baker, B. Bécsy, L. Blecha, A. Bonilla, A. Brazier, P. R. Brook, S. Burke-Spolaor, R. Burnette, R. Case, J. A. Casey-Clyde, M. Charisi, S. Chatterjee, K. Chatziioannou, B. D. Cheeseboro, S. Chen, T. Cohen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, K. Crowter, C. J. Cutler, D. J. D’Orazio, M. E. Decesar, D. Degan, P. B. Demorest, H. Deng, T. Dolch, B. Drachler, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, E. Gardiner, N. Garver-Daniels, P. A. Gentile, K. A. Gersbach, J. Glaser, D. C. Good, K. Gültekin, J. S. Hazboun, S. Hourihane, K. Islo, R. J. Jennings, A. Johnson, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, N. Lewandowska, T. B. Littenberg, T. Liu, J. Luo, R. S. Lynch, C.-P. Ma, D. R. Madison, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, P. M. Meyers, C. M. F. Mingarelli, A. Mitridate, P. Natarajan, C. Ng, D. J. Nice, S. K. Ocker, K. D. Olum, T. T. Pennucci, B. B. P. Perera, P. Petrov, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, J. D. Romano, J. C. Runnoe, S. C. Sardesai, A. Schmiedekamp, C. Schmiedekamp, K. Schmitz, L. Schult, B. J. Shapiro-Albert, X. Siemens, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, J. P. Sun, A. Susobhanan, J. K. Swiggum, J. Taylor, S. R. Taylor, J. E. Turner, C. Unal, M. Vallisneri, S. J. Vigeland, J. M. Wachter, H. M. Wahl, Q. Wang, C. A. Witt, D. Wright, O. Young, and Nanograv Collaboration, The NANOGrav 15 yr Data Set: Constraints on Supermassive Black Hole Binaries from the Gravitational-wave Background, ApJ 952, L37 (2023b), arXiv:2306.16220 [astro-ph.HE] .
  * Wu  _et al._ [2024] Y.-M. Wu, Z.-C. Chen, and Q.-G. Huang, Cosmological interpretation for the stochastic signal in pulsar timing arrays, Science China Physics, Mechanics & Astronomy 67, 240412 (2024).
  * Bécsy  _et al._ [2022] B. Bécsy, N. J. Cornish, and L. Z. Kelley, Exploring realistic nanohertz gravitational-wave backgrounds, The Astrophysical Journal 941, 119 (2022).
  * Agazie  _et al._ [2024] G. Agazie, P. T. Baker, B. Bécsy, L. Blecha, A. Brazier, P. R. Brook, L. Brown, S. Burke-Spolaor, J. A. Casey-Clyde, M. Charisi, _et al._ , The nanograv 15 yr data set: Looking for signs of discreteness in the gravitational-wave background, arXiv preprint arXiv:2404.07020 (2024).
  * Sampson  _et al._ [2015] L. Sampson, N. J. Cornish, and S. T. McWilliams, Constraining the solution to the last parsec problem with pulsar timing, Physical Review D 91, 084055 (2015).
  * Taylor  _et al._ [2017] S. R. Taylor, J. Simon, and L. Sampson, Constraints on the dynamical environments of supermassive black-hole binaries using pulsar-timing arrays, Physical Review Letters 118, 181102 (2017).
  * Afzal  _et al._ [2023] A. Afzal, G. Agazie, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, P. T. Baker, B. Bécsy, J. J. Blanco-Pillado, L. Blecha, K. K. Boddy, A. Brazier, P. R. Brook, S. Burke-Spolaor, R. Burnette, R. Case, M. Charisi, S. Chatterjee, K. Chatziioannou, B. D. Cheeseboro, S. Chen, T. Cohen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, K. Crowter, C. J. Cutler, M. E. Decesar, D. Degan, P. B. Demorest, H. Deng, T. Dolch, B. Drachler, R. von Eckardstein, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, N. Garver-Daniels, P. A. Gentile, K. A. Gersbach, J. Glaser, D. C. Good, L. Guertin, K. Gültekin, J. S. Hazboun, S. Hourihane, K. Islo, R. J. Jennings, A. D. Johnson, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, V. S. H. Lee, N. Lewandowska, R. R. Lino Dos Santos, T. B. Littenberg, T. Liu, D. R. Lorimer, J. Luo, R. S. Lynch, C.-P. Ma, D. R. Madison, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, P. M. Meyers, C. M. F. Mingarelli, A. Mitridate, J. Nay, P. Natarajan, C. Ng, D. J. Nice, S. K. Ocker, K. D. Olum, T. T. Pennucci, B. B. P. Perera, P. Petrov, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, J. D. Romano, S. C. Sardesai, A. Schmiedekamp, C. Schmiedekamp, K. Schmitz, T. Schröder, L. Schult, B. J. Shapiro-Albert, X. Siemens, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, P. Stratmann, J. P. Sun, A. Susobhanan, J. K. Swiggum, J. Taylor, S. R. Taylor, T. Trickle, J. E. Turner, C. Unal, M. Vallisneri, S. Verma, S. J. Vigeland, H. M. Wahl, Q. Wang, C. A. Witt, D. Wright, O. Young, K. M. Zurek, and Nanograv Collaboration, The NANOGrav 15 yr Data Set: Search for Signals from New Physics, ApJ 951, L11 (2023), arXiv:2306.16219 [astro-ph.HE] .
  * Kibble [1976] T. W. Kibble, Topology of cosmic domains and strings, Journal of Physics A: Mathematical and General 9, 1387 (1976).
  * Johnson  _et al._ [2023] A. D. Johnson, P. M. Meyers, P. T. Baker, N. J. Cornish, J. S. Hazboun, T. B. Littenberg, J. D. Romano, S. R. Taylor, M. Vallisneri, S. J. Vigeland, K. D. Olum, X. Siemens, J. A. Ellis, R. van Haasteren, S. Hourihane, G. Agazie, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, L. Blecha, A. Brazier, P. R. Brook, S. Burke-Spolaor, B. Bécsy, J. A. Casey-Clyde, M. Charisi, S. Chatterjee, K. Chatziioannou, T. Cohen, J. M. Cordes, F. Crawford, H. T. Cromartie, K. Crowter, M. E. DeCesar, P. B. Demorest, T. Dolch, B. Drachler, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, N. Garver-Daniels, P. A. Gentile, J. Glaser, D. C. Good, K. Gültekin, R. J. Jennings, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, N. Lewandowska, T. Liu, D. R. Lorimer, R. S. Lynch, C.-P. Ma, D. R. Madison, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, C. M. F. Mingarelli, A. Mitridate, C. Ng, D. J. Nice, S. K. Ocker, T. T. Pennucci, B. B. P. Perera, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, S. C. Sardesai, C. Schmiedekamp, A. Schmiedekamp, K. Schmitz, B. J. Shapiro-Albert, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, A. Susobhanan, J. K. Swiggum, J. E. Turner, C. Unal, H. M. Wahl, C. A. Witt, and O. Young, The NANOGrav 15-year Gravitational-Wave Background Analysis Pipeline, arXiv e-prints , arXiv:2306.16223 (2023), arXiv:2306.16223 [astro-ph.HE] .
  * Taylor [2021] S. R. Taylor,  _Nanohertz gravitational wave astronomy_ (CRC Press, 2021).
  * Taylor  _et al._ [2022] S. R. Taylor, J. Simon, L. Schult, N. Pol, and W. G. Lamb, A parallelized bayesian approach to accelerated gravitational-wave background characterization, Physical Review D 105, 084049 (2022).
  * Freedman  _et al._ [2023] G. E. Freedman, A. D. Johnson, R. van Haasteren, and S. J. Vigeland, Efficient gravitational wave searches with pulsar timing arrays using hamiltonian monte carlo, Physical Review D 107, 043013 (2023).
  * Hourihane  _et al._ [2023] S. Hourihane, P. Meyers, A. Johnson, K. Chatziioannou, and M. Vallisneri, Accurate characterization of the stochastic gravitational-wave background with pulsar timing arrays by likelihood reweighting, Physical Review D 107, 084045 (2023).
  * Anholm  _et al._ [2009] M. Anholm, S. Ballmer, J. D. Creighton, L. R. Price, and X. Siemens, Optimal strategies for gravitational wave stochastic background searches in pulsar timing data, Physical Review D - Particles, Fields, Gravitation and Cosmology 79, 10.1103/PhysRevD.79.084030 (2009).
  * Allen and Romano [1999] B. Allen and J. D. Romano, Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities, Phys. Rev. D 59, 102001 (1999), arXiv:gr-qc/9710117 [gr-qc] .
  * Chamberlin  _et al._ [2015] S. J. Chamberlin, J. D. Creighton, X. Siemens, P. Demorest, J. Ellis, L. R. Price, and J. D. Romano, Time-domain implementation of the optimal cross-correlation statistic for stochastic gravitational-wave background searches in pulsar timing data, Physical Review D 91, 044048 (2015).
  * Taylor  _et al._ [2017] S. R. Taylor, L. Lentati, S. Babak, P. Brem, J. R. Gair, A. Sesana, and A. Vecchio, All correlations must die: Assessing the significance of a stochastic gravitational-wave background in pulsar timing arrays, Phys. Rev. D 95, 042002 (2017), arXiv:1606.09180 [astro-ph.IM] .
  * Cornish and Sampson [2016a] N. J. Cornish and L. Sampson, Towards robust gravitational wave detection with pulsar timing arrays, Physical Review D 93, 104047 (2016a).
  * Romano  _et al._ [2021] J. D. Romano, J. S. Hazboun, X. Siemens, and A. M. Archibald, Common-spectrum process versus cross-correlation for gravitational-wave searches using pulsar timing arrays, Physical Review D 103, 063027 (2021).
  * Allen and Romano [2023] B. Allen and J. D. Romano, Hellings and downs correlation of an arbitrary set of pulsars, Physical Review D 108, 043026 (2023).
  * Allen [2023] B. Allen, Variance of the hellings-downs correlation, Physical Review D 107, 043018 (2023).
  * Sardesai  _et al._ [2023] S. C. Sardesai, S. J. Vigeland, K. A. Gersbach, and S. R. Taylor, Generalized optimal statistic for characterizing multiple correlated signals in pulsar timing arrays, Physical Review D 108, 124081 (2023).
  * Vallisneri  _et al._ [2023] M. Vallisneri, P. M. Meyers, K. Chatziioannou, and A. J. K. Chua, Posterior predictive checking for gravitational-wave detection with pulsar timing arrays. i. the optimal statistic, Phys. Rev. D 108, 123007 (2023).
  * Meyers  _et al._ [2023] P. M. Meyers, K. Chatziioannou, M. Vallisneri, and A. J. K. Chua, Posterior predictive checking for gravitational-wave detection with pulsar timing arrays. ii. posterior predictive distributions and pseudo-bayes factors, Phys. Rev. D 108, 123008 (2023).
  * Vigeland  _et al._ [2018] S. J. Vigeland, K. Islo, S. R. Taylor, and J. A. Ellis, Noise-marginalized optimal statistic: A robust hybrid frequentist-bayesian statistic for the stochastic gravitational-wave background in pulsar timing arrays, Physical Review D 98, 044003 (2018).
  * Siemens  _et al._ [2013] X. Siemens, J. Ellis, F. Jenet, and J. D. Romano, The stochastic background: scaling laws and time to detection for pulsar timing arrays, Classical and Quantum Gravity 30, 224015 (2013).
  * Arzoumanian  _et al._ [2020] Z. Arzoumanian, P. T. Baker, H. Blumer, B. Bécsy, A. Brazier, P. R. Brook, S. Burke-Spolaor, S. Chatterjee, S. Chen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, M. E. DeCesar, P. B. Demorest, T. Dolch, J. A. Ellis, E. C. Ferrara, W. Fiore, E. Fonseca, N. Garver-Daniels, P. A. Gentile, D. C. Good, J. S. Hazboun, A. M. Holgado, K. Islo, R. J. Jennings, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, J. S. Key, N. Laal, M. T. Lam, T. W. J. Lazio, D. R. Lorimer, J. Luo, R. S. Lynch, D. R. Madison, M. A. McLaughlin, C. M. F, D. R. Stinebring, K. Stovall, J. P. Sun, J. K. Swiggum, S. R. Taylor, J. E. Turner, M. Vallisneri, S. J. Vigeland, and C. A. Witt, The NANOGrav 12.5-yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background, ApJ 33, 34 (2020).
  * Hazboun  _et al._ [2023] J. S. Hazboun, P. M. Meyers, J. D. Romano, X. Siemens, and A. M. Archibald, Analytic distribution of the optimal cross-correlation statistic for stochastic gravitational-wave-background searches using pulsar timing arrays, Phys. Rev. D 108, 104050 (2023).
  * Cornish and Sampson [2016b] N. J. Cornish and L. Sampson, Towards robust gravitational wave detection with pulsar timing arrays, Phys. Rev. D 93, 104047 (2016b).
  * Ballmer [2006] S. W. Ballmer, A radiometer for stochastic gravitational waves, Classical and Quantum Gravity 23, S179 (2006), arXiv:gr-qc/0510096 [gr-qc] .
  * Agazie  _et al._ [2023c] G. Agazie, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, P. T. Baker, B. Bécsy, L. Blecha, A. Brazier, P. R. Brook, S. Burke-Spolaor, J. A. Casey-Clyde, M. Charisi, S. Chatterjee, T. Cohen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, K. Crowter, M. E. DeCesar, P. B. Demorest, T. Dolch, B. Drachler, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, E. Gardiner, N. Garver-Daniels, P. A. Gentile, J. Glaser, D. C. Good, K. Gültekin, J. S. Hazboun, R. J. Jennings, A. D. Johnson, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, N. Lewandowska, T. Liu, D. R. Lorimer, J. Luo, R. S. Lynch, C.-P. Ma, D. R. Madison, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, C. M. F. Mingarelli, A. Mitridate, C. Ng, D. J. Nice, S. K. Ocker, K. D. Olum, T. T. Pennucci, B. B. P. Perera, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, J. D. Romano, S. C. Sardesai, A. Schmiedekamp, C. Schmiedekamp, K. Schmitz, L. Schult, B. J. Shapiro-Albert, X. Siemens, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, A. Susobhanan, J. K. Swiggum, S. R. Taylor, J. E. Turner, C. Unal, M. Vallisneri, S. J. Vigeland, H. M. Wahl, C. A. Witt, and O. Young, The NANOGrav 15 yr Data Set: Search for Anisotropy in the Gravitational-wave Background, ApJ 956, L3 (2023c), arXiv:2306.16221 [astro-ph.HE] .
  * Vallisneri [2020] M. Vallisneri, libstempo: Python wrapper for tempo2, Astrophysics Source Code Library, record ascl:2002.017 , ascl:2002.017 (2020).
  * Pol  _et al._ [2021] N. S. Pol, S. R. Taylor, L. Z. Kelley, S. J. Vigeland, J. Simon, S. Chen, Z. Arzoumanian, P. T. Baker, B. Bécsy, A. Brazier, P. R. Brook, S. Burke-Spolaor, S. Chatterjee, J. M. Cordes, N. J. Cornish, F. Crawford, H. Thankful Cromartie, M. E. Decesar, P. B. Demorest, T. Dolch, E. C. Ferrara, W. Fiore, E. Fonseca, N. Garver-Daniels, D. C. Good, J. S. Hazboun, R. J. Jennings, M. L. Jones, A. R. Kaiser, D. L. Kaplan, J. Shapiro Key, M. T. Lam, T. J. W. Lazio, J. Luo, R. S. Lynch, D. R. Madison, A. McEwen, M. A. McLaughlin, C. M. F. Mingarelli, C. Ng, D. J. Nice, T. T. Pennucci, S. M. Ransom, P. S. Ray, B. J. Shapiro-Albert, X. Siemens, I. H. Stairs, D. R. Stinebring, J. K. Swiggum, M. Vallisneri, H. Wahl, C. A. Witt, and Nanograv Collaboration, Astrophysics Milestones for Pulsar Timing Array Gravitational-wave Detection, ApJ 911, L34 (2021), arXiv:2010.11950 [astro-ph.HE] .
  * Ellis  _et al._ [2012] J. Ellis, F. Jenet, and M. McLaughlin, Practical methods for continuous gravitational wave detection using pulsar timing data, The Astrophysical Journal 753, 96 (2012).
  * Cornish and Sesana [2013] N. J. Cornish and A. Sesana, Pulsar timing array analysis for black hole backgrounds, Classical and Quantum Gravity 30, 224005 (2013).
  * Agazie  _et al._ [2023d] G. Agazie, M. F. Alam, A. Anumarlapudi, A. M. Archibald, Z. Arzoumanian, P. T. Baker, L. Blecha, V. Bonidie, A. Brazier, P. R. Brook, S. Burke-Spolaor, B. Bécsy, C. Chapman, M. Charisi, S. Chatterjee, T. Cohen, J. M. Cordes, N. J. Cornish, F. Crawford, H. T. Cromartie, K. Crowter, M. E. Decesar, P. B. Demorest, T. Dolch, B. Drachler, E. C. Ferrara, W. Fiore, E. Fonseca, G. E. Freedman, N. Garver-Daniels, P. A. Gentile, J. Glaser, D. C. Good, K. Gültekin, J. S. Hazboun, R. J. Jennings, C. Jessup, A. D. Johnson, M. L. Jones, A. R. Kaiser, D. L. Kaplan, L. Z. Kelley, M. Kerr, J. S. Key, A. Kuske, N. Laal, M. T. Lam, W. G. Lamb, T. J. W. Lazio, N. Lewandowska, Y. Lin, T. Liu, D. R. Lorimer, J. Luo, R. S. Lynch, C.-P. Ma, D. R. Madison, K. Maraccini, A. McEwen, J. W. McKee, M. A. McLaughlin, N. McMann, B. W. Meyers, C. M. F. Mingarelli, A. Mitridate, C. Ng, D. J. Nice, S. K. Ocker, K. D. Olum, E. Panciu, T. T. Pennucci, B. B. P. Perera, N. S. Pol, H. A. Radovan, S. M. Ransom, P. S. Ray, J. D. Romano, L. Salo, S. C. Sardesai, C. Schmiedekamp, A. Schmiedekamp, K. Schmitz, B. J. Shapiro-Albert, X. Siemens, J. Simon, M. S. Siwek, I. H. Stairs, D. R. Stinebring, K. Stovall, A. Susobhanan, J. K. Swiggum, S. R. Taylor, J. E. Turner, C. Unal, M. Vallisneri, S. J. Vigeland, H. M. Wahl, Q. Wang, C. A. Witt, O. Young, and Nanograv Collaboration, The NANOGrav 15 yr Data Set: Observations and Timing of 68 Millisecond Pulsars, ApJ 951, L9 (2023d), arXiv:2306.16217 [astro-ph.HE] .
  * Isserlis [1918] L. Isserlis, On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables, Biometrika 12, 134 (1918).

##  Appendix A Derivation of the Optimal Statistic’s normalization

In this appendix there are two derivations for the normalization. The first is
for the traditional optimal statistic which assumes a power-law spectral
template, and the second is for the per-frequency optimal statistic which is
agnostic to the spectral template.

###  A1 Traditional Optimal Statistic

In order to calculate the normalization for Equation 8, we set the expected
value of our estimator to the desired quantity,

$\left<\mathcal{N}_{ab}\,\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\bm{\delta
t}_{b}\right>=\Gamma_{ab}A^{2}_{\mathrm{gw}}.$ (38)

To solve this equation, we rearrange the quantities and use two identities.
First, the trace of a scalar is that same scalar, meaning that we can apply
the trace to our estimator. Second, the expectation and trace are both linear
operators that we can freely exchange. After rearranging and adding in the
trace,

$\mathcal{N}_{ab}=\frac{\Gamma_{ab}A^{2}_{\mathrm{gw}}}{\left<\mathrm{tr}\left[\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\bm{\delta
t}_{b}\right]\right>}.$ (39)

From here, we use the cyclic property of the trace (i.e. with matrices
$\bm{A}$, $\bm{B}$, and $\bm{C}$ we can cycle the elements,
$\mathrm{tr}\left[\bm{A}\bm{B}\bm{C}\right]=\mathrm{tr}\left[\bm{B}\bm{C}\bm{A}\right]$).
With this in mind we can cycle the trace and then exchange our trace and
expectation value so that we apply the expectation to only the random
variables in the equation, $\bm{\delta t}_{b}$ and $\bm{\delta t}_{a}^{T}$,

$\mathcal{N}_{ab}=\frac{\Gamma_{ab}A^{2}_{\mathrm{gw}}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\left<\bm{\delta
t}_{b}\,\bm{\delta t}_{a}^{T}\right>\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]}.$
(40)

This allows us to substitute the expected covariance of TOA residuals,
$\left<\bm{\delta t}_{b}\;\bm{\delta
t}_{a}^{T}\right>=\bm{S}_{ba}=\Gamma_{ab}A_{\mathrm{gw}}^{2}\;\hat{\bm{S}}_{ba}$,

$\mathcal{N}_{ab}=\frac{\Gamma_{ab}A^{2}_{\mathrm{gw}}}{\Gamma_{ab}A_{\mathrm{gw}}^{2}\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]}.$
(41)

Canceling the $\Gamma_{ab}A^{2}_{\mathrm{gw}}$ we find our final normalization

$\mathcal{N}_{ab}=\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{-1},$
(42)

implying that our full cross-correlation estimator is then

$\rho_{ab}=\frac{\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\bm{\delta
t}_{b}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]}.$
(43)

###  A2 Per-Frequency Optimal Statistic

The normalization calculation for the PFOS is similar to the traditional OS.
However, we instead start from the rank-reduced form of the pair-wise
estimator in Equation 27. We set our expectation value so that

$\left<\mathcal{N}_{ab}(f_{k})\,\bm{X}_{a}^{T}\tilde{\phi}(f_{k})\bm{X}_{b}\right>=\Gamma_{ab}\,S(f_{k}).$
(44)

Substituting all $\bm{X}$ terms with
$\bm{X}_{a}=\bm{F}^{T}_{a}\bm{P}_{a}^{-1}\bm{\delta t}_{a}$ then solving for
$\mathcal{N}_{ab}(f_{k})$ leads us to

$\mathcal{N}_{ab}(f_{k})=\frac{\Gamma_{ab}\,S(f_{k})}{\left<\bm{\delta
t}^{T}_{a}\bm{P}_{a}^{-1}\bm{F}_{a}\tilde{\phi}(f_{k})\bm{F}^{T}_{b}\bm{P}_{b}^{-1}\bm{\delta
t}_{b}\right>}.$ (45)

Adding in the trace, we cycle it then exchange the trace with the expectation
value to obtain

$\mathcal{N}_{ab}(f_{k})=\frac{\Gamma_{ab}S(f_{k})}{\mathrm{tr}\left[\bm{P}_{a}^{-1}\bm{F}_{a}\tilde{\phi}(f_{k})\bm{F}^{T}_{b}\bm{P}_{b}^{-1}\left<\bm{\delta
t}_{b}\;\bm{\delta t}^{T}_{a}\right>\right]}.$ (46)

Unlike with the traditional OS, we do not assume any particular spectral
shape. Allowing for an arbitrary spectral shape, we write our expected cross-
correlation as a spectrum, $\left<\bm{\delta t}_{b}\;\bm{\delta
t}^{T}_{a}\right>=\Gamma_{ab}\,\bm{F}_{b}\,\phi\,\bm{F}^{T}_{a}$. Substituting
this into our equation then cycling the trace once again leads to

$\mathcal{N}_{ab}(f_{k})=\frac{S(f_{k})}{\mathrm{tr}\left[\bm{F}^{T}_{a}\bm{P}_{a}^{-1}\bm{F}_{a}\tilde{\phi}(f_{k})\bm{F}^{T}_{b}\bm{P}_{b}^{-1}\bm{F}_{b}\phi\right]}.$
(47)

Identifying our matrix substitution for
$\bm{Z}_{a}=\bm{F}^{T}_{a}\bm{P}_{a}^{-1}\bm{F}_{a}$, we can simplify this
equation as

$\mathcal{N}_{ab}(f_{k})=\frac{S(f_{k})}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\phi\right]}.$
(48)

As a final step, we “normalize” our spectrum, meaning that we define a new
quantity

$\Phi(f_{k})=\frac{\phi}{S(f_{k})}.$ (49)

Recall that $\phi$ is a diagonal matrix with elements corresponding to the
Fourier modes of the GWB PSD,
$\phi=\mathrm{diag}(S(f_{1}),S(f_{1}),S(f_{2}),S(f_{2}),...,S(f_{k}),S(f_{k}),...)$.
Also recall that the shape of $\phi$ is ($2k_{\rm max}\times 2k_{\rm max}$)
due to including both the sine and cosine of each frequency. This new
quantity, $\Phi(f_{k})$, represents the spectral shape of the GWB PSD with
unit-value at frequency $f_{k}$. With this substitution we no longer need an
estimate of $S(f_{k})$; rather we need a spectral shape estimate.

Substituting this normalization back into Equation 27, we are left with our
final estimator,

$\rho_{ab}(f_{k})=\frac{\bm{X}^{T}_{a}\tilde{\phi}(f_{k})\bm{X}_{b}}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\Phi(f_{k})\right]}.$
(50)

##  Appendix B Uncertainty on the OS cross-correlation estimators

This appendix contains the derivations for the uncertainty in both the
traditional optimal statistic and the per-frequency optimal statistic through
the variance of the cross-correlation estimators. Both derivations share many
of the same steps for the calculation, and are both calculated under the null
hypothesis of an absence of pulsar-pair correlations. In these cases, the
expectation values of the pair-wise correlated estimators are zero, as was
assumed in the original formulation of the OS [27]. We will denote expectation
values calculated assuming the null hypothesis by using a subscript $0$, so
$\langle\ \rangle_{0}$.

###  B1 Traditional Optimal Statistic

For the traditional optimal statistic, we have the cross-correlation estimator

$\rho_{ab}=\frac{\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\bm{\delta
t}_{b}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]}.$
(51)

Calculating the variance of the OS cross-correlation estimator under the null
hypothesis begins with recalling that the expectation of the individual pair-
wise estimator is zero, so that

$\sigma_{ab,0}^{2}=\left<\rho_{ab}^{2}\right>_{0}-\left<\rho_{ab}\right>^{2}_{0}=\left<\rho_{ab}^{2}\right>_{0}.$
(52)

Substituting for $\rho_{ab}$ in our equation we find

$\sigma_{ab,0}^{2}=\frac{\left<\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\bm{\delta
t}_{b}\,\bm{\delta
t}_{a}^{T}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\bm{\delta
t}_{b}\right>_{0}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}}.$
(53)

Since we have four Gaussian random variables inside our expectation value, we
must isolate these values from the constant matrix products. The first step
involves using a temporary matrix product,
$\bm{Q}_{ab}=\bm{P}_{a}^{-1}\>\hat{\bm{S}}_{ab}\>\bm{P}_{b}^{-1}$. We identify
two instances of this quantity to create

$\sigma_{ab,0}^{2}=\frac{\left<\bm{\delta t}^{T}_{a}\bm{Q}_{ab}\bm{\delta
t}_{b}\,\bm{\delta t}^{T}_{a}\bm{Q}_{ab}\bm{\delta
t}_{b}\right>_{0}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}}$
(54)

From here we move to index notation so that we can isolate our random
variables. The indices $i$, $j$, $k$, and $l$ will represent the individual
TOA residual measurements. This change results in

$\sigma_{ab,0}^{2}=\frac{\sum_{ij,kl}\bm{Q}_{ab,ij}\,\bm{Q}_{ab,kl}\left<\bm{\delta
t}_{a,i}\,\bm{\delta t}_{b,j}\,\bm{\delta t}_{a,k}\,\bm{\delta
t}_{b,l}\right>_{0}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}}$
(55)

We can evaluate the 4th order-moment of zero-mean Gaussian random variables
using Isserlis’ theorem [50]:

$\left<ABCD\right>=\left<AB\right>\left<CD\right>+\left<AC\right>\left<BD\right>+\left<AD\right>\left<BC\right>.$
(56)

With this we write

$\sigma_{ab,0}^{2}=\frac{\sum_{ij,kl}\bm{Q}_{ab,ij}\,\bm{Q}_{ab,kl}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}}\left(\left<\bm{\delta
t}_{a,i}\bm{\delta t}_{b,j}\right>_{0}\left<\bm{\delta t}_{a,k}\bm{\delta
t}_{b,l}\right>_{0}+\left<\bm{\delta t}_{a,i}\bm{\delta
t}_{a,k}\right>_{0}\left<\bm{\delta t}_{b,j}\bm{\delta
t}_{b,l}\right>_{0}+\left<\bm{\delta t}_{a,i}\bm{\delta
t}_{b,l}\right>_{0}\left<\bm{\delta t}_{b,j}\bm{\delta
t}_{a,k}\right>_{0}\right).$ (57)

Next, we identify both the auto-correlation and cross-correlation matrices
inside the expectation values. Since we are assuming the null hypothesis, we
set to zero any cross-correlation quantities and identify auto-correlation
matrices as $\bm{P}_{a,ij}=\left<\bm{\delta t}_{a,i}\bm{\delta
t}_{a,j}\right>_{0}$. This leaves us with

$\sigma_{ab,0}^{2}=\frac{\sum_{ij,kl}\bm{Q}_{ab,ij}\,\bm{Q}_{ab,kl}\,\bm{P}_{a,ik}\,\bm{P}_{b,jl}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}}.$
(58)

Next, we use the fact that $\bm{P}_{a}$ is symmetric to write
$\bm{P}_{a,ik}=\bm{P}_{a,ki}$. When combined with an additional matrix
identity, $\bm{Q}_{ab,kl}=\bm{Q}_{ba,lk}$, we can rearrange the elements of
the equation to be

$\sigma_{ab,0}^{2}=\frac{\sum_{ij,kl}\bm{P}_{a,ki}\,\bm{Q}_{ab,ij}\,\bm{P}_{b,jl}\,\bm{Q}_{ba,lk}}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}}.$
(59)

From here, we move back into matrix notation by identifying the trace of the
matrix products in the numerator, letting us rewrite it as

$\sigma_{ab,0}^{2}=\frac{\mathrm{tr}\left[\bm{P}_{a}\,\bm{Q}_{ab}\,\bm{P}_{b}\,\bm{Q}_{ba}\right]}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}}.$
(60)

Expressing our $\bm{Q}_{ab}$ quantities in terms of $\bm{P}^{-1}_{a}$ and
$\hat{\bm{S}}_{ab}$, and canceling matrix products with their inverse, we are
left with

$\sigma_{ab,0}^{2}=\frac{\mathrm{tr}\left[\hat{\bm{S}}_{ab}\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\right]}{\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{2}},$
(61)

which leads to the final (standard deviation) uncertainty for the cross-
correlation estimators

$\sigma_{ab,0}=\mathrm{tr}\left[\bm{P}_{b}^{-1}\hat{\bm{S}}_{ba}\bm{P}_{a}^{-1}\hat{\bm{S}}_{ab}\right]^{-1/2}.$
(62)

###  B2 Per-Frequency Optimal Statistic

In this subsection, we calculate the uncertainty in the PFOS cross-correlated
PSD estimator

$\rho_{ab}(f_{k})=\frac{\bm{X}^{T}_{a}\tilde{\phi}(f_{k})\bm{X}_{b}}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\Phi(f_{k})\right]}.$
(63)

Much of the calculation is identical to Appendix B1, and as such many
identical steps will be skipped. As was the case with the traditional OS, we
calculate the uncertainty assuming the absence of a correlated-signal, such
that the variance is then given by

$\sigma_{ab,0}(f_{k})^{2}=\left<\rho_{ab,f}^{2}\right>_{0}-\left<\rho_{ab,f}\right>_{0}^{2}=\left<\rho_{ab,f}^{2}\right>_{0}.$
(64)

Substituting for $\rho_{ab}(f_{k})$ leads to

$\sigma_{ab,0}(f_{k})^{2}=\left<\frac{\bm{X}^{T}_{a}\tilde{\phi}(f_{k})\bm{X}_{b}\;\bm{X}^{T}_{a}\tilde{\phi}(f_{k})\bm{X}_{b}}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}_{f}\bm{Z}_{b}\Phi_{f}\right]^{2}}\right>_{0},$
(65)

We then substitute for $\bm{X}_{a}$ and $\bm{X}_{b}$ using
$\bm{X}_{a}=\bm{F}^{T}_{a}\bm{P}_{a}^{-1}\bm{\delta t}_{a}$, and define a
temporary matrix quantity
$\bm{Q}^{\prime}_{ab}(f_{k})=\bm{P}_{a}^{-1}\bm{F}_{a}\tilde{\phi}(f_{k})\bm{F}_{b}^{T}\bm{P}_{b}^{-1}$.
With this, we find

$\sigma_{ab,0}(f_{k})^{2}=\frac{\left<\bm{\delta
t}^{T}_{a}\bm{Q}^{\prime}_{ab}(f_{k})\bm{\delta t}_{b}\,\bm{\delta
t}^{T}_{a}\bm{Q}^{\prime}_{ab}(f_{k})\bm{\delta
t}_{b}\right>_{0}}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\Phi(f_{k})\right]^{2}}.$
(66)

This equation has identical from to Equation 54 from the traditional OS
uncertainty calculation. If we follow through the same steps, we end up in a
similar form as Equation 60 except we now use rank-reduced quantities to
obtain

$\sigma_{ab,0}(f_{k})^{2}=\frac{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\tilde{\phi}(f_{k})\right]}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\Phi(f_{k})\right]^{2}}.$
(67)

Finally, we take the square root to find the standard deviation uncertainty
for the cross-correlated PSD at frequency $f_{k}$,

$\sigma_{ab,0}(f_{k})=\left(\frac{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\tilde{\phi}(f_{k})\right]}{\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\Phi(f_{k})\right]^{2}}\right)^{1/2}.$
(68)

Unlike with the traditional OS, we cannot simplify this equation any further
unless we make an additional assumption. One logical assumption, detailed in
Section III, is that of the narrowband-normalized PFOS. This form of the PFOS
assumes that all GW power comes only from a single frequency. This leads us to
replace $\Phi(f_{k})$ with $\tilde{\phi}(f_{k})$, effectively zeroing the
assumed spectral power in other frequencies. This lets us express the
narrowband-normalized PFOS’ uncertainty as

$\sigma_{ab,0,\mathrm{narrow}}(f_{k})=\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\tilde{\phi}(f_{k})\right]^{-1/2}.$
(69)

##  Appendix C Covariance of pairwise cross-correlation estimates

As PTAs move towards larger signal regimes, the null-hypothesis assumption of
zero covariance between pulsar pair estimators becomes a poor approximation,
and the optimal statistic as a parameter estimator becomes severely biased, as
was shown in subsection II.3. Including the pulsar-pair covariance fixes this
problem. However, the number of elements in this covariance matrix scales as
$N_{\mathrm{pulsars}}^{4}$. As such, this section will focus on deriving more
efficient rank-reduced forms for calculating the elements of this covariance
matrix for both the traditional OS and the PFOS.

###  C1 Traditional Optimal Statistic

With the traditional OS, we assume a power-law form for the GWB PSD. This
assumption allows us to factor out a common GWB amplitude at a reference
frequency of $f_{\mathrm{yr}}$ from our PSD matrix in Equation 6,
$\phi=A^{2}_{\mathrm{gw}}\,\hat{\phi}$, where $\hat{\phi}$ is the unit-
amplitude power-law spectral template.

The elements of the covariance matrix for the traditional OS can be
represented as

$\bm{C}_{ab,cd}=\left<\rho_{ab}\,\rho_{cd}\right>-\left<\rho_{ab}\right>\left<\rho_{cd}\right>.$
(70)

For this calculation, we express our $\rho_{ab}$ quantities in a form that
better translates to the PFOS case later. We write these in a similar manner
to Equation 10 except we opt to factor the normalization and substitute
$\hat{\bm{S}}_{ab}=\bm{F}_{a}\hat{\phi}\bm{F}^{T}_{b}$. With these changes,
the estimator can now be written as

$\rho_{ab}=\mathcal{N}_{ab}\,\bm{\delta
t}^{T}_{a}\,\bm{P}^{-1}_{a}\,\bm{F}_{a}\,\hat{\phi}\,\bm{F}^{T}_{b}\,\bm{P}^{-1}_{b}\,\bm{\delta
t}_{a},$ (71)

where
$\mathcal{N}_{ab}=\mathrm{tr}\left[\bm{Z}_{a}\hat{\phi}\bm{Z}_{b}\hat{\phi}\right]^{-1}$
is the normalization for pulsar pair $ab$, computed in Appendix A1. For
convenience, we now define a temporary matrix quantity

$\bm{Q}_{ab}=\bm{P}^{-1}_{a}\,\bm{F}_{a}\,\hat{\phi}\,\bm{F}^{T}_{b}\,\bm{P}^{-1}_{b},$
(72)

which allows us to rewrite our estimator as

$\rho_{ab}=\mathcal{N}_{ab}\,\bm{\delta t}^{T}_{a}\,\bm{Q}_{ab}\,\bm{\delta
t}_{a}.$ (73)

We now substitute this into our covariance matrix, letting us factor out our
normalization to find

$\displaystyle\bm{C}_{ab,cd}=\mathcal{N}_{ab}\>\mathcal{N}_{cd}\big{[}$
$\displaystyle\left<\bm{\delta t}_{a}^{T}\,\bm{Q}_{ab}\,\bm{\delta
t}_{b}\;\bm{\delta t}_{c}^{T}\,\bm{Q}_{cd}\,\bm{\delta t}_{d}\right>$ (74)
$\displaystyle-\left<\bm{\delta t}_{a}^{T}\,\bm{Q}_{ab}\,\bm{\delta
t}_{b}\right>\left<\bm{\delta t}_{c}^{T}\,\bm{Q}_{cd}\,\bm{\delta
t}_{d}\right>\big{)}.$

To solve this, we will work term by term. Solving for the first term in this
equation requires much of the same process as was done in Appendix B1 and
requires us to move to index notation,

$\displaystyle\left<\bm{\delta t}_{a}^{T}\,\bm{Q}_{ab}\,\bm{\delta
t}_{b}\;\bm{\delta t}_{c}^{T}\,\bm{Q}_{cd}\,\bm{\delta t}_{d}\right>=$ (75)
$\displaystyle\hskip
56.9055pt\sum_{ijkl}\bm{Q}_{ab,ij}\,\bm{Q}_{cd,kl}\left<\bm{\delta
t}_{a,i}\,\bm{\delta t}_{b,j}\,\bm{\delta t}_{c,k}\,\bm{\delta
t}_{d,l}\right>.$

We now apply relationships for the 4th moment for our zero mean, Gaussian
random variables [50]. This will turn our equation into the expectation of
pairs of TOA residuals which we simplify to cross-correlations,
$\left<\bm{\delta t}_{a,i}\,\bm{\delta t}_{b,j}\right>=\bm{S}_{ab,ij}$.
Applying these steps, we end up with

$\left<\bm{\delta t}_{a}^{T}\,\bm{Q}_{ab}\,\bm{\delta t}_{b}\;\bm{\delta
t}_{c}^{T}\,\bm{Q}_{cd}\,\bm{\delta
t}_{d}\right>=\sum_{ijkl}\bm{Q}_{ab,ij}\,\bm{Q}_{cd,kl}\\\
\left[\bm{S}_{ab,ij}\,\bm{S}_{cd,kl}+\bm{S}_{ac,ik}\,\bm{S}_{bd,jl}+\bm{S}_{ad,il}\,\bm{S}_{bc,jk}\right].$
(76)

We can then distribute the $\bm{Q}$ quantities, and apply the following
identities, $\bm{S}_{ab,ij}=\bm{S}_{ba,ji}$, and
$\bm{Q}_{ab,ij}=\bm{Q}_{ba,ji}$, to find traces of matrix products. The
resulting set of traces create

$\displaystyle\left<\bm{\delta t}_{a}^{T}\,\bm{Q}_{ab}\,\bm{\delta
t}_{b}\;\bm{\delta t}_{c}^{T}\,\bm{Q}_{cd}\,\bm{\delta t}_{d}\right>$
$\displaystyle=\mathrm{tr}\left[\bm{S}_{ab}\bm{Q}_{ba}\right]\mathrm{tr}\left[\bm{S}_{cd}\bm{Q}_{dc}\right]$
(77)
$\displaystyle+\mathrm{tr}\left[\bm{S}_{ac}\bm{Q}_{cd}\bm{S}_{db}\bm{Q}_{ba}\right]$
$\displaystyle+\mathrm{tr}\left[\bm{S}_{ad}\bm{Q}_{dc}\bm{S}_{cb}\bm{Q}_{ba}\right].$

Focusing now on the second term from Equation 74, we can use the same
methodology as we did in subsection A1. After adding traces, we cycle the
matrices to then evaluate the expectation values,

$\left<\bm{\delta t}_{a}^{T}\,\bm{Q}_{ab}\,\bm{\delta
t}_{b}\right>\left<\bm{\delta t}_{c}^{T}\,\bm{Q}_{cd}\,\bm{\delta
t}_{d}\right>=\mathrm{tr}\left[\bm{S}_{ba}\,\bm{Q}_{ab}\right]\mathrm{tr}\left[\bm{S}_{dc}\,\bm{Q}_{cd}\right].$
(78)

Note that this term cancels out part of the first, leaving us with a more
compact representation,

$\displaystyle\bm{C}_{ab,cd}=\mathcal{N}_{ab}\>\mathcal{N}_{cd}\big{(}$
$\displaystyle\mathrm{tr}\left[\bm{S}_{ac}\bm{Q}_{cd}\bm{S}_{db}\bm{Q}_{ba}\right]$
(79)
$\displaystyle+\mathrm{tr}\left[\bm{S}_{ad}\bm{Q}_{dc}\bm{S}_{cb}\bm{Q}_{ba}\right]\big{)}.$

In this form, this equation is the same as the one found in Equation $69$ of
Section IV.B of Johnson  _et al._ [22]. However, implementing pulsar-pair
covariance in this form will be computationally expensive due to the
$\bm{Q}_{ab}$ matrix being $(N_{\mathrm{TOA},a}\times N_{\mathrm{TOA},b})$.
Instead we can adopt the rank-reduced matrices to compress our matrix
representations. It also becomes more helpful to identify the unique cases of
$\bm{C}_{ab,cd}$. These three cases are: $0$ matching pulsars
($\bm{C}_{ab,cd}$), $1$ matching pulsar ($\bm{C}_{ab,ac}$), and $2$ matching
pulsars ($\bm{C}_{ab,ab}$).

####  C1.1 0-match case

For the $0$-match case, $\bm{C}_{ab,cd}$, there are no simplifications that we
can make from Equation 79 before substituting for the temporary $\bm{Q}_{ab}$
matrices and the cross-correlation matrices
$\bm{S}_{ab}=A^{2}_{\mathrm{gw}}\Gamma_{ab}\,\bm{F}_{a}\hat{\phi}\bm{F}^{T}_{b}$.
This leaves us with

$\displaystyle\bm{C}_{ab,cd}=\mathcal{N}_{ab}\>\mathcal{N}_{cd}\bigg{(}$
$\displaystyle
A^{4}_{\mathrm{gw}}\,\Gamma_{ac}\Gamma_{bd}\,\mathrm{tr}\left[\bm{F}_{a}\hat{\phi}\bm{F}^{T}_{c}\bm{P}^{-1}_{c}\bm{F}_{c}\hat{\phi}\bm{F}^{T}_{d}\bm{P}^{-1}_{d}\bm{F}_{d}\hat{\phi}\bm{F}^{T}_{b}\bm{P}^{-1}_{b}\bm{F}_{b}\hat{\phi}\bm{F}^{T}_{a}\bm{P}^{-1}_{a}\right]$
(80)
$\displaystyle+A^{4}_{\mathrm{gw}}\,\Gamma_{ad}\Gamma_{bc}\,\mathrm{tr}\left[\bm{F}_{a}\hat{\phi}\bm{F}^{T}_{d}\bm{P}^{-1}_{d}\bm{F}_{d}\hat{\phi}\bm{F}^{T}_{c}\bm{P}^{-1}_{c}\bm{F}_{c}\hat{\phi}\bm{F}^{T}_{b}\bm{P}^{-1}_{b}\bm{F}_{b}\hat{\phi}\bm{F}^{T}_{a}\bm{P}^{-1}_{a}\right]\bigg{)}.$

In this form, we can cycle our matrices inside the trace and identify the
rank-reduced matrix product
$\bm{Z}_{a}=\bm{F}^{T}_{a}\bm{P}^{T}_{a}\bm{F}_{a}$. This lets us simplify our
equation into

$\bm{C}_{ab,cd}=\mathcal{N}_{ab}\,\mathcal{N}_{cd}\bigg{(}A^{4}_{\mathrm{gw}}\,\Gamma_{ac}\Gamma_{bd}\,\mathrm{tr}\left[\bm{Z}_{a}\hat{\phi}\bm{Z}_{c}\hat{\phi}\bm{Z}_{d}\hat{\phi}\bm{Z}_{b}\hat{\phi}\right]+A^{4}_{\mathrm{gw}}\,\Gamma_{ad}\Gamma_{bc}\,\mathrm{tr}\left[\bm{Z}_{a}\hat{\phi}\bm{Z}_{d}\hat{\phi}\bm{Z}_{c}\hat{\phi}\bm{Z}_{b}\hat{\phi}\right]\bigg{)}.$
(81)

####  C1.2 1-match case

For the $1$-match case, $\bm{C}_{ab,ac}$, we can make an additional
simplification to Equation 79 by replacing $\bm{S}_{aa}$ with $\bm{P}_{a}$,
with $\bm{P}_{a}$ potentially containing a contribution from the GWB. This
leaves us with a slightly compressed version,

$\displaystyle\bm{C}_{ab,ac}=\mathcal{N}_{ab}\,\mathcal{N}_{ac}\big{(}$
$\displaystyle\mathrm{tr}\left[\bm{P}_{a}\bm{Q}_{ac}\bm{S}_{cb}\bm{Q}_{ba}\right]$
(82)
$\displaystyle+\mathrm{tr}\left[\bm{S}_{ac}\bm{Q}_{ca}\bm{S}_{ab}\bm{Q}_{ba}\right]\big{)}.$

Following the same strategy as the $0$-match case, we substitute for the
$\bm{Q}_{ab}$ and $\bm{S}_{ab}$ matrices, cycle the traces, and identify the
$\bm{Z}_{a}$ quantities. Additionally we simplify
$\bm{P}_{a}\bm{P}^{-1}_{a}=I$. This leaves

$\displaystyle\bm{C}_{ab,ac}=\mathcal{N}_{ab}\,\mathcal{N}_{ac}\bigg{(}$
$\displaystyle
A^{2}_{\mathrm{gw}}\,\Gamma_{cb}\mathrm{tr}\left[\bm{Z}_{b}\hat{\phi}\bm{Z}_{a}\hat{\phi}\bm{Z}_{c}\hat{\phi}\right]$
(83)
$\displaystyle+A_{\mathrm{gw}}^{4}\,\Gamma_{ac}\Gamma_{ab}\mathrm{tr}\left[\bm{Z}_{b}\hat{\phi}\bm{Z}_{a}\hat{\phi}\bm{Z}_{c}\hat{\phi}\bm{Z}_{a}\hat{\phi}\right]\bigg{)}.$

Note that, despite appearing as though there are four unique $1$-match
scenarios ($\bm{C}_{ab,ac}$, $\bm{C}_{ab,bc}$, $\bm{C}_{ab,ca}$, and
$\bm{C}_{ab,cb}$), they are all mathematically equivalent through index
swapping and transpose identities.

####  C1.3 2-match case

Finally, in the $2$-match case, $\bm{C}_{ab,ab}$, we further simplify Equation
79 to

$\displaystyle\bm{C}_{ab,ab}=\mathcal{N}^{2}_{ab}\big{(}$
$\displaystyle\mathrm{tr}\left[\bm{P}_{a}\bm{Q}_{cd}\bm{P}_{b}\bm{Q}_{ba}\right]$
(84)
$\displaystyle+\mathrm{tr}\left[\bm{S}_{ab}\bm{Q}_{dc}\bm{S}_{ab}\bm{Q}_{ba}\right]\big{)}.$

Again, applying the strategies to rank-reduce, we find more products of
matrices with their inverses which allows us to simplify to,

$\displaystyle\bm{C}_{ab,ab}=\mathcal{N}^{2}_{ab}\bigg{(}$
$\displaystyle\mathrm{tr}\left[\bm{Z}_{b}\hat{\phi}\bm{Z}_{a}\hat{\phi}\right]$
(85)
$\displaystyle+\Gamma_{ab}^{2}\,A_{\mathrm{gw}}^{4}\mathrm{tr}\left[\bm{Z}_{b}\hat{\phi}\bm{Z}_{a}\hat{\phi}\bm{Z}_{b}\hat{\phi}\bm{Z}_{a}\hat{\phi}\right]\bigg{)}.$

####  C1.4 Equivalence with the traditional optimal statistic

Interestingly, we can also see why the traditional OS is effective in the weak
signal regime through these calculations. If one assumes that the effect of
the GWB is negligible, then we can set $A^{2}_{\mathrm{gw}}=0$. This zeroes
the covariance matrix for all off-diagonals and leaves the remaining diagonal
to be

$\bm{C}_{ab,ab}=\mathcal{N}^{2}_{ab}\left(\mathrm{tr}\left[\bm{Z}_{b}\hat{\phi}\bm{Z}_{a}\hat{\phi}\right]\right)=\mathrm{tr}\left[\bm{Z}_{b}\hat{\phi}\bm{Z}_{a}\hat{\phi}\right]^{-1},$
(86)

which is equivalent to the variance of the traditional OS in Equation 19.

###  C2 Per-frequency optimal statistic

Deriving pair-covariance with the PFOS is nearly identical to that of the
traditional OS, with the largest difference being the $\phi$ quantities. We
start by writing our covariance matrix as

$\bm{C}_{ab,cd}=\left<\rho_{ab}(f_{k})\rho_{cd}(f_{k})\right>-\left<\rho_{ab}(f_{k})\right>\left<\rho_{cd}(f_{k})\right>.$
(87)

We then define our per-frequency cross-correlation estimators as

$\rho_{ab}(f_{k})=\mathcal{N}_{ab}(f_{k})\bm{\delta
t}^{T}_{a}\bm{P}_{a}^{-1}\bm{F}_{a}\tilde{\phi}(f_{k})\bm{F}^{T}_{b}\bm{P}^{-1}_{b}\bm{\delta
t}_{b},$ (88)

where
$\mathcal{N}_{ab}(f_{k})=\mathrm{tr}\left[\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{b}\Phi(f_{k})\right]^{-1}$.
If we then use a temporary matrix quantity
$\bm{Q}_{ab}(f_{k})=\bm{P}_{a}^{-1}\bm{F}_{a}\tilde{\phi}(f_{k})\bm{F}^{T}_{b}\bm{P}^{-1}_{b}$,
we can now write our estimator in a format identical to the traditional OS
case,

$\rho_{ab}(f_{k})=\mathcal{N}_{ab}(f_{k})\bm{\delta
t}^{T}_{a}\bm{Q}_{ab}(f_{k})\bm{\delta t}_{b}.$ (89)

With this format, we use the same strategies and methods to solve for the
covariance in each case. The difference come from the way we define $\phi$. As
with the PFOS, we decompose
$\bm{S}_{ab}=S(f_{k})\,\Gamma_{ab}\,\bm{F}_{a}\Phi(f_{k})\bm{F}^{T}_{b}$,
allowing for a variable spectral shape. Running through the same processes, we
get expressions for the following cases.

####  C2.1 0-match case

$\displaystyle\bm{C}_{ab,cd}(f_{k})=\mathcal{N}_{ab}(f_{k})\,\mathcal{N}_{cd}(f_{k})\bigg{(}$
$\displaystyle\Gamma_{ac}\Gamma_{bd}\,S(f_{k})^{2}\,\mathrm{tr}\left[\bm{Z}_{b}\tilde{\phi}(f_{k})\bm{Z}_{a}\Phi(f_{k})\bm{Z}_{c}\tilde{\phi}(f_{k})\bm{Z}_{d}\Phi(f_{k})\right]$
(90)
$\displaystyle+\Gamma_{ad}\Gamma_{bc}\,S(f_{k})^{2}\,\mathrm{tr}\left[\bm{Z}_{b}\tilde{\phi}(f_{k})\bm{Z}_{a}\Phi(f_{k})\bm{Z}_{d}\tilde{\phi}(f_{k})\bm{Z}_{c}\Phi(f_{k})\right]\bigg{)}$

####  C2.2 1-match case

$\displaystyle\bm{C}_{ab,ac}(f_{k})=\mathcal{N}_{ab}(f_{k})\,\mathcal{N}_{ac}(f_{k})\bigg{(}$
$\displaystyle\Gamma_{bc}\,S(f_{k})\,\mathrm{tr}\left[\bm{Z}_{b}\tilde{\phi}(f_{k})\bm{Z}_{a}\tilde{\phi}(f_{k})\bm{Z}_{c}\Phi(f_{k})\right]$
(91)
$\displaystyle+\Gamma_{ac}\Gamma_{ab}\,S(f_{k})^{2}\,\mathrm{tr}\left[\bm{Z}_{b}\tilde{\phi}(f_{k})\bm{Z}_{a}\Phi(f_{k})\bm{Z}_{c}\tilde{\phi}(f_{k})\bm{Z}_{a}\Phi(f_{k})\right]\bigg{)}$

####  C2.3 2-match case

$\displaystyle\bm{C}_{ab,ab}(f_{k})=\mathcal{N}_{ab}(f_{k})^{2}\bigg{(}$
$\displaystyle\mathrm{tr}\left[\bm{Z}_{b}\tilde{\phi}(f_{k})\bm{Z}_{a}\tilde{\phi}(f_{k})\right]$
(92)
$\displaystyle+\Gamma_{ab}^{2}\,S(f_{k})^{2}\,\mathrm{tr}\left[\bm{Z}_{b}\tilde{\phi}(f_{k})\bm{Z}_{a}\Phi(f_{k})\bm{Z}_{b}\tilde{\phi}(f_{k})\bm{Z}_{a}\Phi(f_{k})\right]\bigg{)}$