# Consistency of the Parkes Pulsar Timing Array Signal with a Nanohertz
Gravitational-wave Background

Boris Goncharov<EMAIL_ADDRESS>Gran Sasso Science Institute (GSSI),
I-67100 L’Aquila, Italy  INFN, Laboratori Nazionali del Gran Sasso, I-67100
Assergi, Italy  Eric Thrane School of Physics and Astronomy, Monash
University, Clayton, VIC 3800, Australia  ARC Centre of Excellence for
Gravitational Wave Discovery  Ryan M. Shannon Centre for Astrophysics and
Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC
3122, Australia  ARC Centre of Excellence for Gravitational Wave Discovery
Jan Harms Gran Sasso Science Institute (GSSI), I-67100 L’Aquila, Italy  INFN,
Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy  N. D. Ramesh Bhat
International Centre for Radio Astronomy Research, Curtin University, Bentley,
WA 6102, Australia  George Hobbs  Australia Telescope National Facility,
CSIRO, Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia  Matthew
Kerr Space Science Division, Naval Research Laboratory, Washington, DC
20375-5352, USA  Richard N. Manchester Australia Telescope National Facility,
CSIRO, Space and Astronomy, PO Box 76, Epping, NSW 1710, Australia  Daniel J.
Reardon Centre for Astrophysics and Supercomputing, Swinburne University of
Technology, PO Box 218, Hawthorn, VIC 3122, Australia  ARC Centre of
Excellence for Gravitational Wave Discovery  Christopher J. Russell CSIRO
Scientific Computing, Australian Technology Park, Locked Bag 9013, Alexandria,
NSW 1435, Australia  Xing-Jiang Zhu Advanced Institute of Natural Sciences,
Beijing Normal University, Zhuhai 519087, China  Andrew Zic School of
Mathematical and Physical Sciences, and Research Centre in Astronomy,
Astrophysics and Astrophotonics, Macquarie University, NSW 2109, Australia
Australia Telescope National Facility, CSIRO, Space and Astronomy, PO Box 76,
Epping, NSW 1710, Australia

(Received April 11, 2022; Revised May 27, 2022; Accepted June 8, 2022)

###### Abstract

Pulsar timing array experiments have recently reported strong evidence for a
common-spectrum stochastic process with a strain spectral index consistent
with that expected of a nanohertz-frequency gravitational-wave background, but
with negligible yet non-zero evidence for spatial correlations required for a
definitive detection. However, it was pointed out by the Parkes Pulsar Timing
Array (PPTA) collaboration that the same models used in recent analyses
resulted in strong evidence for a common-spectrum process in simulations where
none is present. In this work, we introduce a methodology to distinguish
pulsar power spectra with the same amplitude from noise power spectra of
similar but distinct amplitudes. The former is the signature of a spatially
uncorrelated pulsar term of a nanohertz gravitational-wave background, whereas
the latter could represent ensemble pulsar noise properties. We test the
methodology on simulated data sets. We find that the reported common process
in PPTA pulsars is indeed consistent with the spectral feature of a pulsar
term. We recommend this methodology as one of the validity tests that the real
astrophysical and cosmological backgrounds should pass, as well as for
inferences about the spatially uncorrelated component of the background.

Gravitational waves, Millisecond pulsars, Pulsar timing method, Bayesian
statistics

††journal: ApJ Letters††software: pypolychord (Handley et al., 2015),
ptmcmcsampler (Ellis & van Haasteren, 2019), enterprise (Ellis et al., 2019),
bilby (Ashton et al., 2019), github.com/bvgoncharov/enterprise_warp,
github.com/bvgoncharov/ppta_dr2_noise_analysis.

##  1 Introduction

Pulsar timing array (PTA) experiments pursue the goal of detecting nanohertz-
frequency gravitational waves through temporal and spatial cross-correlation
of pulse arrival times from millisecond radio pulsars. The primary target
sources of such signals are coalescing supermassive binary black holes
separated by less than $\sim 0.1$ pc. Nanohertz gravitational waves produce
correlations in the timing residuals between the measured arrival times and
the arrival times predicted by the deterministic pulsar timing models (Edwards
et al., 2006). Recent searches for the stochastic gravitational wave
background by NANOGrav111The North American Nanohertz Observatory for
Gravitational Waves (McLaughlin, 2013) (Arzoumanian et al., 2020), PPTA222The
Parkes Pulsar Timing Array (Manchester et al., 2013) (Goncharov et al.,
2021a), EPTA333The European Pulsar Timing Array (Desvignes et al., 2016) (Chen
et al., 2021) and IPTA444The International Pulsar Timing Array (Verbiest et
al., 2016), a consortium of all pulsar timing arrays across the world
(Antoniadis et al., 2022) reported evidence for the “common-spectrum process”,
the same power-law component in Fourier spectra of timing residuals. The
spatial correlations necessary to claim a detection originate from the so-
called Earth term of the gravitational wave background (Hellings & Downs,
1983). The pulsar term of the background arises from the passage of
gravitational waves near the pulsars and only manifests as a spatially
uncorrelated process with the same spectrum of temporal correlations in all of
the pulsars. It is expected that evidence for the common-spectrum process, to
which both of the terms contribute, would precede the detection of the
gravitational wave background (Romano et al., 2021; Pol et al., 2021).
Goncharov et al. (2021a), on the other hand, pointed out that the methodology
employed by Arzoumanian et al. (2020) does not allow to distinguish between
common and similar noise processes in pulsars. So, it is unclear if the recent
searches have detected a bona-fide gravitational wave background.

One might ask why do we need to precisely determine the degree of similarity
in pulsar spectra if the detection mostly depends on spatial correlations
anyway? There are several reasons. Noise processes with similar and not common
noise spectra could arise from e.g. spin noise, stochastic irregularities in
rotation of the pulsars (Shannon & Cordes, 2010). In fact, the spin noise
model discussed in Meyers et al. (2021) suggests a spectral index of timing
residuals to be $\gamma=4$ (where the power spectral density is modelled as
$P(f)\propto f^{-\gamma}$), which would be difficult to distinguish from the
nanohertz gravitational wave background expected from binary supermassive
black holes, $\gamma=13/3$, with current PTAs (see e.g. Figure 13 in Renzini
et al., 2022, to inspect measurement uncertainties). Empirical models for spin
noise in millisecond pulsars predict timing noise having spectral indices and
amplitudes similar to that expected of the gravitational wave background. If
the reported signal in PTAs is not common and thus not of a gravitational-wave
origin, it could have interesting implications for spin noise and hence for
neutron star physics (Melatos & Link, 2014). Furthermore, the amplitude of the
common-spectrum process is in tension with several predictions for the
stochastic background amplitude. Recent work by Izquierdo-Villalba et al.
(2022) suggests that the gravitational wave background of the same strain
amplitude as the common-spectrum process is challenging to produce in theory
given the constraints from the quasar bolometric luminosity functions or the
local black hole mass function. Casey-Clyde et al. (2022) find the local
number density of supermassive binary black holes inferred from the amplitude
of the common-spectrum process, assuming it to originate from the stochastic
background, to be five times larger than theoretical predictions. So, even if
PTAs are observing hints of a gravitational wave background, inferences based
on the common-spectrum process might be contaminated by pulsar-intrinsic
noise, and it is important to clarify to what degree it is true. Such
inferences can be promising because the constraints on the background
amplitude from inter-pulsar correlations lag behind the ones based on spatial
auto-correlations for most pulsar timing arrays (Pol et al., 2021).

Whereas previous analyses of time-correlated noise in pulsars were based on
identifying (e.g. Lentati et al., 2016; Goncharov et al., 2021b) or modelling
(Caballero et al., 2016; Goncharov et al., 2020; Chalumeau et al., 2022) noise
power spectra, here we focus on modelling Bayesian priors using the second
PPTA data release (PPTA-DR2, Kerr et al., 2020). The methodology is described
in Section 2. In Section 3 we outline the results and we summarise the
conclusions in Section 4.

##  2 Methodology

###  2.1 Inference of the common-spectrum process

The timing residuals comprise a number of stochastic processes. Those that are
not due to gravitational waves are considered as noise. Temporally correlated
processes are called red, whereas processes without temporal correlations are
referred to as white. The power spectral density of red processes is usually
assumed to be a power law:

$P(f|A,\gamma)=\frac{A^{2}}{12\pi^{2}}\bigg{(}\frac{f}{\text{yr}^{-1}}\bigg{)}^{-\gamma}\text{yr}^{3},$
(1)

where the amplitude $A$ is in the units of strain amplitude of the stochastic
gravitational wave background at $f=\text{yr}^{-1}$ and $-\gamma$ is a
spectral index. The frequency in Fourier spectra of fluctuations in timing
residuals and concurrently a frequency of a gravitational wave that would have
induced these fluctuations is denoted $f$. To clarify, for white noise $P$ is
a constant, it does not depend on $f$. Red process spectra are modelled at
$n_{c}$ harmonically related frequencies that are multiples of the reciprocal
of the observation span. Values of $n_{c}$ as well as priors for $A$ and
$\gamma$ and the list of noise terms found in PPTA DR2 were published in
Goncharov et al. (2021a). Additional details on the noise models in PPTA DR2
are outlined by Goncharov et al. (2021b). All noise, timing model parameters
and signals of interest are modelled using the multivariate Gaussian
likelihood Lentati et al. (2014); Arzoumanian et al. (2016); Taylor et al.
(2017) and Bayesian posterior sampling. Without accounting for spatial
correlations, the total PTA likelihood is a product of individual pulsar
likelihoods. In particular, for a common-spectrum process with $A_{\text{c}}$
and $\gamma_{\text{c}}$, the total PTA likelihood is:

${\cal
L}(\bm{d}|\bm{\theta},A_{\text{c}},\gamma_{\text{c}})=\prod_{k=1}^{N}{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k},A_{\text{c}},\gamma_{\text{c}}),$ (2)

where $\bm{\theta}=(\bm{\theta}_{1},...,\bm{\theta}_{N})$ are parameters of
models that describe data of individual pulsars
$\bm{d}=(\bm{d}_{1},...,\bm{d}_{N})$, including pulsar-intrinsic “spin” noise
parameters $A_{k},\gamma_{k}$ and parameters of other noise terms which are
not of interest to us for the purpose of this work. Some of these parameters
are marginalized over analytically van Haasteren et al. (2009) and others
numerically555These parameters are fit for but not presented in our results.
Posterior samples for $\bm{\theta}_{k}$ are reweighted in target and proposal
likelihoods, and the total likelihood is generally independent of them unless
otherwise specified., so $\bm{\theta}$ is omitted from the following
equations. Both for the common-spectrum process found in PPTA DR2 and for a
classical model of the stochastic gravitational wave background from circular
supermassive black hole binaries, $\gamma_{\text{c}}=13/3$. We will fix
$\gamma_{\text{c}}$ at this value throughout our analysis. We provide
additional remarks on the data and analysis in Appendix A.

###  2.2 Importance sampling for pulsar timing arrays

When it is computationally challenging to evaluate a likelihood or to include
several measurements in one likelihood, many data analyses resort to the so-
called importance sampling (Chapter 10 in Gelman et al., 1995; Payne et al.,
2019). The idea is that the analysis is first carried out assuming a proposal
distribution – a likelihood or a prior which is easy to evaluate. Next,
posterior samples obtained from this first step are used to evaluate the
target distribution – a likelihood or a prior that represents the model we are
ultimately interested in. Finally, if proposal samples are collected from
subsets of a total data set, they can be combined into a single likelihood
through the procedure known as posterior recycling (Thrane & Talbot, 2019). To
sum up, importance sampling revolves around reweighing of likelihoods and
priors.

Let us represent our signal likelihood given by Equation 2 through target
likelihoods with the common-spectrum process, ${\cal
L}(d_{k}|A_{k,j},\gamma_{k,j},A_{\text{c}},\gamma_{\text{c}})$, and the
proposal likelihood without the common-spectrum process, ${\cal
L}(d_{k}|A_{k,j},\gamma_{k,j})$. Both the proposal and the target likelihoods
will represent individual pulsars, whereas the total signal likelihood
includes contributions from all pulsars:

${\cal L}(\bm{d}|A_{\text{c}})=\prod_{k}^{N}{\cal
Z}(d_{k})\frac{1}{n_{k}}\sum_{j=1}^{n_{k}}\frac{{\cal
L}(d_{k}|A_{k,j},\gamma_{k,j},A_{\text{c}})}{{\cal
L}(d_{k}|A_{k,j},\gamma_{k,j})}.$ (3)

Here, the sum is over the $n_{k}$ fiducial posterior samples $j$ for pulsar
$k$, generated for the proposal distribution. The product is over pulsars. The
posterior samples include $(A_{k},\gamma_{k})$, the amplitude and the spectral
index of the red noise for pulsar $k$. We use these fiducial samples as the
“proposal distribution” in order to explore a more complicated likelihood (the
“target distribution”). The Bayesian evidence $\mathcal{Z}(d_{k})$ is the
integral of the (proposal) likelihood over the prior. Importance sampling is
similar to the factorized likelihood approach (Taylor et al., 2022), where the
amplitude of the common-spectrum process is obtained through the
multiplication of posterior distributions obtained in analyses of individual
pulsar data.

###  2.3 Common versus quasi-common

In the case of a common-spectrum process, which could originate from the
gravitational wave background, nature provides us with one $A_{\text{c}}$ in
all pulsars. In the language of a hierarchical Bayesian statistics (Chapter 5
in Gelman et al., 1995), $A_{\text{c},k}$ for a pulsar $k$ is drawn from a
delta function distribution, described by a hyper-parameter that determines
the position of a delta function distribution along possible values of
$A_{\text{c}}$. The standard null hypothesis is that there is no common-
spectrum process and pulsar data sets are described by individual pulsar
noise. In this work, we propose an alternative null hypothesis where
$A_{\text{c},k}$ are drawn from a Gaussian distribution described by hyper-
parameters $\sigma_{A},\mu_{A}$ and not from the delta function distribution
described by $A_{\text{c}}$. So, in Equation 3 we introduce multiple possible
values of $A_{\text{c},k}$ and marginalize over the likelihood of such a noise
process over the Gaussian666We assume a truncated Gaussian distribution
because, in reality, our measurements are truncated by uniform prior
boundaries that are used to obtain proposal posterior samples. prior
$\pi(A_{\text{qc},k}|\mu_{A},\sigma_{A})$ with hyper-parameters $\mu_{A}$ and
$\sigma_{A}$ being the mean and the standard deviation of the Gaussian prior,
respectively:

${\cal L}(\bm{d}|\mu_{A},\sigma_{A})=\prod_{k}^{N}{\cal
Z}(d_{k})\frac{1}{n_{k}}\sum_{j=1}^{n_{k}}\frac{\int{\cal
L}(d_{k}|A_{k,j},\gamma_{k,j},A_{\text{qc},k})\,\pi(A_{\text{qc},k}|\mu_{A},\sigma_{A})dA_{\text{qc},k}}{{\cal
L}(d_{k}|A_{k,j},\gamma_{k,j},A_{\text{qc},k})}.$ (4)

We provide a derivation in the Appendix C. For $\sigma_{A}=0$,
$\pi(A_{\text{qc},k}|\mu_{A},\sigma_{A})$ reduces to the delta function.
Therefore, the likelihood ${\cal L}(\bm{d}|A_{\text{c}})$ in Equation 3 is a
subset of the likelihood ${\cal L}(\bm{d}|\mu_{A},\sigma_{A})$ in Equation 4.
We refer to noise processes with $\sigma_{A}\neq 0$ as quasi-common, which
means that noise spectra in timing array pulsars are similar but not common.
The measurement of $\sigma_{A}$ and $\mu_{A}$ described in this Section is
also applicable to modelling pulsar spin noise alone, which is broadly-
distributed in $A_{k}$. Instead, here we only infer parameters on the second
common term $A_{\text{qc}}$ to map the measurements directly to the result of
the gravitational wave search with PPTA DR2 (Goncharov et al., 2021a). The
principle of the quasi-common noise model is similar to the dropout analysis
(see Figure 9 in Arzoumanian et al., 2020), it quantifies how constraints on
the common-spectrum process with one of the pulsars are consistent with
constraints from the rest of the pulsars. As Bayes factors, high dropout
factors for pulsars support the common-spectrum process, whereas low dropout
factors that approach zero illuminate pulsars that do not have the same
spectra as the others.

We calculate the integral over $A_{\text{qc},k}$ for every evaluation of the
likelihood ${\cal L}(\bm{d}|\mu_{A},\sigma_{A})$ in Equation 4 the following
way. First, we pre-compute ${\cal
L}(d_{k}|A_{k,j},\gamma_{k,j},A_{\text{qc},k})$ on a grid of $A_{\text{qc},k}$
for $n_{k}$ posterior samples for $N$ pulsars. Next, for every parameter
sample $(\mu_{A},\sigma_{A})$, we evaluate the prior for the grid of
$A_{\text{qc},k}$, multiply it by the pre-computed likelihood ${\cal
L}(d_{k}|A_{k,j},\gamma_{k,j},A_{\text{qc},k})$, and evaluate the integral
over the product numerically. Because parameter estimation is traditionally
performed on the (base-10) logarithm of the red process amplitude, in practice
we measure $\mu_{\log_{10}A}$ and $\sigma_{\log_{10}A}$ instead of $\mu_{A}$
and $\sigma_{A}$.

Figure 1: Comparison of common-spectrum and quasi-common spectrum processes in
the PPTA DR2. Left: marginalized posterior for $\sigma_{\log_{10}A}$. The
value is consistent with zero, which means that the common-spectrum hypothesis
holds for the PPTA data. Vertical dashed line is the upper limit on
$\sigma_{\log_{10}A}$ at 95-% credibility. Right: full posterior for
$\mu_{\log_{10}A}$ and $\sigma_{\log_{10}A}$. Vertical dashed line corresponds
to the measurement of $\log_{10}A$ of the common-spectrum process in PPTA DR2
(Goncharov et al., 2021a). Our inference of the common-spectrum process
amplitude is consistent with the previous measurement considering
generalization of the model and the use of likelihood reweighting. Three
closed lines correspond to the standard $1$-, $2$\- and $3$-$\sigma$ levels.
In both panels, the dash-dotted lines show the measurements performed without
three PSRs with special noise properties discussed in Section 3: J0437$-$4715,
J0711$-$6830, and J1643$-$1224\. Removing these pulsars from the analysis
increases uncertainty levels, as expected.

##  3 Results

We measure $\mu_{\log_{10}A}$ and $\sigma_{\log_{10}A}$ in PPTA DR2 Kerr et
al. (2020), with a particular emphasis on $\sigma_{\log_{10}A}$ which
distinguishes a common-spectrum process from a quasi-common noise process, as
pointed out in Section 2.3. We find Savage-Dickey (1971) natural log Bayes
factor in favor of $\sigma_{\log_{10}A}=0$ against other $\sigma_{\log_{10}A}$
within the uniform prior to be $-0.16$. This means that $\sigma_{\log_{10}A}$
is consistent with zero. We further calculate an upper limit on
$\sigma_{\log_{10}A}<0.44$ at 95% credibility. Thus, we demonstrate that the
data are more consistent with a common-spectrum process hypothesis than the
extended quasi-common process model. We measure
$\mu_{\log_{10}A}=-14.71^{+0.08}_{-0.14}$, which is consistent with
$\log_{10}A=-14.66\pm 0.07$ found in Goncharov et al. (2021a). The result of
parameter estimation is provided in Figure 1. We note that the maximum a
posteriori value of $\sigma_{\log_{10}A}$ differs from zero, which may
represent a noise fluctuation.

We checked that the result is robust to the exclusion of three pulsars with
particular noise properties: PSRs J0437$-$4715, J0711$-$6830 and
J1643$-$1224\. PSRs J1643$-$1224 and J0711$-$6830 show spin noise with unusual
spectral indices and so we consider them population outliers. These two
pulsars also do not contribute significantly to the common-spectrum process
based on the dropout analysis (Fig. 2 in Goncharov et al., 2021a), which is
why it is acceptable to exclude them. PSR J0437$-$4715 shows excess noise that
could, in principle, marginally affect measurement of the common-spectrum
process, and that does affect the inference of spatial correlations (see the
left panel of Figure 3 in Goncharov et al., 2021a).

Additionally, we tested our methodology using simulated data sets. The details
of the simulations can be found in the Appendix B. One important issue that we
identified is that it is essential to include the common-spectrum process in
the proposal likelihood, not only in the target likelihood when testing the
quasi-common noise hypothesis. If it is excluded, the measurement of
$\sigma_{\log_{10}A}$ will always peak at zero, even if data sets were
simulated with e.g. $\sigma_{\log_{10}A}=2$. This is because the proposal
likelihood does not sufficiently match the target likelihood. Non-inclusion of
the separate common-spectrum process term will prevent accumulating posterior
support in two areas of the parameter space that correspond to the spin noise
and the common-spectrum process, respectively. Nevertheless, our methodology
is very sensitive to sub-threshold red-noise contributions, which is generally
true for all Bayesian inference that combines data from multiple measurements
(e.g., Goncharov et al., 2020). In particular, even though two power law
processes – (quasi-)common-spectrum process and pulsar-intrinsic spin noise –
are not individually resolved in single pulsars, Bayesian inference with the
total data set distinguishes them when models match the simulation.

##  4 Conclusions

We have shown that a pulsar timing array data set, in our case PPTA DR2,
contains a red process that has spectra consistent with stochastic time-series
realizations with the single power-law amplitude by measuring $\sigma_{A}$
consistent with zero. Note, this conclusion is (a) made under the assumption
that the red process power-law index is $13/3$, which is supported by the PPTA
data, but (b) not yet applicable to all available PTA data sets. This
discrepancy can potentially be resolved through further data combinations, for
example by the IPTA. The identified common-spectrum process in the PPTA DR2
could therefore be the spatially-uncorrelated component of the stochastic
gravitational wave background. We expect the variance in measured pulsar
spectra from spatial correlations of the pulsar term to exceed that from
different time-series realizations of the Earth term, and thus the uncertainty
in $\sigma_{A}$ to be dominated by the pulsar term. A more detailed
investigation of the contribution of both terms to $\sigma_{A}$ is a subject
of a follow-up study.

The methodology of previous nanohertz gravitational wave searches might have
led to identifying a common-spectrum process when none is present in the data
(Goncharov et al., 2021a). Specifically, incorrect conclusions can be caused
by a mismatch between uniform noise priors on $A$ and $\gamma$ in the models
and the clustering of these parameters in real data. In particular, for the
case of a single power-law red noise term in pulsars, evidence for the common-
spectrum process disappears when the distribution of true parameters of such
noise processes matches the Bayesian priors. The demonstration of the above
two points is to be provided by A. Zic et al. (in prep.) simulations where the
apparent common-spectrum process arises in a variety of scenarios that only
contain pulsar-intrinsic spin noise, whereas the spurious emergence of spatial
correlations in such data sets is very unlikely. Determining that the
identified noise process is quasi-common could either mean that it is not the
gravitational wave background or that the assumed uniform prior distributions
for pulsar-intrinsic noise need to be replaced by more realistic priors. The
latter may inform about stochastic irregularities in rotational properties of
neutron stars and thus about neutron star physics (Melatos & Link, 2014).

Neither our test nor any other test will be able to confirm the common-
spectrum hypothesis because it is a matter of measuring the width of the
distribution, and we would need an infinitely small measurement uncertainty to
rule out all $\sigma_{A}$ other than zero. Thus, the techniques we introduce
allow to perform a consistency test and to infer common noise properties, but
cannot be used to detect the stochastic background. Once the gravitational
wave background is detectable through spatial correlations, further modelling
of the priors of pulsar-intrinsic spin noise with the techniques we outlined
will be required to disentangle these terms. Further simulation study will be
useful to test the fidelity of the presented methods under different
conditions.

We would also like to point out that it is challenging to extend the
methodology to allow $\gamma$ to vary and thus to measure $\mu_{\gamma}$ and
$\sigma_{\gamma}$ simultaneously with $\mu_{A}$ and $\sigma_{A}$ because the
accuracy of numerical integration on a grid will decline for a higher-
dimensional problem and the computational burden will significantly increase.
We foresee several other modifications of the outlined methodology. For
example, one could model the distribution of common-process amplitudes $A$
linearly (rather than logarithmically), as well as to fix the amplitude and
fit for the distribution of power-law indices. Moreover, one could apply our
approach to fitting the distribution of pulsar-intrinsic noise amplitudes and
power-law indices, which would result in larger $\sigma_{A}$ and
$\sigma_{\gamma}$. Furthermore, in case of finding stronger evidence for
quasi-common noise with $\sigma_{A}>0$, it is possible to test a range of
functional forms outlining the distribution of $A$ other than the Gaussian
distribution we have assumed.

From the perspective of future validation of the detection of the stochastic
background, the tests we propose are complementary to other work. The dropout
analysis in Arzoumanian et al. (2020) is developed to find outliers in the
common-spectrum process. Taylor et al. (2013) introduced interpolant-based
modeling for spatial correlations, which is important to ensure that pulsars
exhibit precisely the Hellings-Downs correlations and not something else.
Johnson et al. (2022) examine a bias from using a finite number of pulsars.
Taylor et al. (2022) develop a factorized likelihood approach for cross-
validation of measurements between subsets of timing array pulsars. Romero-
Shaw et al. (2022) review several other tests to spot incorrectly specified
models. Because, as outlined by Goncharov et al. (2021b), the background noise
in pulsar timing arrays is often not white and Gaussian as assumed by current
models and simulations, we suggest to study the effect of this noise on
gravitational wave searches.

## Acknowledgements

We thank Matthew Miles, Paul Baker, Stephen Taylor and Sarah Vigeland for
useful comments. This work has been carried out using the Parkes Pulsar Timing
Array, which is part of the International Pulsar Timing Array. Murriyang, the
Parkes radio telescope is part of the Australia Telescope, which is funded by
the Commonwealth Government for operation as a National Facility managed by
CSIRO. This paper includes archived data obtained through the CSIRO Data
Access Portal (data.csiro.au). BG is supported by the Italian Ministry of
Education, University and Research within the PRIN 2017 Research Program
Framework, n. 2017SYRTCN. RMS acknowledges support through Australian Research
Council Future Fellowship FT190100155. Part of this work was undertaken as
part of the Australian Research Council Centre of Excellence for Gravitational
Wave Discovery (CE170100004). Work at NRL is supported by NASA.

##  Appendix A Notes on data and analysis

We emphasise that the sources of noise, especially, the red-noise terms,
should be modelled, to ensure the correctness of the analysis. The second data
release of the PPTA contains several sources of red noise which were found and
identified in Goncharov et al. (2021b) based their attribution to telescope
observing band or system as well as on their radio frequency dependence.
Goncharov et al. (2021b) also identified non-stationary sources of noise which
affect red noise measurements if not modelled. Incorrect modelling of pulsar
timing model parameters such as spin frequency derivatives or instrumental
phase jumps can also appear as red noise. Pulsar timing models for PPTA DR2 in
coordination with the noise analysis were performed by Reardon et al. (2021).
The data set and the code is available at
github.com/bvgoncharov/ppta_dr2_noise_analysis.

Parameter estimation and Bayesian evidence evaluation for individual pulsars
is performed with nested sampling (Skilling, 2006) implemented in pypolychord
by Handley et al. (2015). White noise parameter estimation is performed with
the ptmcmcsampler (Ellis & van Haasteren, 2019). Pulsar likelihoods are
modelled using the code enterprise (Ellis et al., 2019) and linked to
pypolychord with bilby (Ashton et al., 2019), using the code available at
github.com/bvgoncharov/enterprise_warp.

##  Appendix B Validating the general quasi-common noise model

We performed a study using a simulated data set to confirm the validity of the
importance sampling method as well as the quasi-common spectrum process model,
for which the common-spectrum process is a limiting case. We note that the
common-spectrum process hypothesis for a given pulsar is described by a
likelihood from Equation 2, approximated with the importance sampling by
Equation 3. This likelihood is then generalized in Equation 4 to represent the
quasi-common spectrum process hypothesis.

We simulated a data set with 26 pulsars, as in PPTA DR2, observed for 2555
days, with the pulse arrival time errors of $\sigma_{\text{ToA}}=0.1~{}\mu s$
and the white noise with variance proportional to $\sigma_{\text{ToA}}$. We
then simulate two red noise processes described by $\log_{10}A$ and $\gamma$,
as we expect from the stochastic gravitational wave background and pulsar
noise. The simulations are based on 30 frequencies, which corresponds to the
Fourier basis to model red noise in previous analyses. One red noise process
for each pulsar is drawn from the truncated Normal distribution
$\mathcal{N}(\mu,\sigma)$ with $\mu_{\log_{10}A}=-16.3$, $\mu_{\gamma}=5$,
$\sigma_{\log_{10}A}=2$, $\sigma_{\gamma}=2$. Another red noise process, the
quasi-common noise, has a fixed $\gamma=13/3$ and $\log_{10}A$ is drawn from
the truncated Normal distribution with $\mu_{\log_{10}A}=-13.3$ and
$\sigma_{\log_{10}A}=0.5$. Both Normal distributions are truncated to the
edges of uniform priors used in Goncharov et al. (2021a) to avoid edge
effects. The values were chosen so that pulsar-intrinsic red noise was below
the detection threshold of some pulsars, as in PPTA DR2 where only 9 pulsars
of 26 show evidence for spin noise. Moreover, quasi-common noise realizations
have simulated $A$ of the order of $10^{-14}$ \- $10^{-13}$, as six
realizations of the spin noise. The same number of pulsars in PPTA DR2 have
the inferred common process amplitude and the spin noise amplitude of the same
order of magnitude. The results of parameter estimation on the simulated data
set are presented in Figure 2. Both the $\mu_{\log_{10}A}$ and the
$\sigma_{\log_{10}A}$ are consistent with the simulated values. Note, we
tested that this data set shows strong evidence for the common-spectrum
process when using the priors from Goncharov et al. (2021a), and yet with our
new generalized model we correctly infer that it is not the case because
$\sigma_{\log_{10}A}\neq 0$. More precisely, we find
$\log\mathcal{B}^{\text{CP}}_{\varnothing}>13.5$, where $\varnothing$ is the
noise-only null hypothesis that includes only white and red noise, as per the
simulations.

The noise in real data is more sophisticated than in the simulations, but it
is not necessary to represent the whole complexity of the data to demonstrate
the validity of the approach. We defer exhaustive tests of the methodology to
future work. Our method is applicable to any data set under the assumption
that physically relevant red noise processes (from interstellar propagation
effects, instrumental noise, etc.) are separated from the common-spectrum
process by including them in the models. Moreover, we also find that the
current simulation study is robust to the choice of priors for pulsar-
intrinsic noise, with the only caveat being that the proposal likelihood
should include all red noise terms in pulsars to allow parameters of both
processes to be sampled. We tested that the simulation works for a reduced
case with only one red noise process per pulsar, correctly recovering position
and width of a Gaussian distribution of pulsar red noise parameters. We also
trialed it for the case where the maximum probability density of the
distributions of pulsar spin noise amplitudes ($\mu_{\log_{10}A}=-14.3$,
$\sigma_{\log_{10}A}=1.3$) and quasi-common noise amplitudes
($\mu_{\log_{10}A}=-13.8$, $\sigma_{\log_{10}A}=0.4$) are similar and the
distributions have a broader overlap.

Figure 2: Searches for a quasi common-spectrum process in a simulated data
set. Left: marginalized posterior for $\sigma_{\log_{10}A}$. The value is
inconsistent with zero, which means that the signal that would have been
interpreted as the common-spectrum process with the standard methods is only a
quasi-common spectrum process. Based on this result, the gravitational wave
origin of the signal can be ruled out. The simulated value is represented by
the vertical dotted line and cumulative $1$-$\sigma$ credible levels are
represented by a shaded region. Right: full posterior on $\mu_{\log_{10}A}$
and $\sigma_{\log_{10}A}$. The vertical and the horizontal lines correspond to
the simulated values. Three closed lines correspond to the standard $1$-,
$2$\- and $3$-$\sigma$ levels.

##  Appendix C Derivation of the quasi-common process likelihood

Whereas Equation 2 models contributions from the common-spectrum process to be
represented by one $A_{\text{c}}$ in all pulsars, we can generalize it to
model different $A_{\text{c}}$ in pulsars. This way, assuming fixed $\gamma$,
Equation 2 transforms into

${\cal L}(\bm{d}|\bm{\theta},A_{c,1},...,A_{c,N})=\prod_{k=1}^{N}{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k},A_{c,k}).$ (C1)

Marginalizing over all possible values of $A_{c,k}$ in pulsars over the prior
distribution these parameters are drawn from, which we model as a Gaussian
distribution with $\mu$ and $\sigma$, the likelihood becomes

${\cal L}(\bm{d}|\bm{\theta},\mu,\sigma)=\prod_{k=1}^{N}{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k},\mu,\sigma)=\prod_{k=1}^{N}\int{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k},A_{c,k})\pi(A_{c,k}|\mu,\sigma)dA_{c,k}.$ (C2)

Next, let us marginalize over “nuisance” parameters $\bm{\theta}$,

${\cal L}(\bm{d}|\mu,\sigma)=\prod_{k=1}^{N}\int{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k},\mu,\sigma)\pi(\bm{\theta}_{k})d\bm{\theta}_{k}.$
(C3)

Let us then multiply the equation by unity and expand the unity via the
proposal likelihood $\mathcal{L}(\bm{d}_{k}|\bm{\theta}_{k})$,

${\cal
L}(\bm{d}|\mu,\sigma)=\prod_{k=1}^{N}\int\frac{\mathcal{L}(\bm{d}_{k}|\bm{\theta}_{k})}{\mathcal{L}(\bm{d}_{k}|\bm{\theta}_{k})}{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k},\mu,\sigma)\pi(\bm{\theta}_{k})d\bm{\theta}_{k}.$
(C4)

Representing one proposal likelihood through the posterior
$\mathcal{P}(\bm{\theta}_{k}|\bm{d}_{k})$, the evidence
$\mathcal{Z}(\bm{d}_{k})$, and the prior $\pi(\bm{\theta}_{k})$, we obtain

${\cal
L}(\bm{d}|\mu,\sigma)=\prod_{k=1}^{N}\int\mathcal{Z}(\bm{d}_{k})\mathcal{P}(\bm{\theta}_{k}|\bm{d}_{k})\bigg{(}\frac{{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k},\mu,\sigma)}{\mathcal{L}(\bm{d}_{k}|\bm{\theta}_{k})}\bigg{)}d\bm{\theta}_{k},$
(C5)

where the prior canceled out. Next, we approximate the integral for a $k$’th
pulsar with a sum over $n_{k}$ posterior samples, $\bm{\theta}_{k,j}$,
providing

${\cal
L}(\bm{d}|\mu,\sigma)=\prod_{k=1}^{N}\mathcal{Z}(\bm{d}_{k})\frac{1}{n_{k}}\sum_{j=1}^{n_{k}}\frac{{\cal
L}(\bm{d}_{k}|\bm{\theta}_{k,j},\mu,\sigma)}{\mathcal{L}(\bm{d}_{k}|\bm{\theta}_{k,j})},$
(C6)

where, omitting $\bm{\theta}_{k,j}$,

${\cal L}(\bm{d}_{k}|\mu,\sigma)=\int{\cal
L}(\bm{d}_{k}|A_{c,k})\pi(A_{c,k}|\mu,\sigma)dA_{c,k}.$ (C7)

The approximation of an integral via the sum of posterior samples is explained
for various applications in Hogg & Foreman-Mackey (2018) and Thrane & Talbot
(2019).

## References

  * Antoniadis et al. (2022) Antoniadis, J., Arzoumanian, Z., Babak, S., et al. 2022, MNRAS, 510, 4873, doi: 10.1093/mnras/stab3418
  * Arzoumanian et al. (2016) Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al. 2016, ApJ, 821, 13, doi: 10.3847/0004-637X/821/1/13
  * Arzoumanian et al. (2020) Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020, ApJ, 905, L34, doi: 10.3847/2041-8213/abd401
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, ApJS, 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Caballero et al. (2016) Caballero, R. N., Lee, K. J., Lentati, L., et al. 2016, MNRAS, 457, 4421, doi: 10.1093/mnras/stw179
  * Casey-Clyde et al. (2022) Casey-Clyde, J. A., Mingarelli, C. M. F., Greene, J. E., et al. 2022, ApJ, 924, 93, doi: 10.3847/1538-4357/ac32de
  * Chalumeau et al. (2022) Chalumeau, A., Babak, S., Petiteau, A., et al. 2022, MNRAS, 509, 5538, doi: 10.1093/mnras/stab3283
  * Chen et al. (2021) Chen, S., Caballero, R. N., Guo, Y. J., et al. 2021, MNRAS, 508, 4970, doi: 10.1093/mnras/stab2833
  * Desvignes et al. (2016) Desvignes, G., Caballero, R. N., Lentati, L., et al. 2016, MNRAS, 458, 3341, doi: 10.1093/mnras/stw483
  * Dickey (1971) Dickey, J. M. 1971, The Annals of Mathematical Statistics, 42, 204.  http://www.jstor.org/stable/2958475
  * Edwards et al. (2006) Edwards, R. T., Hobbs, G. B., & Manchester, R. N. 2006, MNRAS, 372, 1549, doi: 10.1111/j.1365-2966.2006.10870.x
  * Ellis & van Haasteren (2019) Ellis, J., & van Haasteren, R. 2019, PTMCMCSampler: Parallel tempering MCMC sampler package written in Python, Astrophysics Source Code Library, record ascl:1912.017.  http://ascl.net/1912.017
  * Ellis et al. (2019) Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T. 2019, ENTERPRISE: Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE, Astrophysics Source Code Library, record ascl:1912.015.  http://ascl.net/1912.015
  * Gelman et al. (1995) Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. 1995, Bayesian data analysis (Chapman and Hall/CRC)
  * Goncharov et al. (2020) Goncharov, B., Zhu, X.-J., & Thrane, E. 2020, MNRAS, 497, 3264, doi: 10.1093/mnras/staa2081
  * Goncharov et al. (2021a) Goncharov, B., Shannon, R. M., Reardon, D. J., et al. 2021a, ApJ, 917, L19, doi: 10.3847/2041-8213/ac17f4
  * Goncharov et al. (2021b) Goncharov, B., Reardon, D. J., Shannon, R. M., et al. 2021b, MNRAS, 502, 478, doi: 10.1093/mnras/staa3411
  * Handley et al. (2015) Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015, MNRAS, 453, 4384, doi: 10.1093/mnras/stv1911
  * Hellings & Downs (1983) Hellings, R. W., & Downs, G. S. 1983, ApJ, 265, L39, doi: 10.1086/183954
  * Hogg & Foreman-Mackey (2018) Hogg, D. W., & Foreman-Mackey, D. 2018, ApJS, 236, 11, doi: 10.3847/1538-4365/aab76e
  * Izquierdo-Villalba et al. (2022) Izquierdo-Villalba, D., Sesana, A., Bonoli, S., & Colpi, M. 2022, MNRAS, 509, 3488, doi: 10.1093/mnras/stab3239
  * Johnson et al. (2022) Johnson, A. D., Vigeland, S. J., Siemens, X., & Taylor, S. R. 2022, arXiv e-prints, arXiv:2201.10657.  https://arxiv.org/abs/2201.10657
  * Kerr et al. (2020) Kerr, M., Reardon, D. J., Hobbs, G., et al. 2020, PASA, 37, e020, doi: 10.1017/pasa.2020.11
  * Lentati et al. (2014) Lentati, L., Alexander, P., Hobson, M. P., et al. 2014, MNRAS, 437, 3004, doi: 10.1093/mnras/stt2122
  * Lentati et al. (2016) Lentati, L., Shannon, R. M., Coles, W. A., et al. 2016, MNRAS, 458, 2161, doi: 10.1093/mnras/stw395
  * Manchester et al. (2013) Manchester, R. N., Hobbs, G., Bailes, M., et al. 2013, PASA, 30, e017, doi: 10.1017/pasa.2012.017
  * McLaughlin (2013) McLaughlin, M. A. 2013, Classical and Quantum Gravity, 30, 224008, doi: 10.1088/0264-9381/30/22/224008
  * Melatos & Link (2014) Melatos, A., & Link, B. 2014, MNRAS, 437, 21, doi: 10.1093/mnras/stt1828
  * Meyers et al. (2021) Meyers, P. M., O’Neill, N. J., Melatos, A., & Evans, R. J. 2021, MNRAS, 506, 3349, doi: 10.1093/mnras/stab1952
  * Payne et al. (2019) Payne, E., Talbot, C., & Thrane, E. 2019, Phys. Rev. D, 100, 123017, doi: 10.1103/PhysRevD.100.123017
  * Pol et al. (2021) Pol, N. S., Taylor, S. R., Kelley, L. Z., et al. 2021, ApJ, 911, L34, doi: 10.3847/2041-8213/abf2c9
  * Reardon et al. (2021) Reardon, D. J., Shannon, R. M., Cameron, A. D., et al. 2021, MNRAS, 507, 2137, doi: 10.1093/mnras/stab1990
  * Renzini et al. (2022) Renzini, A. I., Goncharov, B., Jenkins, A. C., & Meyers, P. M. 2022, Galaxies, 10, 34, doi: 10.3390/galaxies10010034
  * Romano et al. (2021) Romano, J. D., Hazboun, J. S., Siemens, X., & Archibald, A. M. 2021, Phys. Rev. D, 103, 063027, doi: 10.1103/PhysRevD.103.063027
  * Romero-Shaw et al. (2022) Romero-Shaw, I. M., Thrane, E., & Lasky, P. D. 2022, arXiv e-prints, arXiv:2202.05479.  https://arxiv.org/abs/2202.05479
  * Shannon & Cordes (2010) Shannon, R. M., & Cordes, J. M. 2010, ApJ, 725, 1607, doi: 10.1088/0004-637X/725/2/1607
  * Skilling (2006) Skilling, J. 2006, Bayesian Analysis, 1, 833 , doi: 10.1214/06-BA127
  * Taylor et al. (2013) Taylor, S. R., Gair, J. R., & Lentati, L. 2013, Phys. Rev. D, 87, 044035, doi: 10.1103/PhysRevD.87.044035
  * Taylor et al. (2017) Taylor, S. R., Lentati, L., Babak, S., et al. 2017, Phys. Rev. D, 95, 042002, doi: 10.1103/PhysRevD.95.042002
  * Taylor et al. (2022) Taylor, S. R., Simon, J., Schult, L., Pol, N., & Lamb, W. G. 2022, arXiv e-prints, arXiv:2202.08293.  https://arxiv.org/abs/2202.08293
  * Thrane & Talbot (2019) Thrane, E., & Talbot, C. 2019, PASA, 36, e010, doi: 10.1017/pasa.2019.2
  * van Haasteren et al. (2009) van Haasteren, R., Levin, Y., McDonald, P., & Lu, T. 2009, MNRAS, 395, 1005, doi: 10.1111/j.1365-2966.2009.14590.x
  * Verbiest et al. (2016) Verbiest, J. P. W., Lentati, L., Hobbs, G., et al. 2016, MNRAS, 458, 1267, doi: 10.1093/mnras/stw347