# Measuring the Hubble Constant with Double Gravitational Wave Sources in
Pulsar Timing

Casey McGrath ,1 2 3 Daniel J. D’Orazio ,4 and Jolien Creighton ,5
1Center for Space Sciences and Technology, University of Maryland, Baltimore
County, Baltimore, MD 21250, USA
2Gravitational Astrophysics Lab, NASA/GSFC, Greenbelt, MD 20771, USA
3Center for Research and Exploration in Space Science and Technology II,
NASA/GSFC, Greenbelt, MD 20771, USA
4Niels Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17,
DK-2100 Copenhagen, Denmark
5Center for Gravitation, Cosmology and Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA
casey.d.mcgrath@nasa.govdaniel.dorazio<EMAIL_ADDRESS>0000-0002-6155-3501 0000-0002-1271-6247 0000-0003-3600-2406

(Accepted 2022 September 09. Received 2022 September 09; in original form 2022
August 12)

###### Abstract

Pulsar timing arrays (PTAs) are searching for gravitational waves from
supermassive black hole binaries (SMBHBs). Here we show how future PTAs could
use a detection of gravitational waves from individually resolved SMBHB
sources to produce a purely gravitational wave-based measurement of the Hubble
constant. This is achieved by measuring two separate distances to the same
source from the gravitational wave signal in the timing residual: the
luminosity distance $D_{L}$ through frequency evolution effects, and the
parallax distance $D_{\mathrm{par}}$ through wavefront curvature (Fresnel)
effects. We present a generalized timing residual model including these
effects in an expanding universe. Of these two distances, $D_{\mathrm{par}}$
is challenging to measure due to the pulsar distance wrapping problem, a
degeneracy in the Earth-pulsar distance and gravitational wave source
parameters that requires highly precise, sub-parsec level, pulsar distance
measurements to overcome. However, in this paper we demonstrate that combining
the knowledge of two SMBHB sources in the timing residual largely removes the
wrapping cycle degeneracy. Two sources simultaneously calibrate the PTA by
identifying the distances to the pulsars, which is useful in its own right,
and allow recovery of the source luminosity and parallax distances which
results in a measurement of the Hubble constant. We find that, with optimistic
PTAs in the era of the Square Kilometre Array, two fortuitous SMBHB sources
within a few hundred Mpc could be used to measure the Hubble constant with a
relative uncertainty on the order of 10 per cent.

###### keywords:

gravitational waves – quasars: supermassive black holes – pulsars: general –
cosmological parameters

††pubyear: 2022††pagerange: Measuring the Hubble Constant with Double
Gravitational Wave Sources in Pulsar Timing–D

##  1 Introduction

Pulsar timing arrays (PTAs) are currently searching for gravitational waves
with frequencies $\mathcal{O}$(1 to 100 nHz) produced by supermassive black
hole binaries (SMBHBs) at the centers of coalescing galaxies. These
experiments time highly regular millisecond pulsars across our Galaxy and look
for irregularities with the arrival times of the pulsar’s pulses. A continuous
gravitational wave source is expected to cause these arrival times to
periodically drift in and out of synchronization with a reference clock. While
the first gravitational wave detection using this experimental method may be a
stochastic background of unresolved sources (see, however, Kelley et al.,
2018), these experiments are expected to become more sensitive to the point
where loud individual sources will be resolved and parameter estimation will
be able to extract the source parameters from the data.

Of notable interest is the ability of a gravitational wave-based experiment to
measure cosmological parameters such as the Hubble constant $H_{0}$. These new
gravitational wave-based measurements are important in helping us to resolve
the current tension between experimental measurements of $H_{0}$ (The LIGO,
Virgo, 1M2H, Dark Energy Camera GW-EM, & DES Collaborations, 2017; The LIGO,
Virgo, & KAGRA Collaborations, 2021; Feeney et al., 2019, and references
therein). In this study we are motivated by the recent work of D’Orazio & Loeb
(2021) and McGrath & Creighton (2021) (hereafter DL21 and MC21). DL21
demonstrated how a measurement of the Hubble constant could be made from a
purely gravitational wave-based method, by measuring the source’s luminosity
distance $D_{L}$ and it’s comoving distance $D_{c}$. In this approach $D_{L}$
is recovered from frequency evolution in the timing residual signal, and
$D_{c}$ comes from probing the curvature of the wavefront across the Earth-
pulsar baseline of the PTA experiment. In this paper we show that more
generally, the distance measured from the curvature of the wavefront is the
“parallax distance” $D_{\mathrm{par}}$, which is equivalent to $D_{c}$ in a
flat universe. MC21 generalized the current gravitational wave timing residual
models by classifying these wavefront curvature effects into the “Fresnel
regime,” and they studied how well this new distance parameter could be
recovered for different PTA constraints and source parameters. Both of these
studies were strongly motivated by and synthesized the previous work of Deng &
Finn (2011) and Corbin & Cornish (2010).

For comparison, the “standard sirens” approach to measuring $H_{0}$ is a
hybrid technique requiring a source observable via both gravitational wave and
electromagnetic messengers (Schutz, 1986; Holz & Hughes, 2005). Here the
luminosity distance $D_{L}$ is measured via gravitational waves from a
chirping source, while the redshift $z$ is measured from the host galaxy
through electromagnetic observations. It has been showed that future PTAs may
also be able to make measurements of $H_{0}$ through this technique (Wang et
al., 2022). Combining these measurements of $D_{L}$ and $z$ allows an
inference of $H_{0}$ to be made. As an example of a purely gravitational wave-
based measurement of $H_{0}$, gravitational wave signals from binary neutron
star or neutron star black hole mergers could be used to obtain both
measurements of the source’s luminosity distance $D_{L}$ and redshift $z$
(Messenger & Read, 2012; Ghosh et al., 2022; Shiralilou et al., 2022). This
would require well constrained knowledge of the neutron star equation of state
from numerous detections, which may be possible with future Einstein Telescope
and Cosmic Explorer-era gravitational wave observatories. The approach
presented in this paper is also purely gravitational wave-based, but here we
trade the redshift measurement for a second distance measurement made from the
Fresnel wavefront curvature effects, and we do not require any intrinsic
knowledge of rest-frame source properties.

In this study we take the models developed in MC21 and further generalize them
to a cosmologically expanding universe. We apply the same Bayesian framework
and methods described in that study, in order to predict how well future PTA
experiments may be able to measure the Hubble constant. When detecting a
single SMBHB source, we find that in order to recover the parallax distance
parameter (and hence measure $H_{0}$), we require highly accurate measurements
of the distances to the pulsars in our array. But crucially, we find that when
detecting two SMBHB sources simultaneously, no prior knowledge on the pulsar
distances is required in order to recover $D_{\mathrm{par}}$ and $H_{0}$.
Multiple simultaneous continuous wave SMBHBs have not been widely investigated
within PTA research (some studies include Babak & Sesana, 2012; Petiteau et
al., 2013; Qian et al., 2022), therefore this work helps us to motivate a
particular insight that multiple sources can provide over just a single
source. The results presented in this work also demonstrate how the same
methods used to measure $H_{0}$ do so through improved measurements of the
pulsar distances, which is an important result in its own right. In this paper
we focus on the $H_{0}$ measurement, but an additional paper (blueMcGrath,
D’Orazio, & Creighton, in preparation) will provide greater details on the
pulsar distance recovery.

This paper is outlined as follows. Section 2 presents the theoretical
background and generalization of our timing residual models, which shows how
the Hubble constant enters the pulsar timing model. Section 3 explains the
methods used to estimate system parameters and $H_{0}$ from mock observations.
Section 4 presents the primary obstacle to a practical $H_{0}$ measurement,
the pulsar distance wrapping problem, and our solution. Section 5 presents the
main results, simulated measurements of the Hubble constant for different mock
PTAs. Finally, conclusions and future directions are presented in Section 6.

##  2 Theoretical Background

###  2.1 Fresnel and Frequency Evolving Regimes

To detect and extract information from gravitational waves using PTAs, one
must model the pulse-time-of-arrival deviations induced by a gravitational
wave passing through the array. Such gravitational wave-induced timing
residual models can be classified based on two physical qualities/assumptions
of the model: the shape of the incoming gravitational wavefront, and its
frequency evolution. Following MC21, Figure 1 categorizes four distinct model
regimes. Models that assume no wavefront curvature over the Earth-pulsar
baseline make the “plane-wave” assumption (I), while models which account for
the first-order curvature terms are classified as “Fresnel” models (II). Then,
if the frequency of the gravitational wave is assumed constant over the
thousands of years it takes a photon to travel from the pulsar to the Earth,
the model is “monochromatic” (A); otherwise it is in a “frequency evolving”
(B) regime.

Figure 1: A classification summary of the four gravitational wave regimes,
increasing in generality from left-to-right and top-to-bottom. The frequency
evolution regime reduces to the monochromatic regime in the large coalescence
time limit ($\Delta\tau_{c}\rightarrow\infty$), and the Fresnel regime reduces
to plane-wave regime in either the small Fresnel number limit ($F\rightarrow
0$) or the natural plane-wave limit ($\hat{r}\cdot\hat{p}\rightarrow\pm 1$).
Importantly, in a cosmologically expanding universe the frequency evolution
regimes allow for the direct measurement of the source luminosity distance
$D_{L}$, and the Fresnel regimes allow for the direct measurement of the
source parallax distance $D_{\mathrm{par}}$. In the fully general IIB model,
both of these distances can be independently measured, and therefore through
equation 7 a value of $H_{0}$ can be inferred.

The result is four distinct timing residual models, increasing in generality
from left-to-right, and top-to-bottom. While simpler models are mathematically
and computationally simpler and less expensive to calculate, they may produce
inaccurate predictions of the timing residual if the source is producing
gravitational waves which fall into a more complicated regime. The appropriate
limits wherein one model approaches another are governed by the coalescence
time $\Delta\tau_{c}$ (equation 10) of the source, the Fresnel number $F$
(equation 5), and the relative source-pulsar-Earth orientation
($\hat{r}\cdot\hat{p}$; see equation 11). Therefore considering the relative
size of these quantities is crucial in deciding when one model’s prediction
will fail in comparison to another.

All four models are provided and studied explicitly in MC21, which critically
assumed a flat, static universe. However, DL21 show that in a cosmologically
expanding universe, two distinct distances appear in the gravitational wave-
induced pulsar timing residual model – the luminosity distance $D_{L}$ and the
comoving coordinate distance $D_{c}$. Accounting for frequency evolution
effects (row B, Figure 1) in the timing model allows the direct measurement of
$D_{L}$, while accounting for Fresnel effects (column II, Figure 1) allows the
direct measurement of the parallax distance $D_{\mathrm{par}}$, which reduces
to $D_{c}$ in a flat universe. The goal of this paper is to remove the
assumption that the universe is static on cosmological scales, re-derive the
formulae of MC21 in a cosmologically expanding universe, and study
implications for recovering $H_{0}$. We provide a general list of assumptions
made in this work in Appendix C.

###  2.2 Incorporating Cosmological Expansion and the Hubble Constant

In a flat static universe we have eight gravitational wave source parameters
for a circular binary SMBHB system in the Newtonian regime,
$\\{R,\theta,\phi,\iota,\psi,\theta_{0},\mathcal{M},\omega_{0}\\}$. These are
the Earth-source comoving coordinate distance (for a flat static universe),
sky angles, orientation Euler angles (inclination and polarization), initial
phase (which also can be interpreted as the third orientation Euler angle),
system chirp mass
$\mathcal{M}\equiv\left(m_{1}m_{2}\right)^{3/5}/\left(m_{1}+m_{2}\right)^{1/5}$
(for component masses $m_{1}$ and $m_{2}$), and initial (angular) orbital
frequency. The pulsar parameters are $\\{L,\theta_{p},\phi_{p}\\}$, which are
the Earth-pulsar distance, and pulsar sky angles.

In our cosmological framework we divide space-time into two types of frames: a
“global cosmological frame” where the background metric is
Friedmann–Lemaître–Robertson–Walker (FLRW), and “local frames” where the
background metric is Minkowski (see Appendix C, assumption 3). The
gravitational waves (as described by the metric perturbation below) are
generated in the local source frame, propagate over the cosmological frame,
and reach the Milky Way Galaxy where we assume they are in the local observer
frame. Therefore the cosmological effects of an expanding universe only need
to be considered during the gravitational wave’s propagation between the local
source and observer frames.

We begin by assuming that the FLRW metric describes our global space-time
background between the Earth and our gravitational wave source:

$\displaystyle\mathrm{d}s^{2}$
$\displaystyle=-c^{2}\mathrm{d}t^{2}+a^{2}(t)\left[\frac{\mathrm{d}r^{2}}{1-kr^{2}}+r^{2}\mathrm{d}\Omega^{2}\right],$
$\displaystyle=a^{2}(\eta)\left[-c^{2}\mathrm{d}\eta^{2}+\mathrm{d}\chi^{2}+r(\chi)^{2}\mathrm{d}\Omega^{2}\right]$
(1)

where
$\mathrm{d}\Omega^{2}\equiv\mathrm{d}\theta^{2}+\sin^{2}(\theta)\mathrm{d}\phi^{2}$,
$a$ is the universe scale factor, and $k$ is the space-time curvature
constant. The coordinates $(t,r,\theta,\phi)$ are the comoving coordinates,
while the version shown in the second line is made through the change to
conformal time $\mathrm{d}\eta\equiv\frac{\mathrm{d}t}{a}$ (i.e.
$\eta\equiv\int^{t}\frac{\mathrm{d}t^{\prime}}{a}$) and the spatial coordinate
transformation:

$r(\chi)\equiv\begin{cases}\frac{1}{\sqrt{k}}\sin\left(\sqrt{k}\chi\right)&,\quad
k>0\\\ \chi&,\quad k=0\\\
\frac{1}{\sqrt{|k|}}\sinh\left(\sqrt{|k|}\chi\right)&,\quad k<0.\end{cases}$
(2)

Here we are using the convention that the curvature parameter carries units of
$\left[\mathrm{length}^{-2}\right]$, where $k>0$, $k=0$, or $k<0$ for closed,
flat, and open spatial geometries, respectively. This implies that our time-
dependent scale factor $a$ is unitless. We also choose to normalize
$a(t_{0})\equiv 1$ at the present day time $t_{0}$, but we will often write it
explicitly in our derivations here. Therefore, the Gaussian curvature of a
spatial slice of our universe at the present time is $k$.

In this framework, without an independent measurement of the source’s redshift
(either from electromagnetic observations, or as we will explain, from effects
of the curvature of the wavefront itself), the redshift parameter $z$ cannot
be disentangled from the source rest frame $R$ (now the transverse comoving
distance; see Appendix A), $\mathcal{M}$, and $\omega_{0}$ parameters (see for
example, Maggiore, 2008). The result is that we can only observe the
luminosity distance $D_{L}=(1+z)R$, a redshifted system chirp mass
$\mathcal{M}_{\mathrm{obs}}\equiv(1+z)\mathcal{M}$, and a redshifted orbital
frequency parameter $\omega_{0,\mathrm{obs}}=\frac{\omega_{0}}{(1+z)}$.

As the gravitational wave propagates through the FLRW universe, the wavefront
is traced by the retarded time. Therefore in order to understand time delay
effects of the wave arriving at our pulsar compared to the time delay effects
of the wave arriving at the Earth, we must calculate the expression for the
retarded time of the wave at a field position $\vec{x}$ for a source at the
location $\vec{x}^{\prime}$ (see for example, Caldwell, 1993). Crucially, if
we make the assumption that the $a$ does not evolve appreciably over the time
it takes a photon to travel from our pulsar to the Earth (see assumption 2),
then
$\eta_{\mathrm{ret}}-\eta\approx\frac{1}{a(t_{0})}\left(t_{\mathrm{ret}}-t\right)$
and we can write:

$\displaystyle\eta_{\mathrm{ret}}=\eta-\frac{\left|\vec{x}^{\prime}-\vec{x}\right|}{c}$
$\displaystyle\therefore\quad t_{\mathrm{ret}}\approx
t-a(t_{0})\frac{\left|\vec{x}^{\prime}-\vec{x}\right|}{c}.$ (3)

Next we must calculate $\left|\vec{x}^{\prime}-\vec{x}\right|$ in the three
possible spatial geometries of our universe. The details of this are shown in
Appendix B, and the result is that regardless of the background curvature of
the spatial slices of our cosmology, the generalized expression for the
retarded time evaluated with the field position set to the pulsar’s location
is:

$t_{\mathrm{ret}}\left(t,\vec{x}=L\hat{p}\right)\ \ \approx\ \ t\ \
-\underbrace{\frac{D_{c}}{c}}_{\mathrm{``Far}\ \mathrm{Field"}}\ +\ \
\underbrace{\left(\hat{r}\cdot\hat{p}\right)a(t_{0})\frac{L}{c}}_{\mathrm{``Plane-
Wave"}}\ \ -\ \
\underbrace{\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)a^{2}(t_{0})\frac{L}{c}\frac{L}{D_{\mathrm{par}}}}_{\mathrm{``Fresnel"}}\
\ +\ \ \ldots$ (4)

The Fresnel term gives the first-order curvature of the physical wavefront,
and introduces an additional time delay to the arrival time of the wavefront
as predicted by the plane-wave approximation. A useful proxy for determining
when this term in the expansion will contribute a significant correction to
the plane-wave approximation is the Fresnel number (see MC21 and DL21 for more
discussion of this quantity):

$\displaystyle F$
$\displaystyle\equiv\frac{L^{2}}{\lambda_{\mathrm{gw,obs}}D_{\mathrm{par}}}\approx
0.0003\left(\frac{\omega_{0,\mathrm{obs}}}{1\text{
nHz}}\right)\left(\frac{L}{1\text{ kpc}}\right)^{2}\left(\frac{100\text{
Mpc}}{D_{\mathrm{par}}}\right),$ (5)
$\displaystyle\mathrm{and}\quad\lambda_{\mathrm{gw,obs}}\equiv\frac{\pi
c}{\omega_{0,\mathrm{obs}}}\approx(30.5\text{ pc})\left(\frac{1\text{
nHz}}{\omega_{0,\mathrm{obs}}}\right).$ (6)

There are two notable results that come out of the expression in equation 4.
First, we find that two cosmological distance measurements to the same source
appear in the “Fresnel” term in the expansion – the line-of-sight comoving
distance $D_{c}$, and the parallax distance $D_{\mathrm{par}}$. In practice,
however, we can only measure $D_{\mathrm{par}}$ from our timing residual
models. This is because in addition to $\frac{D_{c}}{c}$ appearing as the
first term in the expression for the retarded time measured at the pulsar, it
also appears in the expression for the retarded time measured at the Earth
(for the Earth, just set $L=0$ in equation 4). Therefore, when the rest of our
timing residual model is worked out (see for example the IIB model in
equations 17 and 18), we can simply choose the fiducial time
$t_{0,\mathrm{obs}}=-\frac{D_{c}}{c}$ and all dependence on $D_{c}$ vanishes
(this point is discussed in MC21). Second, is that in principle we could also
choose to include the cosmological curvature constant $k$ as a parameter to
attempt to directly measure it from the gravitational wave signal. Equation 40
provides the connection between $D_{\mathrm{par}}$ and $D_{L}$, and therefore
would introduce both $z$ and $k$ as additional model parameters. However,
since the Fresnel corrections are already smaller order corrections and $k$ is
currently understood to be very close to zero, we will restrict our attention
in this work to a geometrically flat universe and assume $k=0$.

Therefore, working under the assumption of a geometrically flat universe, the
relationships given in Appendix A provide us the connection to the Hubble
constant:

$k=0\quad\longrightarrow\quad\begin{cases}D_{L}=(1+z)D_{\mathrm{par}},\\\\[4.0pt]
H_{0}=\frac{c}{R}\int^{z}_{0}\frac{dz^{\prime}}{E(z^{\prime})}=\frac{c}{D_{\mathrm{par}}}\int^{D_{L}/D_{\mathrm{par}}-1}_{0}\frac{\mathrm{d}z^{\prime}}{E(z^{\prime})}.\end{cases}$
(7)

By procuring a measurement of both the distances $D_{L}$ and
$D_{\mathrm{par}}$ from our pulsar timing model (see Section 2.3.1), we can
directly measure the source’s redshift. By then assuming values of the
cosmological density parameters which appear in the Hubble function $E$ (see
equation 32), or by using the small redshift approximation (see the discussion
around equation 20), we can directly measure the Hubble constant $H_{0}$. In
this work we simulate a flat $\Lambda$CDM universe using the density
parameters values $\Omega_{r}=0$, $\Omega_{M}=0.3$, $\Omega_{\Lambda}=0.7$,
and $\Omega_{k}=0$, and a Hubble constant value of $H_{0}=70\ \mathrm{km}\
\mathrm{s}^{-1}\ \mathrm{Mpc}^{-1}$.

###  2.3 Constructing the Gravitational Wave Timing Residual

We now derive the response of a PTA to the passing of a gravitational wave in
the most general regime (IIB of Figure 1), in an expanding universe.
Throughout, the notation and conventions used in this paper mirror those in
MC21. The main pedantic distinction is the “obs” subscript on many of the
quantities, simply to remind the reader that unlike in MC21, here there is a
difference between source frame quantities and “observer” frame quantities.

The metric perturbation produced by a circular binary system can be written as
(Maggiore, 2008; Creighton & Anderson, 2011):

$\begin{cases}h_{+}(t_{\mathrm{obs}})&\equiv-h(t_{\mathrm{obs}})\cos\big{(}2\Theta(t_{\mathrm{obs}})\big{)},\\\
h_{\times}(t_{\mathrm{obs}})&\equiv-h(t_{\mathrm{obs}})\sin\big{(}2\Theta(t_{\mathrm{obs}})\big{)},\\\
h(t_{\mathrm{obs}})&\equiv\frac{4(G\mathcal{M}_{\mathrm{obs}})^{5/3}}{c^{4}D_{L}}\omega(t_{\mathrm{obs}})^{2/3},\\\
h_{0,\mathrm{obs}}&\equiv\frac{4(G\mathcal{M}_{\mathrm{obs}})^{5/3}}{c^{4}D_{L}}\omega_{0,\mathrm{obs}}^{2/3},\end{cases}$
(8)

for ‘plus’ and ‘cross’ polarizations, and where $t_{\mathrm{obs}}$ denotes the
observer’s clock time. The angular phase and frequency functions
$\Theta(t_{\mathrm{obs}})$ and $\omega(t_{\mathrm{obs}})$ are defined by the
monochromatic regime (A) or the frequency evolving regime (B), depending on
which model we desire:

$\displaystyle\underset{(A)}{\text{Monochromatism}}\hskip 14.22636pt$
$\displaystyle\begin{cases}\begin{tabular}[]{l
l}$\omega(t_{\mathrm{obs}})$&$=\omega_{0,\mathrm{obs}}$ ,\\\
$\Theta(t_{\mathrm{obs}})$&$=\theta_{0}+\omega_{0,\mathrm{obs}}(t_{\mathrm{obs}}-t_{0,\mathrm{obs}})$
,\end{tabular}\end{cases}$ (9) $\displaystyle\underset{(B)}{\text{Frequency
Evolution}}\quad$ $\displaystyle\begin{cases}\begin{tabular}[]{l l l
r}$\omega(t_{\mathrm{obs}})$&$=\omega_{0,\mathrm{obs}}\left[1-\frac{t_{\mathrm{obs}}-t_{0,\mathrm{obs}}}{\Delta\tau_{c,\mathrm{obs}}}\right]^{-3/8}$
,&&\\\
$\Theta(t_{\mathrm{obs}})$&$=\theta_{0}+\theta_{c}\left[1-\left(\frac{\omega(t_{\mathrm{obs}})}{\omega_{0,\mathrm{obs}}}\right)^{-5/3}\right]$
,&&\\\
$\Delta\tau_{c,\mathrm{obs}}$&$\equiv\frac{5}{256}\left(\frac{c^{3}}{G\mathcal{M}_{\mathrm{obs}}}\right)^{5/3}\frac{1}{\omega_{0,\mathrm{obs}}^{8/3}}$&$=(1+z)\Delta\tau_{c}$
,&\\\
$\theta_{c}$&$=\frac{1}{32}\left(\frac{c^{3}}{G\mathcal{M}_{\mathrm{obs}}\omega_{0,\mathrm{obs}}}\right)^{5/3}$&$\equiv\frac{8}{5}\Delta\tau_{c,\mathrm{obs}}\omega_{0,\mathrm{obs}}$&$=\frac{8}{5}\Delta\tau_{c}\omega_{0}$
,\end{tabular}\end{cases}$ (10)

where $t_{0,\mathrm{obs}}$ here (and below) denotes the fiducial time for the
model. The frequency evolution regime is governed by the physically
significant quantities $\Delta\tau_{c,\mathrm{obs}}$ which is the “observed
coalescence time,” and $\theta_{c}$ which is the “coalescence angle” (the
total angle swept out in the source’s orbital plane before the system
coalesces). A more in-depth discussion of these quantities and the physical
assumptions in these regimes is given in MC21.

The strength of the gravitational wave as it reaches the path that the pulsar
photons take between the pulsar and the Earth also depends on the relative
Earth-pulsar-source geometrical alignment. With the Earth at the center of the
coordinate system, the source frame orientation vectors and the Earth-to-
pulsar vector are defined by:

$\begin{cases}\begin{tabular}[]{c l c r
l}$\hat{r}$&$\equiv\big{[}\sin(\theta)\cos(\phi),$&$\sin(\theta)\sin(\phi),$&$\cos(\theta)\big{]}$
\ ,&(Earth-to-gravitational wave source unit vector)\\\
$\hat{\theta}$&$\equiv\big{[}\cos(\theta)\cos(\phi),$&$\cos(\theta)\sin(\phi),$&$\sin(\theta)\big{]}$
\ ,&(transverse plane basis vector)\\\
$\hat{\phi}$&$\equiv\big{[}-\sin(\phi),$&$\cos(\phi),$&$0\big{]}$ \
,&(transverse plane basis vector)\\\
$\hat{p}$&$\equiv\big{[}\sin(\theta_{p})\cos(\phi_{p}),$&$\sin(\theta_{p})\sin(\phi_{p}),$&$\cos(\theta_{p})\big{]}$
\ ,&(Earth-to-pulsar unit vector)\end{tabular}\end{cases}$ (11)

from which we define the following generalized polarization tensors:

$\displaystyle\begin{cases}e^{\hat{r}+}_{ij}&\equiv\hat{\theta}_{i}\hat{\theta}_{j}-\hat{\phi}_{i}\hat{\phi}_{j},\\\
e^{\hat{r}\times}_{ij}&\equiv\hat{\phi}_{i}\hat{\theta}_{j}+\hat{\theta}_{i}\hat{\phi}_{j},\\\
E^{\hat{r}+}_{ij}&\equiv\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)e^{\hat{r}+}_{ij}+\sin(2\psi)e^{\hat{r}\times}_{ij}\right],\\\
E^{\hat{r}\times}_{ij}&\equiv\cos(\iota)\left[-\sin(2\psi)e^{\hat{r}+}_{ij}+\cos(2\psi)e^{\hat{r}\times}_{ij}\right].\end{cases}$
(12)

This notation groups all of the geometrical orientation and location angles
into the definition of the polarization tensor, in order to keep them from
mixing the plus and cross metric perturbations. With these we then define the
following antenna patterns which decide the detector’s sensitivity to the
gravitational wave based on the Earth-pulsar-source geometrical alignment:

$\displaystyle\begin{cases}\begin{tabular}[]{l l
l}$f^{+}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}e^{\hat{r}+}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\frac{\hat{p}^{i}\hat{p}^{j}\hat{\theta}_{i}\hat{\theta}_{j}-\hat{p}^{i}\hat{p}^{j}\hat{\phi}_{i}\hat{\phi}_{j}}{\left(1-\hat{r}\cdot\hat{p}\right)}\quad=\frac{\left(\hat{p}\cdot\hat{\theta}\right)^{2}-\left(\hat{p}\cdot\hat{\phi}\right)^{2}}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\\\\[3.0pt]
$f^{\times}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}e^{\hat{r}\times}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\frac{\hat{p}^{i}\hat{p}^{j}\hat{\phi}_{i}\hat{\theta}_{j}+\hat{p}^{i}\hat{p}^{j}\hat{\theta}_{i}\hat{\phi}_{j}}{\left(1-\hat{r}\cdot\hat{p}\right)}\quad=\frac{2\left(\hat{p}\cdot\hat{\theta}\right)\left(\hat{p}\cdot\hat{\phi}\right)}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\\\\[3.0pt]
$F^{+}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}E^{\hat{r}+}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)f^{+}+\sin(2\psi)f^{\times}\right]$
,\\\\[3.0pt]
$F^{\times}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\times}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\cos(\iota)\left[-\sin(2\psi)f^{+}+\cos(2\psi)f^{\times}\right]$
.\end{tabular}\end{cases}$ (13)

The final gravitational waveform in the transverse-traceless gauge along the
$\hat{r}$-axis can now be expressed as:

$h^{TT\hat{r}}_{ij}=E^{\hat{r}+}_{ij}h_{+}+E^{\hat{r}\times}_{ij}h_{\times}=E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times].$
(14)

In the local observer frame the gravitational waves will affect the timing of
local pulsars. The derivation at this point remains unchanged to what was
shown in MC21 (see, for example, Section A2 of that paper), and we again
compute the effect that the gravitational waves have on the timed signals of
these pulsars. The gravitational wave-induced fractional shift of the pulsar’s
period $T$ is:

$\frac{\Delta
T^{\mathrm{GW}}}{T}(t_{\mathrm{obs}})\approx\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{obs}}}_{t_{\mathrm{obs}}-\frac{L}{c}}\frac{\partial
h_{\textsc{A}}\Big{(}t_{\mathrm{obs,ret}}\left(t^{\prime}_{\mathrm{obs}},\vec{x}\right)\Big{)}}{\partial
t^{\prime}_{\mathrm{obs}}}\Bigg{\rvert}_{\vec{x}=\vec{x}_{0}\left(t^{\prime}_{\mathrm{obs}}\right)}\mathrm{d}t^{\prime}_{\mathrm{obs}},$
(15)

where $t_{\mathrm{obs}}$ is the time a pulsar’s photon is observed arriving at
Earth, $t_{\mathrm{obs,ret}}$ is the retarded time of the gravitational wave,
and $\vec{x}_{0}(t_{\mathrm{obs}})$ is the spatial path of the photon between
the pulsar and the Earth. Finally the gravitational wave-induced timing
residual is the integrated fractional period shift due to the gravitational
wave over the observation time. Conceptually this is the difference between
the observed and expected time-of-arrival of a pulsar’s pulse (Creighton &
Anderson, 2011; Schneider, 2015; Maggiore, 2018):

$\mathrm{Res}^{\mathrm{GW}}(t_{\mathrm{obs}})=\int\frac{\Delta
T^{\mathrm{GW}}}{T}\left(t^{\prime}_{\mathrm{obs}}\right)\mathrm{d}t^{\prime}_{\mathrm{obs}}=\int\frac{T_{\mathrm{obs}}\left(t^{\prime}_{\mathrm{obs}}\right)-T}{T}\mathrm{d}t^{\prime}_{\mathrm{obs}}=\
\mathrm{Obs}(t_{\mathrm{obs}})-\mathrm{Exp}(t_{\mathrm{obs}}).$ (16)

####  2.3.1 The Fresnel, Frequency Evolution Model (Regime IIB)

Section 3 of MC21 details the four gravitational wave-induced timing residual
model regimes IA-IIB, characterized by frequency evolution and curvature
effects (see Figure 1). Generalizing to a flat cosmologically expanding
universe does not change the derivation behind those four models, but it does
change the interpretation of some of the model parameters.

Specifically, cosmological redshift causes the luminosity distance $D_{L}$ to
become the parameter that appears in the amplitude of the metric perturbation
(equation 8), while similarly, we now recognize that the chirp mass and
orbital frequency parameters are measured in the observer frame of reference,
and are no longer equivalent to their values in the source frame (although the
source frame values can be obtained if the source redshift $z$ is recovered
through parameter estimation and equation 7). Importantly, the parallax
distance $D_{\mathrm{par}}$ (i.e. the transverse comoving distance $D_{M}$, or
comoving distance, $D_{c}$, in our flat universe; see equation 40) now enters
via the retarded time calculated in a flat expanding universe (equation 4).

Hence, the most general IIB model becomes:

$\displaystyle\mathrm{Res}^{\mathrm{GW}}(t_{\mathrm{obs}})$
$\displaystyle=\frac{F^{\textsc{A}}}{4}\left[\frac{h_{\textsc{A}}\Big{(}\omega_{0E},\Theta_{E}-\frac{\pi}{4}\Big{)}}{\omega_{0E}}-\frac{h_{\textsc{A}}\Big{(}\overline{\omega}_{0P},\overline{\Theta}_{P}-\frac{\pi}{4}\Big{)}}{\overline{\omega}_{0P}}\right]\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(17) $\displaystyle\underset{\left(t_{0,\mathrm{obs}}\ =\
-\frac{D_{c}}{c}\right)}{\mathrm{where}}\quad\begin{cases}\omega_{0E}&\equiv\omega_{0,\mathrm{obs}},\\\
\overline{\omega}_{0P}&\equiv\omega_{0,\mathrm{obs}}\left[1+\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}+\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{D_{\mathrm{par}}}}{\Delta\tau_{c,\mathrm{obs}}}\right]^{-3/8},\\\
\theta_{0E}&=\theta_{0},\\\
\overline{\theta}_{0P}&=\theta_{0}+\theta_{c}\left(1-\left[1+\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}+\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{D_{\mathrm{par}}}}{\Delta\tau_{c,\mathrm{obs}}}\right]^{5/8}\right),\\\
\Theta_{E}&=\theta_{0E}+\omega_{0E}t_{\mathrm{obs}},\\\
\overline{\Theta}_{P}&=\overline{\theta}_{0P}+\overline{\omega}_{0P}t_{\mathrm{obs}},\end{cases}$
(18)

along with equations 8, 10, 11, 12, and 13. Written this way the timing
residual is the sum of an “Earth term” (subscripted with an “E”) and a “pulsar
term” (subscripted with a “P”). Note that the overbar notation here is used
for the same distinguishing purpose as in MC21, to remind the reader that this
model has not been analytically derived, but rather proposed. MC21 explains
why this model is physically motivated.

##  3 Methods

###  3.1 Model Parametrization

In this cosmological framework, the gravitational wave timing residual in the
full IIB regime is governed by the source parameters:
$\vec{s}=\\{D_{L},D_{\mathrm{par}},\theta,\phi,\iota,\psi,\theta_{0},\mathcal{M}_{\mathrm{obs}},\omega_{0,\mathrm{obs}}\\}$.
Following the explanation given in MC21 we swap the $D_{L}$ parameter in our
model with the “observed Earth term timing residual amplitude” parameter,
defined as:

$A_{E,\mathrm{obs}}\equiv\frac{h_{0,\mathrm{obs}}}{4\omega_{0,\mathrm{obs}}}=\frac{\left(G\mathcal{M}_{\mathrm{obs}}\right)^{5/3}}{c^{4}D_{L}\omega_{0,\mathrm{obs}}^{1/3}}\approx(140\text{
ns})\left(\frac{\mathcal{M}_{\mathrm{obs}}}{10^{9}M_{\odot}}\right)^{5/3}\left(\frac{100\text{
Mpc}}{D_{L}}\right)\left(\frac{1\text{
nHz}}{\omega_{0,\mathrm{obs}}}\right)^{1/3}.$ (19)

Because present day capabilities on measuring the distances to pulsars often
cannot constrain them to better than the order of 100 pc (Arzoumanian et al.,
2018; Deller et al., 2019), which is larger than our gravitational wave
wavelengths (see equation 6), we include all of the pulsar distances as free
parameters in our model (Corbin & Cornish, 2010; Lee et al., 2011, ; MC21).
Therefore we can divide the model parameters into source parameters and pulsar
distance parameters $\vec{X}=\left[\vec{s},\vec{L}\right]$.

To measure $H_{0}$, one measures $D_{\mathrm{par}}$, and also
$A_{E,\mathrm{obs}}$, $\mathcal{M}_{\mathrm{obs}}$, $\omega_{0,\mathrm{obs}}$,
which together give $D_{L}$ via equation 19. Then combining $D_{L}$ and
$D_{\mathrm{par}}$ in equation 7 gives $H_{0}$. Note that this does require
knowledge of the density parameters that go into $H_{0}$ (equation 31). In
practice we take a simplified approach by using the small redshift
approximation to replace the parallax distance parameter $D_{\mathrm{par}}$
with $H_{0}$. This does not qualitatively change the calculation, but allows
computational simplicity for this proof-of-principle study and is relevant for
the source distances for which Fresnel effects are prominent (; ). For a flat
$k=0$ universe and for $z\ll 1$, the Hubble Law is
$z\approx\frac{H_{0}}{c}D_{L}$, so we can combine equations 40 and 19 to
write:

$\displaystyle D_{\mathrm{par}}$
$\displaystyle\approx\left[\frac{c^{4}A_{E,\mathrm{obs}}\omega_{0,\mathrm{obs}}^{1/3}}{\left(G\mathcal{M}_{\mathrm{obs}}\right)^{5/3}}+\frac{H_{0}}{c}\right]^{-1}\qquad\qquad\qquad(\text{for
}k=0,\ z\ll 1),$ (20)

Note that if one were to use the equally valid small redshift approximation
$z\approx\frac{H_{0}}{c}D_{\mathrm{par}}$, then our result for $H_{0}$ would
differ by $\mathcal{O}(z^{2})$ ($\lesssim$ a few percent for the fiducial
cases in this study).

This approach is particularly useful when working with the two source problem
described in Section 4.2, as equation 20 will reduce the dimensionality of the
model by one ($D_{\mathrm{par,1}}$ and $D_{\mathrm{par,2}}$ are both replaced
by $H_{0}$), therefore giving us the direct joint posterior on $H_{0}$ from
the two sources. In the full calculation, one must make an estimate of the
joint posterior on $H_{0}$ (for example, with a kernel density estimate) using
the measured parameters from both sources. The low-redshift approach builds
into our model the additional knowledge that $H_{0}$ is a constant
irrespective of the source we are detecting. But again the approximation is
only valid for $z\ll 1$.

As a final note, unless otherwise stated, for the studies presented in this
work we chose to inject the following default values for the gravitational
wave source parameters: $D_{\mathrm{par}}=100$ Mpc, $D_{L}=102.35$ Mpc,
$\mathcal{M}_{\mathrm{obs}}=10^{9}\ \mathrm{M}_{\odot}$,
$\omega_{0,\mathrm{obs}}=30$ nHz, and angular parameters $\\{\theta,\ \phi,\
\iota,\ \psi,\ \theta_{0}\\}=\\{\frac{3\pi}{7},\ \frac{5\pi}{3},\
\frac{2\pi}{5},\ \frac{\pi}{7},\ 1\\}$.

###  3.2 Likelihood and Priors

Following MC21, for this work we assume that the timing residual data is the
sum of the underlying gravitational wave-induced residual plus some random
noise, that is
$\overrightarrow{\mathrm{Res}}=\overrightarrow{\mathrm{Res}}^{\mathrm{GW}}+\vec{N}$.
For simplicity the noise we model is white noise, and each data point
collected has some uncertainty uncorrelated between observations/pulsars, that
is the timing covariance matrix
$\mn@boldsymbol{\Sigma}=\mathrm{diag}\left(\sigma_{1}^{2},\sigma_{2}^{2},\ldots,\sigma_{d}^{2}\right)$.
Therefore the likelihood function we propose and use in this work is:

$\displaystyle\mathcal{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{X}\right)$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{d}\mathrm{det}(\mn@boldsymbol{\Sigma})}}\exp{\left[-\frac{1}{2}\bigg{(}\overrightarrow{\mathrm{Res}}-\overrightarrow{\mathrm{Res}}^{\mathrm{GW}}\left(\vec{X}\right)\bigg{)}^{T}\mn@boldsymbol{\Sigma}^{-1}\bigg{(}\overrightarrow{\mathrm{Res}}-\overrightarrow{\mathrm{Res}}^{\mathrm{GW}}\left(\vec{X}\right)\bigg{)}\right]},$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{d}}\sigma^{d}}\prod^{d}_{a=1}\exp{\left[-\frac{1}{2}\frac{1}{\sigma_{a}^{2}}\left(\mathrm{Res}_{a}-\mathrm{Res}^{\mathrm{GW}}_{a}\left(\vec{X}\right)\right)^{2}\right]}.$
(21)

Here $\overrightarrow{\mathrm{Res}}$ is the $d$-dimensional timing residual
data (with dimension equal to the number of pulsars times the number of
observations per pulsar, indexed by $a$),
$\overrightarrow{\mathrm{Res}}^{\mathrm{GW}}$ is our residual model (for the
IIB model, see Section 2.3.1) for every data point, and the model parameters
are contained in the vector $\vec{X}$. Log-parameters are used for the non-
angular parameters
$\left\\{\log_{10}\left(A_{E,\mathrm{obs}}/A_{E,\mathrm{obs},\mathrm{true}}\right),\log_{10}\left(D_{\mathrm{par}}/D_{\mathrm{par,true}}\right),\log_{10}\left(\mathcal{M}_{\mathrm{obs}}/\mathcal{M}_{\mathrm{obs,\mathrm{true}}}\right),\right.$
$\left.\log_{10}\left(\omega_{0,\mathrm{obs}}/\omega_{0,\mathrm{obs},\mathrm{true}}\right)\right\\}$,
as well as $\log_{10}\left(H_{0}/H_{0,\mathrm{true}}\right)$. For our priors,
we required all of the physical parameters be non-negative
($A_{E,\mathrm{obs}}$, $\mathcal{M}_{\mathrm{obs}}$,
$\omega_{0,\mathrm{obs}}$, $D_{L}$, $D_{\mathrm{par}}$, $L$, $H_{0}$), and
that $D_{\mathrm{par}}\leq D_{L}$ (which assumes $z>0$ in equation 40). For
the angular parameters we placed the general boundaries: $0<\theta,\ \iota,\
\theta_{0}<\pi$, and $0<\phi,\ \psi<2\pi$.

###  3.3 Mock Observations

In order to create the mock timing residual data for this study, we generate
different mock PTAs, namely “fiducial PTAs” and “Square Kilometre Array
(SKA)-era PTAs.” The details for constructing these arrays follow below, but
each serves a distinct purpose in this study. While the fiducial PTA is
physically motivated, it is a simpler construction, and is primarily used
throughout this work to look for various trends and scaling laws in our
predicted measurement capabilities. On the other hand, the various SKA PTAs
are meant to represent more realistic hypothetical future timing arrays, and
are therefore used to simulate more realistic results.

Note that for simplicity, in this work we do not add noise $\vec{N}$ on top of
our gravitational wave-induced timing residuals. Therefore the simulated
timing residuals are purely the gravitational wave timing residual component,
which is why in the MCMC results shown (Figures 2, 4, and 11) the posteriors
are centered on the true injected parameters. Finally, for all of our
simulations in this paper, we use an observation time of 10 years, and a
timing cadence of 30 observations per year.

####  3.3.1 A Fiducial PTA

To gauge dependence of parameter recovery on pulsar number and timing
precision we choose a simple fiducial PTA. The fiducial PTA is realistically
motivated by using a simple Milky Way Galaxy structure model to create a
density distribution of pulsars $n_{\mathrm{psr}}$ (Schneider, 2015) to sample
from,

$n_{\mathrm{psr}}\propto\exp\left(-\frac{r}{h_{R}}\right)\Bigg{[}\exp\left(-\frac{|z|}{h_{\mathrm{thin}}}\right)+0.02\exp\left(-\frac{|z|}{h_{\mathrm{thick}}}\right)\Bigg{]}\times\exp\left(-\frac{1}{2}\left[\left(\frac{x}{R}\right)^{2}+\left(\frac{y}{R}\right)^{2}+\left(\frac{z}{R}\right)^{2}\right]\right)H(r-r_{\mathrm{min}}),$
(22)

where the center of the coordinate system is at the Earth’s position, the
Galactic center is at $x_{\mathrm{MW}}\approx 8$ kpc, and
$r\equiv\sqrt{(x-x_{\mathrm{MW}})^{2}+y^{2}}$. Here $h_{R}$ is the scale
length of the Galactic disc ($h_{R}\approx 3.5$ kpc), $h_{\mathrm{thin}}$ is
the scale height of the thin disc ($h_{\mathrm{thin}}\approx 0.325$ kpc), and
$h_{\mathrm{thick}}$ is the scale height of the thick disc
($h_{\mathrm{thick}}\sim 1.5$ kpc). The first term in this expression creates
the distribution of stars in the Milky Way, while the second term simply
places preference on stars within a Gaussian ball centered on Earth. The
motivation for including this second term is simply the idea that we will
likely be more sensitive to timing pulsars within some scale distance $R$ from
the Earth. We also placed a minimum distance $r_{\mathrm{min}}$ on this sphere
(hence the Heaviside function). For our fiducial PTA we chose $R=5$ kpc and
$r_{\mathrm{min}}=0.1$ kpc. An example of 1000 pulsars generated from this
distribution is shown in Figure 9.

####  3.3.2 Simulated SKA-era PTAs

In addition to the fiducial PTA, we use the PSRPOPPY package (Bates et al.,
2014) to generate a population of millisecond pulsars (MSPs) that could be
discovered and used for timing in the SKA era. We begin by generating a
population of $3\times 10^{4}$ MSPs in the galaxy (motivated by Smits et al.,
2009) and simulate an SKA-like survey which detects $8893$ of these with
signal-to-noise (SNR) above 9. We separate the detected MSPs into 20 radial
bins out to the furthest detected pulsars at $\sim 18$ kpc (with the closest
at $\sim 0.2$ kpc). To simulate which MSPs will be best for high-precision
timing, we sort the MSPs by SNR in each radial bin and retain the top $p\%$
for the final PTA. In the outer three bins which contain $<10$ MSPs we take
only the highest SNR pulsar, each of which have SNRs $\geq 20$. Depending on
the value of $p$, the final PTAs contain between $\sim 40$ and 1800 pulsars
(See Appendix D for further details).

Because our SKA simulation allows us to estimate detection SNRs and spin
periods for the simulated SKA-era population of MSPs, we can additionally
estimate a timing uncertainty for the $a^{\mathrm{th}}$ MSP in the array. To
do so, we assume that the pulse time-of-arrival uncertainty is proportional to
the radiometer noise (e.g. Verbiest & Shaifullah, 2018),

$\sigma_{a}\propto\frac{P_{a}}{\sqrt{\Delta T_{a}}\mathrm{SNR}_{a}},$ (23)

for integration time $\Delta T$, pulse period $P$, and detection signal-to-
noise SNR. While this represents only part of the noise budget, it allows a
study where pulsars have a more realistic spread of timing uncertainties than
in the fiducial PTAs. We then consider two scenarios, one where the
integration time is constant for all pulsars ($\Delta T_{a}=\mathrm{cst.}$),
and one where the observation strategy is intelligently chosen to boost the
signal from lower SNR MSPs ($\Delta T_{a}\propto\mathrm{SNR}^{-1}_{a}$),

$\sigma_{a}=\sigma_{\min}\left(\frac{P_{a}}{1\mathrm{ms}}\right)\begin{cases}\left(\frac{\mathrm{SNR}_{a}}{\mathrm{SNR}_{\max}}\right)^{-1},&\Delta
T_{a}=\mathrm{cst.}\\\
\left(\frac{\mathrm{SNR}_{a}}{\mathrm{SNR}_{\max}}\right)^{-1/2},&\Delta
T_{a}\propto\mathrm{SNR}^{-1}_{a},\end{cases}$ (24)

where we choose a characteristic, best-case timing uncertainty
$\sigma_{\mathrm{min}}=\\{1,10\\}$ ns and enforce
$\sigma_{a}\geq\sigma_{\mathrm{min}}$. While somewhat arbitrary, we assume an
$\mathrm{SNR}_{\max}=10^{4}$, corresponding to an SNR threshold above which
pulsars with $1$ ms pulse periods can be timed to $\sigma_{\mathrm{min}}$.
Mock array pulsar distributions and timing uncertainties are shown in Appendix
D.

###  3.4 Fisher Matrix and MCMC

To recover best-fit source and pulsar parameter values we carry out both
Fisher matrix and Markov Chain Monte Carlo (MCMC) analyses which compare our
model and mock timing residuals through the likelihood (equation 21).

A Fisher matrix is a useful tool for quick parameter estimation. Calculating
the inverse of the Fisher matrix gives the estimated parameter covariance
matrix $\textbf{{F}}=\textbf{{C}}^{-1}$ for the experiment. Therefore finding
a model’s inverse Fisher matrix can help to tell us which parameters are
covariant with each other, and roughly how well we might expect to recover
each model parameter given our experimental set-up. Computationally, this is
efficient and can quickly allow us to test many different experiments. Fisher-
based surveys are useful for quickly searching large parts of parameter space
for trying to understand what types of sources would result in good parameter
estimation. But it does require that one assume a more robust search (such as
with an MCMC analysis) could successfully identify the true mode from any
potential secondary modes. This limitation is due to the fact that the Fisher
matrix is only meant to approximate the shape of the posterior near the
maximum likelihood. If for example the posterior is multimodal (see Section
4), then the Fisher matrix will not capture this behavior. For a more
extensive (and computationally expensive) targeted study of the posterior, we
use MCMC.

From our likelihood equation 21 and the definition of the Fisher matrix as
$\textsf{F}_{ij}\equiv\Big{\langle}\left(\frac{\partial}{\partial
X_{i}}\ln{\mathcal{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{X}_{\mathrm{true}}\right)}\right)\left(\frac{\partial}{\partial
X_{j}}\ln{\mathcal{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{X}_{\mathrm{true}}\right)}\right)\Big{\rangle}$
we have:

$\textsf{F}_{ij}=\sum_{a}\frac{1}{\sigma_{a}^{2}}\left(\frac{\partial\mathrm{Res}^{\mathrm{GW}}_{a}\left(\vec{X}_{\mathrm{true}}\right)}{\partial
X_{i}}\right)\left(\frac{\partial\mathrm{Res}^{\mathrm{GW}}_{a}\left(\vec{X}_{\mathrm{true}}\right)}{\partial
X_{j}}\right).$ (25)

In some instances, namely when using the Fresnel models, inverting the Fisher
matrix cannot be done accurately because inclusion of the pulsar distance
parameters introduces too much uncertainty and covariance amongst the other
parameters (this typically happens with a single gravitational wave source).
In these instances we add uncorrelated Gaussian pulsar distance priors to the
Fisher matrix (Wittman, unpublished notes111Wittman D., no date, Fisher Matrix
for Beginners, UC Davis, http://wittman.physics.ucdavis.edu/Fisher-matrix-
guide.pdf.; Iacovelli et al., 2022). See Section 5 of MC21 for more details on
adding these priors, and for a discussion of why this problem is more
prominent in the Fresnel regimes.

Since all parameters in a Fisher matrix approximation are normally
distributed, we often use the coefficient of variation (CV), the predicted
fractional error on a given parameter, as a useful way of quantifying the
measurability of a given parameter when performing a Fisher matrix analysis:

$\mathrm{CV}_{x}\equiv\frac{\text{Standard Deviation of }x}{\text{Expectation
Value of }x}=\begin{cases}\frac{\sigma}{\mu},&\quad\text{Normal Distribution
in $x$}\\\\[8.0pt]
\sqrt{e^{\sigma^{2}\ln(10)^{2}}-1}.&\quad\text{Log}_{10}\text{-Normal
Distribution in $x$}\end{cases}$ (26)

Here $\mu$ and $\sigma$ are the distributions’ parameters. If the parameter
$x$ is normally distributed then $\mu$ and $\sigma$ are also the
distribution’s mean and standard deviation, and if the parameter is log-
normally distributed then CV only depends on the log-normal $\sigma$
parameter. A resulting CV of the order of unity or larger suggests the
parameter would not be measureable given the source and experiment, while
smaller CVs suggest better parameter recovery.

We use the PYTHON emcee package (Foreman-Mackey et al., 2013) to carry out
MCMC sampling of the posterior. We use differential evolution jump proposals
with the number of walkers set between 5 and 7 times the number of parameter
dimensions, and the results shown in this paper ran between $350~{}000$ to
$500~{}000$ iterations. We target our searches by initializing the MCMC
walkers in very small Gaussian “balls” about the true injected parameters, as
well as around nearby theoretical secondary modes of interest. For the
covariance matrix of these Gaussian balls we use the inverse of the Fisher
matrix prediction, typically multiplied by some overall scale factor. This
choice helps the MCMC simulations more quickly find and explore the local
distributions around the true parameters.

While our MCMC simulations help provide proof-of-principle for the main ideas
in this paper, future work should explore using different MCMC samplers to
perform consistency checks on the work presented here, such as with parallel
tempering or nested samplers (Samajdar et al., 2022). Further investigation of
custom jump proposals should also be implemented to try and improve the
exploration of walkers from one posterior mode to another.

##  4 A Two-Source Solution to the Pulsar Distance Wrapping Problem

The pulsar distance wrapping problem creates a challenge for parameter
estimation in this work. As explained in detail in MC21, due to the way the
pulsar distance $L$ affects the phase of the pulsar term (equation 18),
increasing or decreasing the pulsar distance by specific $\Delta L$ values
will cause the phase to wrap around the interval $\left[0,2\pi\right)$,
resulting in the same timing residual for multiple pulsar distances.

However, we show here that the wrapping problem is prominent primarily when
there is only one gravitational wave source in the signal. When there are two
or more sources simultaneously creating residuals in every pulsar, the
degeneracy of the wrapping cycle is broken. This has a significant effect on
our ability to recover the parallax distance $D_{\mathrm{par}}$ and
subsequently the Hubble constant, largely because breaking the wrapping
degeneracy can strongly constrain the pulsar distances.

###  4.1 Classifying the Error Envelope (Single Source)

More formally, in the monochromatic regimes IA and IIA, the wrapping cycle
degeneracy arises because the timing residual is identical for $L\pm\Delta
L_{n}$:

$\Delta
L_{n}\equiv\begin{cases}n\frac{\lambda_{\mathrm{gw}}}{\left(1-\hat{r}\cdot\hat{p}\right)},&\quad\text{(Plane-
Wave, monochromatic)}\\\\[10.0pt]
-\left(\frac{D_{\mathrm{par}}}{\left(1+\hat{r}\cdot\hat{p}\right)}+L\right)+\sqrt{\left(\frac{D_{\mathrm{par}}}{\left(1+\hat{r}\cdot\hat{p}\right)}+L\right)^{2}+n\frac{2\lambda_{\mathrm{gw}}D_{\mathrm{par}}}{\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)}},&\quad\text{({Heuristic}
Fresnel, monochromatic)}\end{cases}$ (27)

for $n\in\mathbb{Z}$.222See MC21 and note that in equation 27 we replace
$R\rightarrow D_{\mathrm{par}}$ to fit with our generalization described in
Section 2. Hence, the wrapping cycle distance $\Delta L$ is a scaled version
of the gravitational wave wavelength $\lambda_{\mathrm{gw}}$, which includes
geometric factors dependent on the relative positions of the pulsar and the
source. Most notably, the more aligned a source and pulsar are on the sky
(within 90°, $0<\hat{r}\cdot\hat{p}<1$), the larger the wrapping cycle
becomes, and vice versa as the source becomes more anti-aligned (more than
90°, $-1<\hat{r}\cdot\hat{p}<0$).

Consider the problem now in terms of parameter estimation using the likelihood
function equation 21. If we work entirely in one of the monochromatic regimes
IA or IIA, meaning we use one of these models to both calculate the timing
residuals of our injected source $\overrightarrow{\mathrm{Res}}$ and recover
the parameters $\overrightarrow{\mathrm{Res}}\left(\vec{X}\right)$, then the
likelihood function is perfectly multimodal at every $L\pm\Delta L_{n}$ for
every pulsar in the PTA (see the left-hand panel of Figure 2). Therefore in
either of these regimes, we cannot identify the true pulsar distance with this
likelihood since all modes have equal probability.

1 Source IA Model
(Plane-Wave, Monochromatic)

1 Source IB Model
(Plane-Wave, Frequency Evolution)

Figure 2:  Corner plots demonstrating the pulsar distance error envelope from
an MCMC simulation of a PTA with a single gravitational wave source (default
parameters). The results are compared using the (non-Fresnel regime) IA versus
IB models (MC21) for the exact same PTA and source. The full PTA contained 40
pulsars from the fiducial PTA (see Figure 11), with a uniform timing
uncertainty of $\sigma=0.2$ ns. This corner plot singles out three pulsars and
displays the modal overlap criterion ratios from equation 28 above each 1D
histogram. The bottom row shows the (un-normalized) log-posterior values.
Left-hand panel: For the monochromatic IA model, each mode receives the same
support in the likelihood function at every wrapping cycle. Right-hand panel:
For the IB model, we see how frequency evolution helps us to weaken the
wrapping cycle degeneracy by reducing support in the likelihood at secondary
wrapping cycles modes. Note the samples in this panel have been cleaned, such
that all samples with log-posterior $\leq-100$ have been removed from the data
set (see the discussion of cleaning in Section 4.3).

In the frequency evolution regimes IB and IIB this wrapping cycle does not
formally exist as the frequencies are now time-dependent. While this
technically breaks the degeneracy, there still exists some support at
$L\pm\Delta L_{n}$. That is, the IB and IIB likelihood functions are still
multimodal at each pulsar’s true distance modulo its wrapping cycle distance
(see the right-hand panel of Figure 2). However, as moving out a wrapping
cycle no longer returns the same exact timing residual as the timing residual
at the true pulsar distance, these modes become less and less probable for
every additional cycle away from the true distance.

Additionally, every mode will have some width. If these modes are close enough
to each other (relative to the size of their widths), then the modes will
blend together. Corbin & Cornish identified this blending of the pulsar
distance wrapping modes in their model as an “error envelope” which forms
about the true pulsar distance. The effect is that the true distance will be
further buried within the uncertainty surrounding the secondary modes, making
it much harder to measure in practice (see for example the top-most 1D
histogram in both panels of Figure 2). However, Corbin & Cornish did not fully
classify the properties of the error envelope. Namely the geometric factors in
the wrapping cycle (equation 27) are of crucial importance in determining the
separation between the true mode and the secondary modes.

Figure 3: Here we visualize the modal overlap criterion used to analyse the
pulsar wrapping cycle degeneracy (equations 28-30). Left-hand panel: The
“modal overlap” criterion between the primary and secondary modes in the
pulsar distance wrapping problem is defined as the ratio of the separation of
the primary and secondary modes ($\Delta L_{1}$ for a single source, and
$\Delta L^{\mathrm{common}}_{1}=|\Delta L_{n}^{\mathrm{S1}}+\Delta
L_{m}^{\mathrm{S2}}|/2$ for two sources) to the width of the modes
($\sigma_{L,\mathrm{Fisher}})$. We classify values of this ratio $<4$ as
“overlapping”, between $4$ and $8$ as “intermediate”, and $>8$ as “resolved.”
Right-hand panel: Individual likelihoods for two sources S1 and S2 will agree
on the location of the primary mode in the pulsar distance, but as long as the
sources are located on different parts of the sky and/or have different
frequencies, they will have different sets of pulsar distance wrapping cycles
(equation 27). Therefore the joint likelihood for both sources will find
support at the distances where there is some common modal overlap between
$\Delta L_{n}^{\mathrm{S1}}$ and $\Delta L_{m}^{\mathrm{S2}}$.

Therefore we construct a useful way of classifying the amount of “modal
overlap” between the primary (thick blue line) and secondary modes (dashed
orange lines). Overlap between the first and second mode will depend on two
factors: the width of each mode, and the separation between modes. If each
mode is a Gaussian of approximately the same width, then the ratio of mode
separation to mode width provides a numerical quantification of the amount of
overlap (see the left-hand panel of Figure 3). We define such a modal overlap
criterion by calculating the ratio of the first wrapping cycle distance
$\Delta L_{1}$ (i.e. the distance between our primary and secondary modes) to
the approximated Fisher width $\sigma_{L,\mathrm{Fisher}}$ of each mode:

$\frac{\Delta L_{1}}{\sigma_{L,\mathrm{Fisher}}}\sim\begin{cases}<4\
,&\text{Overlapping}\\\ 4-8\ ,&\text{Intermediate}\\\ >8\
.&\text{Resolved}\end{cases}$ (28)

Here $\sigma_{L,\mathrm{Fisher}}$ is the uncertainty of a given pulsar
distance as predicted by the Fisher matrix analysis (that is, the standard
deviation of a pulsar distance $L$ calculated from the inverse Fisher matrix,
$\textbf{{C}}=\textbf{{F}}^{-1}$). Conceptually, the values we choose here are
motivated such that a secondary mode “overlapping” with the primary mode will
have it’s peak within $4\sigma_{L,\mathrm{Fisher}}$ of the true value, and
“resolved” when it’s peak is greater than $8\sigma_{L,\mathrm{Fisher}}$ from
the true value. Figure 2 shows an example of three pulsars from a simulation
with overlapping, intermediate, and resolved modes clearly visible. A complete
version is given in Figure 11.

Note that in general the Fisher uncertainty predictions are only valid around
the true mode, and therefore cannot accurately forecast the multimodal
behavior of the pulsar distances. However, in practice we observe that the
secondary and primary modes have approximately the same width. Therefore, we
use $\sigma_{L,\mathrm{Fisher}}$ as our proxy in defining equation 28. This
effectively assumes that the behavior of the likelihood function near each
wrapping cycle will look the same as the behavior of the likelihood function
near the true mode.

In practice, this specific modal overlap criterion also works best if no prior
knowledge is placed on the pulsar distances. This is partly because
$\sigma_{L,\mathrm{Fisher}}$ would be affected by a prior (see discussion in
Section 3.4), while $\Delta L_{1}$ would not. In our MCMC simulations, we find
that when a very constricting prior is placed on the pulsar distances (such as
a prior of the order of the wrapping cycle distance
$\sim\mathcal{O}\left(\Delta L_{1}\right)$), then the location of the
secondary modes in the posterior can shift away from $\Delta L_{1}$ and closer
to the true mode. This results in a change in the numerical ranges for which
we find the overlapping, intermediate, and resolved modes in equation 28.

One other important caveat to note is that this criterion also works best if
there are no secondary modes in the other source parameters. Consider equation
27 for the wrapping cycle modes. In the case of the plane-wave, monochromatic
wrapping cycle distance, $\Delta L_{n}=\Delta
L_{n}\left(\theta,\phi,\omega_{0,\mathrm{obs}}\right)$. If the likelihood has
secondary modes in the frequency parameter $\omega_{0,\mathrm{obs}}$, then
those secondary mode solutions can produce their own wrapping cycle modes,
which we observed in several MCMC tests. Specifically, for some high mass
systems $\mathcal{M}_{\mathrm{obs}}>10^{9}\ \mathrm{M}_{\odot}$, secondary
modes (of lower probability) form in the $\mathcal{M}_{\mathrm{obs}}$ and
$\omega_{0,\mathrm{obs}}$ parameters, which then change the locations of the
observed secondary modes in pulsar distances (meaning our criterion becomes
less accurate in this case). In such cases, these new modes may be the result
of stronger frequency evolution effects within the model, and should be kept
in mind when performing these types of analyses.

Hence, equation 28 provides a new and more general criterion for quantifying
the wrapping problem. While this criterion relies on the Fisher predicted
uncertainty $\sigma_{L,\mathrm{Fisher}}$, which does not offer an immediately
intuitive connection to the system parameter values, equation 28 can still be
calculated for all pulsars in a PTA without requiring a full MCMC simulation.
Therefore we can quickly predict for any pulsar in our PTA whether or not that
pulsar will have an error envelope. The more pulsars in the PTA which are
predicted to have resolved modes, the more likely it will be that we can
obtain a measurement of the desired parallax distance $D_{\mathrm{par}}$ and
hence the Hubble constant.

###  4.2 Breaking the Wrapping Cycle Degeneracy With Two Sources

Now consider two gravitational wave sources leaving their combined signal in
the timing residuals of every pulsar. With a single source, every pulsar has a
specific wrapping cycle distance (equation 27) which primarily depends on the
frequency of the source and the angular sky separation between the source and
the pulsar. But with two sources, every pulsar will have two wrapping cycle
distances. The key idea is that these two sets of wrapping cycle distances
will not be the same as long as the two sources differ in either frequency
and/or angular sky position. Therefore, the joint likelihood for both sources
should only contain secondary modes in an individual pulsar’s distance at
common multiples of both the wrapping cycle distances. In principle adding
even more sources would further break the degeneracy, because secondary modes
should then only form in the joint likelihood at common multiples of all
sources. We restrict our attention to two SMBHB sources in this work, and
leave it open to future work to consider more sources.

Consider the right-hand panel of Figure 3. The joint likelihood should find
some support at common distances where the two separate wrapping cycle
distance uncertainty modes overlap. The amount of support at those common
distances will depend on how strong the modal overlap is. And if there is
strong modal overlap at a common distance, we can then check if that common
mode will overlap with the primary mode to create an error envelope (as we did
in equation 28). Therefore we can define two new criteria. The first criterion
predicts the “common mode overlap strength,” and the second criterion
classifies the “common modal overlap” with the primary mode.

For the common mode overlap strength criterion we choose to write:

$\frac{\left|\Delta L_{n}^{\mathrm{S1}}-\Delta
L_{m}^{\mathrm{S2}}\right|}{\sigma_{L,\mathrm{Fisher}}}\sim\begin{cases}0-2\
,&\text{``strong'' support}\\\ 2-4\ ,&\text{``weak'' support}\\\ \end{cases},$
(29)

where $\Delta L_{n}^{\mathrm{S1}}$ and $\Delta L_{m}^{\mathrm{S2}}$ are the
$n$th and $m$th wrapping cycle distances of source 1 and source 2,
respectively. As an example, if this quantity equaled zero then the $n$th and
$m$th modes would perfectly overlap, hence there would be strong support for a
mode here. At a value of two, the point in between the $n$th and $m$th modes
(the average distance) would be $1\sigma_{L,\mathrm{Fisher}}$ from each peak,
and at a value of four, the point in between the $n$th and $m$th modes would
be $2\sigma_{L,\mathrm{Fisher}}$ from each peak. We then expect that the
common wrapping cycle distance is at approximately the average distance
between these individual modes, that is $\Delta
L^{\mathrm{common}}_{q}\equiv\frac{\Delta L_{n}^{\mathrm{S1}}+\Delta
L_{m}^{\mathrm{S2}}}{2}$, for the $q$th common wrapping cycle (between the
$n$th and $m$th wrapping cycles of sources 1 and 2). Note that when
calculating positive or negative common wrapping cycles, both $\Delta
L^{\mathrm{S1}}_{n}$ and $\Delta L^{\mathrm{S2}}_{m}$ should be the same sign.

Note that for the purpose of our studies, we will only be searching for common
modes where we think there is strong support, as we define in the criterion in
equation 29. The reason for this is that “weak” common modes will not have
nearly as much support in the likelihood function. So even if the modal
overlap between a weak common mode and the primary mode is such that it would
produce an error envelope (using equation 30), we expect that this envelope
would not contribute significantly to widening the uncertainty about the true
pulsar distance.

For the common modal overlap criterion we write the following expression in
analogy to equation 28 (see also the left-hand panel of Figure 3):

$\frac{\Delta
L^{\mathrm{common}}_{1}}{\sigma_{L,\mathrm{Fisher}}}\sim\begin{cases}<4\
,&\text{Overlapping}\\\ 4-8\ ,&\text{Intermediate}\\\ >8\
.&\text{Resolved}\end{cases}$ (30)

This condition checks if the uncertainty around the 1st common wrapping cycle
will be close enough to the uncertainty around the true distance as to blend
those uncertainties into a greater envelope. Once again, in order to make this
prediction we are assuming as a proxy that the width of the true mode
predicted by the Fisher matrix for the joint source likelihood also
approximates the secondary modes at the individual source wrapping cycles.

The key point is that with more than one source, it is more likely that common
modes in the pulsar distance will form further away from the true distance
mode. This means that the uncertainties surrounding these common modes will be
less likely to blend with the uncertainty surrounding the true distance to
form an error envelope, hence resulting in better recovery of the pulsar
distances. Furthermore, recall that in the regimes IB and IIB, frequency
evolution alone nominally breaks the wrapping cycle degeneracy (Section 4.1).
Therefore two sources with significant frequency evolution will further
constrain the pulsar distances. All of this coupled with prior knowledge on
our pulsar distances thanks to electromagnetic observations provide three
separate means of localizing the pulsars in our PTA. In application to our
Hubble constant measurement, better pulsar distance measurements mean we can
better recover the source parallax distance $D_{\mathrm{par}}$ from the
Fresnel effects, which is what we show in Section 5.

###  4.3 Targeted MCMC Results for Two Sources

Next we run a targeted MCMC simulation to test our discussion from Section 4
of how two sources break the wrapping cycle degeneracy. For two gravitational
wave sources, our likelihood function (equation 21) now has
$\overrightarrow{\mathrm{Res}}^{\mathrm{GW}}=\overrightarrow{\mathrm{Res}}^{\mathrm{GW}}_{1}+\overrightarrow{\mathrm{Res}}^{\mathrm{GW}}_{2}$
for sources “1” and “2.” Decomposing this joint likelihood function, we can
write it as the product of the likelihood function of just source 1, the
likelihood of just source 2, and their cross terms. Fortunately adding a
second source does not double our parameter space, since both sources will
share the pulsar distance parameters $\vec{L}$. We simply double the number of
source parameters, so
$\vec{X}=\vec{X}_{1}\cup\vec{X}_{2}=\left[\vec{s}_{1},\vec{s}_{2},\vec{L}\right]$.
Furthermore, using the small redshift approximation (equation 20) to replace
the parallax distance $D_{\mathrm{par,1}}$ and $D_{\mathrm{par,2}}$ parameters
with $H_{0}$ reduces the model parameter dimensionality by one and gives us
the joint posterior recovery of $H_{0}$.

In order to test the ideas in Section 4.2, our methodological approach was to
inject balls of walkers in specific areas of parameter space in order to
strategically initialize them (as mentioned in Section 3.4). All walkers were
initialized near the true injected source parameters, but the pulsar distance
parameters were initialized distinctly. One ball of walkers was placed around
the true mode pulsar distance parameter values ($L_{\mathrm{true}}$); two
balls of walkers were placed at the source one wrapping cycles
$L_{\mathrm{true}}+\Delta L^{\mathrm{S1}}_{-1,+1}$; two balls of walkers were
placed at the source two wrapping cycles $L_{\mathrm{true}}+\Delta
L^{\mathrm{S2}}_{-1,+1}$; and finally, two balls of walkers were placed at the
first predicted strong (such that equation 29 was $\leq 2$) common wrapping
cycles $L_{\mathrm{true}}+\Delta L^{\mathrm{common}}_{-1,+1}$.

This was done because we wanted to see if indeed the walkers placed near the
original wrapping cycles of the individual sources would no longer find any
meaningful support in posterior, as they did in Figure 2. The walkers near the
new predicted common wrapping cycles were placed there in order to see if any
support in the posterior would be found at those locations. However, depending
on how far away those common modes are from the true mode we did not
necessarily expect to find meaningful support their either. This is because we
know that when frequency evolution is included in the model, secondary modes
further away from the true mode become less probable (see again, the right-
hand panel of Figure 2).

The result of a targeted MCMC search is shown in Figure 4, which is identical
to the set-up in Figure 2, but now with the inclusion of a second SMBHB source
(and using the full 2 Source IIB model). For this example, our common modal
overlap criterion predicts that all 40 pulsars would be resolved (again given
that for each pulsar we look at the first common mode where the common mode
overlap strength criterion equation 29 is $\leq 2$). Compare this to the
single source IB model results in Figure 11, where only 10 pulsars were
resolved. Therefore the addition of a second source changed this prediction
such that the remaining 30 pulsars should be resolvable. Additionally, the
predicted $\mathrm{CV}_{\mathrm{H0}}=0.257$. Therefore our expectation was
that we would find no meaningful support in the posterior at the secondary
modes as compared to the primary mode, resulting in a measurement of the
Hubble constant. The left-hand panel of Figure 4 shows that the posterior has
sharp peaks at the true parameter values with broad bases surrounding them.

But if we consider the bottom row of the left-hand panel, which shows the (un-
normalized) log-posterior values of these MCMC results, we also see that the
samples span a large range of parameter space with very low probability.
Therefore if we try “cleaning” this data by simply removing all samples with
log-posterior values below a certain threshold (in this case, we cleaned all
values for which the log-posterior $\leq-2\times 10^{6}$), then we obtain the
results shown in the right-hand panel. This reveals that walkers near the true
mode very closely traced the Fisher matrix prediction, walkers away from the
true mode had very low-probability support, and that there were no clearly
defined modal peaks – which starkly contrasts the results found in Figure 2.
Therefore these results support our conclusion about how two sources can help
break the wrapping cycle degeneracy.

2 Source IIB Model
(Fresnel, Frequency Evolution)

“Uncleaned” Samples

“Cleaned” Samples

Figure 4: A targeted MCMC simulation of the recovery of $H_{0}$ for two
sources (using the IIB model). For simplicity, both sources are identical
except for being located on different parts of the sky. Source 1 is at the
default location (and matches the single source in Figure 2), while source 2
has sky angles $\theta_{2}=\frac{\pi}{3}$ rad and $\phi_{2}=\frac{\pi}{5}$
rad. Here 40 pulsars from our fiducial PTA are timed with $\sigma=0.2$ ns
uncertainty (same as in Figure 2). Both panels are of the same simulation,
showing just the recovery of the two sources’ distances ($D_{L}$ and
$D_{\mathrm{par}}$), and the joint recovery of $H_{0}$. The bottom row shows
the (un-normalized) log-posterior values. For each pulsar, if we look at the
first common mode where the common mode overlap strength criterion equation 29
is $\leq 2$, then for this given source and PTA our common modal overlap
criterion predicts that there will be 0 overlapping, 0 intermediate, and 40
resolved pulsar error envelopes. For this reason, our prediction before even
running the MCMC is that we should not find any strong support at secondary
modes, making this is an ideal candidate for the $H_{0}$ measurement, which
for this set of source and PTA parameters has a predicted
$\mathrm{CV}_{H_{0}}=0.257$. Indeed, this is what we see here – unlike in
Figure 2, our walkers did not seem to find any high probability secondary
modes. Left-hand panel: For the full simulation we initialized multiple
“balls” of walkers as described in Section 4.3. The bottom row shows many
orders of magnitude difference in the values of the log-posterior for the
walkers centered on the true mode versus the walkers away from the true mode.
If we ran the MCMC simulation long enough, we would expect that these low-
probability walkers would eventually find and settle on the true mode. Right-
hand panel: This is the same corner plot, but now with all values of the log-
posterior $\leq-2\times 10^{6}$ removed from the data set. This “cleaned” data
reveals the sharply peaked true values in parameter space, which closely agree
with the Fisher predictions.

It is important to point out that we have shown proof-of-principle in the
ability to disentangle secondary posterior modes from the primary mode. The
caveat, as we see here is that the true mode is sharply peaked, and therefore
may be difficult to find in a blind search over parameter space. Other
techniques may be needed to find these true modes, such as parallel tempering
MCMC, nested sampling, or different jump proposals to improve sampling
efficiency. We leave this open to further investigation in future studies.

##  5 Measurement of the Hubble Constant

###  5.1 Comparing $H_{0}$ Measurements With One versus Two Gravitational Wave
Sources

MC21 and DL21 found that with a single gravitational wave source we must know
a priori the pulsar distances to sub-gravitational wavelength uncertainties in
order to avoid the pulsar distance wrapping problem. Only then can
$D_{\mathrm{par}}$ be recovered with sufficient accuracy to measure the Hubble
constant through equation 7 (exact) or equation 20 (approximate). An example
showing the recovery of $D_{L}$ compared to $D_{\mathrm{par}}$ for a single
source and assuming pulsar distance uncertainties of
$\sigma_{L}=0.1\lambda_{\mathrm{gw}}$ is shown in Figure 5. However, this
level of precision in the pulsar distance measurement is far beyond current
capabilities for all but the most nearby pulsars, making this a very
challenging measurement restricted to gravitational wave sources within $\leq
100$ Mpc (DL21 and MC21).

Figure 5: An example comparison of the uncertainties on the $D_{L}$ and
$D_{\mathrm{par}}$ parameters using a Fisher analysis for a single
gravitational wave source (using the IIB model and default parameters). For
this simulation, we set uniform uncertainties $\sigma=10$ ns and
$\sigma_{L}=0.1\lambda_{\mathrm{gw}}$. This PTA contains 100 pulsars from our
fiducial PTA. The orange region shows the uncertainty range on the luminosity
distance, while the green region shows the uncertainty range on the parallax
distance. For reference, the Fresnel number (equation 5) for a pulsar at the
average pulsar distance $L_{\mathrm{avg}}$ for this particular array is
indicated on the top axis. Generally, the luminosity distance is easier to
measure from the frequency evolution, as compared to the parallax distance
from the Fresnel corrections, which for this example, is measurable for
$D_{\mathrm{par}}\lesssim 400$ Mpc.

The top panel of Figure 6 shows a direct comparison of measuring $H_{0}$ using
one versus two gravitational wave sources. With one gravitational wave source,
when performing a Fisher matrix analysis the only way to accurately invert the
Fisher matrix requires that we first add in pulsar distance priors which would
constrain the uncertainty on the pulsar distances down to the order of the
wrapping cycle (see Section 3.4). Otherwise, the Fisher matrix is ill-
conditioned, due to the strong covariances introduced from the pulsar
distances and parallax distance parameters. So for a single source, the blue
line in the top panel of Figure 6 shows the recovery of $H_{0}$ as a function
of the pulsar distance prior, which is a constraint applied to all of the
pulsars in the PTA uniformly. We see that $H_{0}$ can only be recovered in
this example for distance priors smaller than a few per cent of the
gravitational wave wavelength.

With two gravitational wave sources, we can achieve the same level of accuracy
in $H_{0}$ recovery without pulsar distance prior knowledge. The orange lines
in the top panel of Figure 6 demonstrate this for both a “good” and “poor” sky
location. Simply having two sources at favorable sky separations and
frequencies results in recovery of the Hubble constant better than a single
source with sub-gravitational wavelength prior knowledge. Even unfavorable
source sky positions still result in better recovery than with a single
source.

The bottom panel of Figure 6 demonstrates $H_{0}$ recovery for all sky
realizations of the relative position of the second source, assuming the first
source is located at the position of the red star. In the limit that the two
sources become perfectly aligned with each other on the sky, the recovery of
$H_{0}$ becomes effectively impossible. This is a rather interesting result,
because from our previous discussions in Sections 4.1 and 4.2 we originally
predicted that as long as the source frequency of both sources were different
(even for the same sky position), the wrapping cycle degeneracy should break.
Mathematically this still happens, but when testing various scenarios we found
that in the case of perfect source sky alignment the model parameters for our
PTA and sources needed to be highly favorable, making it very unlikely that
such a circumstance would occur naturally.

Hence, we see that the Fisher matrix analysis demonstrates an important
result, namely that two sources remove the need for prohibitively precise
pulsar distance measurements a priori in order to measure the Hubble constant.
Our one source simulation only obtained values of $\mathrm{CV}_{H_{0}}\sim
0.1$ when all pulsars distances were known with precision
$\sigma_{L}<0.01\lambda_{\mathrm{gw}}$, which is far beyond our current
capabilities. However, two sources without any prior pulsar distance knowledge
resulted in $H_{0}$ measurements ranging from a percent to 10’s of percent,
dependent on the relative sky positions of the two sources.

Figure 6: Fisher analysis predictions of the recovery of $H_{0}$ (using the
IIB model). Both sources have $\mathcal{M}_{\mathrm{obs}}=2\times 10^{9}\
\mathrm{M}_{\odot}$, default parameters otherwise, and the PTA contains 100
pulsars from the fiducial PTA with $\sigma=1$ ns timing uncertainty. Top
panel: A comparison of the recovery of the Hubble constant using one vs. two
gravitational wave sources. For one source, the recovery of $H_{0}$ is plotted
as a function of the required pulsar distance prior knowledge, set uniformly
on all pulsars in the PTA. For two sources, no prior knowledge on the pulsar
distances is assumed. Instead, recovery of $H_{0}$ is shown as a function of
the second source’s frequency. Both changing the second source’s frequency as
well as its sky position relative to the first source breaks the pulsar
distance wrapping cycle degeneracy, so here we show two cases: one where the
second source is in a “good” sky position ($\theta_{2}=\frac{\pi}{3}$ rad,
$\phi_{2}=\frac{2\pi}{3}$ rad), and one where it is in a “poor” sky position
($\theta_{2}=\frac{\pi}{3}$ rad, $\phi_{2}=\frac{5\pi}{3}$ rad). Both of these
locations are marked on the sky plot in the bottom panel – these two positions
lie on the same latitude, but on different parts of the sky. Bottom panel:
$H_{0}$ recovery with two sources, as a function of the second source’s sky
position. Here we see a degenerate region centered on the first source (red
star) and reflected about the galactic plane. In the limit that the two
sources become perfectly aligned the recovery of the Hubble constant becomes
very difficult, even when the two sources have different frequencies.

Figure 7 shows Fisher analysis surveys of the $H_{0}$ measurability in terms
of intrinsic source chirp mass and frequency, as well as PTA characteristics
like timing accuracy and number of pulsars. Not surprisingly, we measure
$H_{0}$ best from sources with high chirp mass and frequency, since these
produce strong frequency evolution effects in the signal. For reference, we
include a rough estimate of the NANOGrav 11 yr continuous wave strain upper
limit, $h_{0,\mathrm{11yr}}\approx 10^{-6}\frac{\omega_{0}}{\pi(1\
\mathrm{Hz})}$ (see fig. 3 of Aggarwal et al., 2019), and indicate the 1 kyr
coalescence time $\Delta\tau_{c}$ contour. These two lines give a sense of
what part of parameter space is most interesting. The lightly shaded region
where $\Delta\tau_{c}<1$ kyr begins to break the original model assumptions,
where frequency evolution becomes very significant (for reference, see
assumption xiv of MC21). The NANOGrav upper limit suggests that sources above
the line would have already been seen (out to $D_{L}=102.35$ Mpc, which is
where the distance is fixed in this plot) in the data if they existed. This
leaves us with a portion of parameter space for two sources with
$\mathcal{M}_{\mathrm{obs}}\lesssim 2\times 10^{9}\ \mathrm{M}_{\odot}$ and
frequencies $\omega_{0,\mathrm{obs}}\gtrsim 20$ nHz where we may be able to
find sources that could recover a measurement of $H_{0}$ in the future.

Figure 7: Fisher analysis predictions for the recovery of $H_{0}$ (using the
IIB model) with two sources, as a function of the observed chirp mass and
orbital frequency of the source, and as a function PTA timing uncertainty and
number of pulsars. For simplicity, both sources are assumed to be identical
(signified by the “$1/2$” notation) except for being located on different
parts of the sky. Source 1 is at the default location, while source 2 has sky
angles $\theta_{2}=\frac{\pi}{3}$ rad and $\phi_{2}=\frac{\pi}{5}$ rad. Left-
hand panel: Here 100 pulsars from our fiducial PTA are timed with $\sigma=1$
ns uncertainty. Hypothetical sources with $\Delta\tau_{c}<1$ kyr have been
shaded out as these represent highly chirping systems (see the general
discussion of assumptions in MC21). Additionally the results of Aggarwal et
al. (2019) are used to include a rough estimate of NANOGrav’s 11 yr continuous
wave strain upper limit, $h_{0,\mathrm{11yr}}\approx
10^{-6}\frac{\omega_{0}}{\pi(1\ \mathrm{Hz})}$, as reference. Right-hand
panel: Here both sources have $\mathcal{M}_{\mathrm{obs}}=10^{10}\
\mathrm{M}_{\odot}$. From our fiducial PTA which has constant timing
uncertainty, we see that the recovery of $H_{0}$ scales as
$\mathrm{CV}_{H_{0}}\propto\sigma/\sqrt{N_{p}}$. For variable timing
uncertainty, the SKA PTA no longer achieves the same improvement by adding
more pulsars, but rather now scales as
$\mathrm{CV}_{H_{0}}\propto\sigma_{\mathrm{min}}$ at large $N_{p}$.

The right-hand panel of Figure 7 shows that timing precision makes a
significant difference in improving the measurement of $H_{0}$. For the
fiducial PTAs, we gain an order of magnitude improvement in
$\mathrm{CV}_{H_{0}}$ for the same improvement in $\sigma$. Additionally, we
see simply adding more fiducial PTA pulsars (randomly drawn from our
distribution equation 22) improves the entire network’s ability to recover
$H_{0}$. Overall, recovery of $H_{0}$ for the fiducial PTAs with constant
timing uncertainties scales as $\mathrm{CV}_{H_{0}}\sim\sigma/\sqrt{N_{p}}$.

Also plotted in Figure 7 are the $\mathrm{CV}_{H_{0}}$ recovery uncertainties
found from various SKA-era PTAs, with more realistic timing uncertainties (see
Section 3.3.2). The upright triangles (optimistic) represent PTAs for which
$\sigma\propto\mathrm{SNR}^{-1/2}$ while the upside-down triangles
(pessimistic) represent PTAs for which $\sigma\propto\mathrm{SNR}^{-1}$
(equation 24). Grey (black) markers represent PTAs for which the best-timed
pulsars have $\sigma=1$ ns (10 ns). For small number arrays, the scaling with
$N_{p}$ approaches the $1/\sqrt{N_{p}}$ scaling of the constant $\sigma$ case.
For larger $N_{p}$ this relation limits towards a constant value. This is due
to the way in which we add pulsars to the array, in ranked order of SNR per
radial bin, and to the limited number of high SNR pulsars at large distance.
As $N_{p}$ increases, we exhaust the supply of high-SNR pulsars at large
distance from Earth so while $N_{p}$ increases, $N_{p}$ at large distance and
high SNR effectively does not, and these are the pulsars which mostly strongly
constrain $D_{\mathrm{par}}$ and hence $H_{0}$. Given these specific SKA-era
PTAs, we find that arrays composed of $\sim 200$ pulsars timed to a best case
precision of $\sigma_{\min}=1$ ns (10 ns) could measure the Hubble constant to
within $2\%$ ($20\%$) precision using the two-source method described here.

###  5.2 Prospects for Future SKA-like PTAs

To further explore the prospects of using future SKA-like PTAs to measure the
Hubble constant, we examine the possibility that a detectable SMBHB source
exists. Because closer sources within a few hundred Mpcs will have larger
Fresnel effects, and hence are more likely to produce a measurement of
$H_{0}$, we take an observational approach and consider the catalogue of PTA
constraints on gravitational waves from massive SMBHBs in galaxies within 600
Mpc generated from the NANOGrav 11 yr data set in Arzoumanian et al.
(2021).333We used the data from Table 3 – the “Mass,” “Dist,” and “$q_{95}$
(10.0)” columns. This includes 216 nearby galaxies for which PTA upper limits
can constrain the mass ratios of putative SMBHBs at a given gravitational wave
frequency (we consider the $f_{\mathrm{gw}}=10$ nHz, circular orbit
constraints). In the left-hand panel of Figure 8 we plot these limiting
putative SMBHBs in chirp mass and distance space. Using a best-case, 438
pulsar, SKA-era PTA of Section 3.3.2 ($5\%$ of all detected MSPs), we overlay
contours of the corresponding $H_{0}$ recoverability $\mathrm{CV}_{H_{0}}$.
For this optimistic SKA-era PTA, there are 10’s of possible nearby sources.

The right-hand panel of Figure 8 explores this further by taking each putative
SMBHB (red dots) from the left-hand panel, assuming it has a twin at a
different position on the sky, and computing the $H_{0}$ recovery uncertainty
$\mathrm{CV}_{H_{0}}$. The cumulative distribution reveals how increasing the
size of our PTA increases the number of twin SMBHBs for which $H_{0}$ recovery
at a given precision is possible. With 438 pulsars or more, we find 10’s of
systems with $\mathrm{CV}_{H_{0}}\leq 0.5$, and over $30\%$ of all of these
systems procure $\mathrm{CV}_{H_{0}}\leq 0.9$. Over $5\%$ (i.e. more than
$\sim 12$) of these putative systems result in $\mathrm{CV}_{H_{0}}\leq 0.4$,
and there are a few cases where we recover $\mathrm{CV}_{H_{0}}\sim 0.1$.

Hence, an $\mathcal{O}(10\%)$, purely gravitational wave $H_{0}$ measurement
is possible if two of these optimal putative SMBHBs exist. Note that already
for the ‘pessimistic,’ $\sigma_{\min}=1$ ns SKA-era arrays (see the right-hand
panel of Figure 7), the overlap between putative sources and
$\mathrm{CV}_{H_{0}}\leq 1$ in Figure 8 becomes marginal. Hence these best-
case PTAs are likely required.

Figure 8:  Left-hand panel: Fisher analysis predictions of the recovery of
$H_{0}$ (using the IIB model) with two sources using 438 ($5\%$) of our
“optimistic” SKA-era pulsars with $\sigma_{\mathrm{min}}=1$ ns, as a function
of the observed chirp mass and parallax distance of the source. In the contour
plot, for simplicity, both SMBHB sources are assumed to be identical (hence
the “$1/2$” notation) except for their sky location. Source 1 is at the
default location, while source 2 has sky angles $\theta_{2}=\frac{\pi}{3}$ rad
and $\phi_{2}=\frac{\pi}{5}$ rad. Additionally all sources are simulated with
$f_{\mathrm{gw}}=10$ nHz (i.e. $\omega_{0,\mathrm{obs}}\approx 31.4$ nHz).
Plotted as red dots on top of the contours, we compare 216 hypothetical SMBHB
sources within 600 Mpc published in Arzoumanian et al. (2021). This gives us a
sense of the quality required of an SKA-era PTA in order to make the $H_{0}$
measurement. For the same number of pulsars, the pessimistic version of this
PTA resulted in a marginal number of sources with $CV_{\mathrm{H_{0}}}<1$.
Right-hand panel: For each of the 216 SMBHB sources indicated on the left-hand
panel, we simulate two identical SMBHBs with those mass and distance values
(and with the indicated SKA-era PTAs), and bin them based on their level of
recovery of $H_{0}$ in a cumulative distribution. Both SMBHB sources are
assumed to be identical with the same caveats as the left-hand panel. As
expected from the right-hand panel of Figure 7, adding more pulsars into the
PTA improves the recovery of $H_{0}$ in more two-system scenarios. For the
PTAs with $>173$ pulsars, we begin to find double sources where
$\mathrm{CV}_{H_{0}}\leq 0.3$ (comparable to what is shown in Figure 4). These
plots do not rule out the possibility of a measurement of $H_{0}$ with a
future SKA PTA, and instead suggests that there do exist hypothetical SMBHB
candidates wherein this measurement may be viable.

##  6 Conclusions

In this work we have shown that future PTA experiments could make purely
gravitational wave-based measurements of the Hubble constant. This is made
possible by accounting for the Fresnel curvature effects in the wavefront
across the Earth-pulsar baseline. By using the fully general Fresnel frequency
evolution timing residual model, we can obtain two separate distance
measurements to the source: the luminosity distance $D_{L}$ (from the
frequency evolution effects) and the parallax distance $D_{\mathrm{par}}$
(from the Fresnel effects).

The measurement of the parallax distance is particularly challenging due to
the pulsar distance wrapping problem. Unless the distances to the pulsars in
our array can be measured to sub-gravitational wavelength ($\sim$sub-parsec)
precision, this measurement cannot be made with a single SMBHB source.
However, we demonstrate that two SMBHB sources detected simultaneously break
the wrapping cycle degeneracy, allowing a viable scenario for measuring the
Hubble constant. Two or more sources will actively calibrate the PTA by
helping to identify the true pulsar distances from the degenerate secondary
modes. The ability to better resolve the pulsar distances using this method
will be explored in a follow-up study (blueMcGrath, D’Orazio, & Creighton, in
preparation). Additionally, this method may better improve the recovery of the
luminosity distance $D_{L}$ (to both sources), even at higher redshifts where
the Fresnel effects may be too small to measure $D_{\mathrm{par}}$. It would
be interesting for future studies to investigate how significant this
improvement is, compared to measuring $D_{L}$ of just a single source, since
the luminosity distance by itself would still be useful for other
astrophysical and cosmological studies.

For the two source problem, we explored the measurement of $H_{0}$ for a wide
range of systems and PTAs. We developed new Fisher matrix-based criteria which
can quickly predict the ability of a particular source-PTA set-up to measure
$H_{0}$, and facilitate follow-up with MCMC verification. While our MCMC
techniques are less complicated than other samplers for PTA problems, we find
they support our proof-of-principle calculations, and we leave it to future
work to improve upon methods for efficiently extracting the model parameters
(see for example, Samajdar et al., 2022).

With fiducial PTAs, with constant timing uncertainty $\sigma$ and $N_{p}$
pulsars, the precision in our simulated $H_{0}$ measurement scales as
$\mathrm{CV}_{H_{0}}\sim\sigma/\sqrt{N_{p}}$. With more realistic SKA-era
PTAs, we observe a similar scaling for arrays with less than a few hundred
pulsars and a saturation in measurement precision for larger arrays, where
$\mathrm{CV}_{H_{0}}\sim\sigma_{\mathrm{min}}$. We also find that hypothetical
two-source SMBHB systems, based on the NANOGrav 11 yr data in Arzoumanian et
al. (2021), suggest that SKA-era PTAs may be capable of measuring $H_{0}$ to
as low as $\mathrm{CV}_{H_{0}}\sim\mathcal{O}\left(0.1\right)$. This will
likely require fortuitous SMBHB sources and “best-case” PTAs (such that
$\sigma_{\mathrm{min}}\sim 1$ ns to $10$ ns) containing several hundred
pulsars. While we have not considered the possibility of more than two sources
here, this will likely improve measurement prospects.

In conclusion, the existence of two individually resolvable gravitational wave
signals from inspiraling SMBHBs allows a purely gravitational wave-based
measurement of the Hubble constant with PTAs in the SKA-era and beyond.
Current and future $H_{0}$ measurements can achieve higher precision than what
we envision here on that time-scale; for example, the standard sirens approach
using binary neutron star mergers with electromagnetic counterparts could
reach a 1% level $H_{0}$ measurement by the 2030s (Chen et al., 2018).
However, as is apparent from the current Hubble tension, precision does not
guarantee agreement between measurement methods. The value of the method
discussed here, using only the gravitational messenger and geometry, is that
it can provide a novel and independent measurement of $H_{0}$.

## Acknowledgements

This material is based upon work supported by the National Aeronautics and
Space Administration (NASA) under award number 80GSFC21M0002. This work was
also supported by the National Science Foundation (NSF) PHY-1430284 [through
the North American Nanohertz Observatory for Gravitational Waves (NANOGrav’s)
Physics Frontier Center], PHY-1912649, and PHY-2207728, and by UW-Milwaukee’s
computational resources PHY-1626190. DJD received funding from the European
Union’s Horizon 2020 research and innovation programme under Marie Sklodowska-
Curie grant agreement No. 101029157, and from the Danish Independent Research
Fund through Sapere Aude Starting Grant No. 121587. We thank the referee, Neil
Cornish, for a constructive report.

## Data Availability

Calculations in this paper were performed using code developed by the authors,
and the data are available on reasonable request to the authors.

## References

  * Aggarwal et al. (2019) Aggarwal K., et al., 2019, ApJ, 880, 116
  * Arzoumanian et al. (2018) Arzoumanian Z., et al., 2018, ApJS, 235, 37
  * Arzoumanian et al. (2021) Arzoumanian Z., et al., 2021, ApJ, 914, 121
  * Babak & Sesana (2012) Babak S., Sesana A., 2012, Phys. Rev. D, 85, 044034
  * Bates et al. (2014) Bates S. D., Lorimer D. R., Rane A., Swiggum J., 2014, MNRAS, 439, 2893
  * Caldwell (1993) Caldwell R. R., 1993, Phys. Rev. D, 48, 4688
  * Chen et al. (2018) Chen H.-Y., Fishbach M., Holz D. E., 2018, Nature, 562, 545
  * Corbin & Cornish (2010) Corbin V., Cornish N. J., 2010, arXiv e-prints, p. arXiv:1008.1782
  * Cordes & Chernoff (1997) Cordes J. M., Chernoff D. F., 1997, ApJ, 482, 971
  * Creighton & Anderson (2011) Creighton J. D. E., Anderson W. G., 2011, Gravitational-Wave Physics and Astronomy: An Introduction to Theory, Experiment and Data Analysis. Wiley-VCH
  * D’Orazio & Loeb (2021) D’Orazio D. J., Loeb A., 2021, Phys. Rev. D, 104, 063015
  * Deller et al. (2019) Deller A. T., et al., 2019, ApJ, 875, 100
  * Deng & Finn (2011) Deng X., Finn L. S., 2011, MNRAS, 414, 50
  * Feeney et al. (2019) Feeney S. M., Peiris H. V., Williamson A. R., Nissanke S. M., Mortlock D. J., Alsing J., Scolnic D., 2019, Physical Review Letters, 122, 061105
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Ghosh et al. (2022) Ghosh T., Biswas B., Bose S., 2022, arXiv e-prints, p. arXiv:2203.11756
  * Hogg (1999) Hogg D. W., 1999, arXiv e-prints, pp astro–ph/9905116
  * Holz & Hughes (2005) Holz D. E., Hughes S. A., 2005, ApJ, 629, 15
  * Iacovelli et al. (2022) Iacovelli F., Mancarella M., Foffa S., Maggiore M., 2022, arXiv e-prints, p. arXiv:2207.02771
  * Kelley et al. (2018) Kelley L. Z., Blecha L., Hernquist L., Sesana A., Taylor S. R., 2018, MNRAS, 477, 964
  * Lee et al. (2011) Lee K. J., Wex N., Kramer M., Stappers B. W., Bassa C. G., Janssen G. H., Karuppusamy R., Smits R., 2011, MNRAS, 414, 3251
  * Lorimer et al. (2006) Lorimer D. R., et al., 2006, MNRAS, 372, 777
  * Maggiore (2008) Maggiore M., 2008, Gravitational Waves: Theory and Experiments. Vol. 1, Oxford University Press
  * Maggiore (2018) Maggiore M., 2018, Gravitational Waves: Astrophysics and Cosmology. Vol. 2, Oxford University Press
  * McGrath & Creighton (2021) McGrath C., Creighton J., 2021, MNRAS,
  * Messenger & Read (2012) Messenger C., Read J., 2012, Phys. Rev. Lett., 108, 091101
  * Petiteau et al. (2013) Petiteau A., Babak S., Sesana A., de Araújo M., 2013, Phys. Rev. D, 87, 064036
  * Qian et al. (2022) Qian Y.-Q., Mohanty S. D., Wang Y., 2022, Phys. Rev. D, 106, 023016
  * Samajdar et al. (2022) Samajdar A., et al., 2022, arXiv e-prints, p. arXiv:2205.04332
  * Schneider (2015) Schneider P., 2015, Extragalactic Astronomy and Cosmology: An Introduction, 2nd edn. Springer, p. 177–192, doi:10.1007/978-3-642-54083-7
  * Schutz (1986) Schutz B. F., 1986, Nature, 323, 310
  * Shiralilou et al. (2022) Shiralilou B., Raaijmakers G., Duboeuf B., Nissanke S., Foucart F., Hinderer T., Williamson A., 2022, arXiv e-prints, p. arXiv:2207.11792
  * Smits et al. (2009) Smits R., Kramer M., Stappers B., Lorimer D. R., Cordes J., Faulkner A., 2009, A&A, 493, 1161
  * The LIGO, Virgo, & KAGRA Collaborations (2021) The LIGO, Virgo, & KAGRA Collaborations 2021, arXiv e-prints, p. arXiv:2111.03604
  * The LIGO, Virgo, 1M2H, Dark Energy Camera GW-EM, & DES Collaborations (2017) The LIGO, Virgo, 1M2H, Dark Energy Camera GW-EM, & DES Collaborations 2017, Nature, 551, 85
  * Verbiest & Shaifullah (2018) Verbiest J. P. W., Shaifullah G. M., 2018, Classical and Quantum Gravity, 35, 133001
  * Wang et al. (2022) Wang L.-F., Zhang G.-P., Shao Y., Zhang X., 2022, arXiv e-prints, p. arXiv:2201.00607
  * Weinberg (1972) Weinberg S., 1972, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity. New York: Wiley

##  Appendix A Cosmological Distances

For a gravitational wave propagating from our source (at comoving coordinate
distance $R$ and redshift $z$) to the Earth (“infalling”), we define the
following “boundary” conditions:

Source: $\displaystyle\quad r=R,\quad z^{\prime}=z,\quad
t^{\prime}=t_{\mathrm{ret}},$ Earth: $\displaystyle\quad r=0,\quad
z^{\prime}=0,\quad t^{\prime}=t,$

where our reference frame is centered on the Earth. From the Friedmann
equations and the definition of the redshift in terms of the present day
$t_{0}$ scale factor, $1+z=\frac{a(t_{0})}{a(t)}$, we explicity list the
following important quantities (Weinberg, 1972; Hogg, 1999; Schneider, 2015):

Expansion Rate $\displaystyle\qquad H$ $\displaystyle\equiv\ $
$\displaystyle\left(\frac{\dot{a}(t)}{a(t)}\right)=H_{0}E,$ (31) Dimensionless
Hubble Function $\displaystyle\qquad E$ $\displaystyle\equiv\ $
$\displaystyle\left[\frac{\Omega_{r}}{a^{4}(t)}+\frac{\Omega_{m}}{a^{3}(t)}+\Omega_{\Lambda}+\frac{\Omega_{k}}{a^{2}(t)}\right]^{1/2},$
$\displaystyle=\ $
$\displaystyle\left[\frac{(1+z)^{4}}{a^{4}(t_{0})}\Omega_{r}+\frac{(1+z)^{3}}{a^{3}(t_{0})}\Omega_{m}+\Omega_{\Lambda}+\frac{(1+z)^{2}}{a^{2}(t_{0})}\Omega_{k}\right]^{1/2},$
(32) Curvature Density $\displaystyle\qquad\Omega_{k}$ $\displaystyle\equiv\ $
$\displaystyle
a^{2}(t_{0})\left[1-\frac{\Omega_{r}}{a^{4}(t_{0})}-\frac{\Omega_{m}}{a^{3}(t_{0})}-\Omega_{\Lambda}\right]=-kD_{H}^{2},$
(33) Hubble Distance $\displaystyle\qquad D_{H}$ $\displaystyle\equiv\ $
$\displaystyle\frac{c}{H_{0}},$ (34) Line-of-Sight Comoving Distance
$\displaystyle\qquad D_{c}$ $\displaystyle\equiv\ $ $\displaystyle
D_{H}\int^{z}_{0}\frac{\mathrm{d}z^{\prime}}{E},$ $\displaystyle=\ $
$\displaystyle
ca(t_{0})\int^{t}_{t_{\mathrm{ret}}}\frac{\mathrm{d}t^{\prime}}{a},$ (35)
Parallax Distance $\displaystyle\qquad D_{\mathrm{par}}$ $\displaystyle\equiv\
$ $\displaystyle a(t_{0})\frac{R}{\sqrt{1-kR^{2}}},$ (36)

where $H_{0}$ is the Hubble constant and $\Omega_{r}$, $\Omega_{m}$, and
$\Omega_{\Lambda}$ are the radiation, matter, and vacuum density parameters,
respectively. For completeness here we will continue to write the present day
scale factor explicitly, but typical convention is to normalize this to
$a(t_{0})\equiv 1$. The two forms of the line-of-sight comoving distance in
equation 35 are made through the change of variables $t\rightarrow z$ using
the definition of $z$ and equation 31.

The luminosity distance of our source is defined by considering the energy
flux $F_{\mathrm{obs}}$ (in this case, of our gravitational waves) as measured
by the observer (Maggiore, 2008). In a cosmologically expanding universe the
observed energy $E_{\mathrm{obs}}$ is redshifted compared to the energy that
was emitted in the source frame $E_{\mathrm{s}}$, such that
$E_{\mathrm{s}}=(1+z)E_{\mathrm{obs}}$. Additionally, time dilation means the
time measured by the observer and source clocks are related by
$t_{\mathrm{s}}=\frac{t_{\mathrm{obs}}}{(1+z)}$. Finally, from the FLRW metric
equation 1, the flux at time $t$ will be spread over a total area of
$A(t)=4\pi a^{2}(t)R^{2}$. This means that we can write:

$F_{\mathrm{obs}}(t)\ \ \equiv\ \ \frac{L_{\mathrm{obs}}}{A(t)}\ \ =\ \
\frac{L_{\mathrm{s}}}{(1+z)^{2}4\pi a^{2}(t)R^{2}}\ \ =\ \
\frac{L_{\mathrm{s}}}{4\pi\Big{[}(1+z)^{2}a^{2}(t)R^{2}\Big{]}}\ \ \equiv\ \
\frac{L_{\mathrm{s}}}{4\pi D_{L}^{2}},$ (37)

where in the final equality we define our luminosity distance $D_{L}$ such
that the observed flux matches our standard notion of flux, but as a function
of the source frame luminosity $L_{\mathrm{s}}$. Therefore at the present
time, we have:

$D_{L}\equiv(1+z)a(t_{0})R.$ (38)

Next we draw the connection between the comoving coordinate distance $R$, the
luminosity distance $D_{L}$, and the line-of-sight comoving distance $D_{c}$.
Integrating the metric equation 1 along the radial
($\mathrm{d}\theta=\mathrm{d}\phi=0)$, null ($\mathrm{d}s=0$), infalling path
of the gravitational wave, and using equation 35 gives us:

$\displaystyle-\int^{0}_{R}\frac{\mathrm{d}r}{\sqrt{1-kr^{2}}}$
$\displaystyle=\int^{t}_{0}\frac{c}{a(t^{\prime})}\mathrm{d}t^{\prime},$
$\displaystyle=\frac{D_{c}}{a(t_{0})},$

The solution of the left-hand integral depends on the sign of the curvature
constant $k$ (or alternatively, the density parameter $\Omega_{k}$). Solving
this integral we now have:

$R\ \ \equiv\ \ D_{M}\ \ =\ \
\left.\begin{cases}\frac{1}{\sqrt{k}}\sin\left(\frac{\sqrt{k}}{a(t_{0})}D_{c}\right),&\quad
k>0\\\\[4.0pt] \frac{1}{a(t_{0})}D_{c},&\quad k=0\\\\[4.0pt]
\frac{1}{\sqrt{|k|}}\sinh\left(\frac{\sqrt{|k|}}{a(t_{0})}D_{c}\right),&\quad
k<0\end{cases}\right\\}\ \ =\ \
\left.\begin{cases}\frac{D_{H}}{\sqrt{|\Omega_{k}|}}\sin\left(\frac{\sqrt{|\Omega_{k}|}}{a(t_{0})D_{H}}D_{c}\right),&\quad\Omega_{k}<0\\\\[4.0pt]
\frac{1}{a(t_{0})}D_{c},&\quad\Omega_{k}=0\\\\[4.0pt]
\frac{D_{H}}{\sqrt{\Omega_{k}}}\sinh\left(\frac{\sqrt{\Omega_{k}}}{a(t_{0})D_{H}}D_{c}\right),&\quad\Omega_{k}>0\end{cases}\right\\}\
\ =\ \ \frac{D_{L}}{a(t_{0})(1+z)},$ (39)

where $D_{M}$ is introduced as the “transverse comoving distance” following
Hogg, which is just a re-naming of the comoving coordinate distance. Thus
depending on the global universe geometry, the relationship between the
various distances $R$, $D_{L}$, $D_{c}$, and $D_{M}$ are given by equation 39.

One final and very important distance for our consideration is the “parallax
distance” $D_{\mathrm{par}}$ in equation 36 (Weinberg, 1972). Substituting
equation 39 into equation 36 lets us write:

$D_{\mathrm{par}}\ \ \equiv\ \ \frac{a(t_{0})R}{\sqrt{1-kR^{2}}}\ \ =\ \
\frac{a(t_{0})D_{M}}{\sqrt{1-kD_{M}^{2}}}\ \ =\ \
\left.\begin{cases}\frac{a(t_{0})}{\sqrt{k}}\tan\left(\frac{\sqrt{k}}{a(t_{0})}D_{c}\right),&\quad
k>0\\\\[4.0pt] D_{c},&\quad k=0\\\\[4.0pt]
\frac{a(t_{0})}{\sqrt{|k|}}\tanh\left(\frac{\sqrt{|k|}}{a(t_{0})}D_{c}\right),&\quad
k<0\end{cases}\right\\}\ \ =\ \
\frac{D_{L}}{(1+z)\sqrt{1-\frac{kD_{L}^{2}}{a^{2}(t_{0})(1+z)^{2}}}}.$ (40)

Notably in a flat universe (with $a(t_{0})=1$ normalization) we have
$D_{\mathrm{par}}=R=D_{M}=D_{c}=D_{L}/(1+z)$.

As a final note, the closed universe case needs further special consideration,
as the parallax distance not only diverges, but can take on negative values.
For closed geometry, the expression for $D_{\mathrm{par}}$ in equations 36 and
40 has a $\pm$ sign in front of the expression. Conceptually, consider a point
at the north pole of a globe generating some wave. As the wave propagates away
from the source, in the northern hemisphere two points on the wavefront
(constant latitude) would produce a positive parallax measurement
$D_{\mathrm{par}}>0$ back to the source (the wavefront curves towards the
source and away from the direction of propagation). At the equator, however,
all points along the wavefront will experience a plane-wave. Therefore there
is no well-defined parallax distance for points along the equator, and
mathematically we see that the value of $D_{\mathrm{par}}\rightarrow\infty$
(i.e. at $R=R_{\mathrm{max}}=\frac{1}{\sqrt{k}}$, or equivalently
$\chi=\frac{1}{\sqrt{k}}\frac{\pi}{2}$ through equation 2).

Then interestingly, in the southern hemisphere two points on the wavefront
(constant latitude) would now produce a negative parallax measurement
$D_{\mathrm{par}}$ back to the source (the wavefront curves away the source
and towards the direction of propagation – towards the source antipode). Using
equation 2 to replace $R$ with $\chi$ in equation 40, and including the $\pm$
sign out front, we write:

$D_{\mathrm{par}}=\pm\frac{a(t_{0})\sin\left(\sqrt{k}\chi\right)}{\sqrt{k}\sqrt{1-\sin^{2}\left(\sqrt{k}\chi\right)}}=\pm\frac{a(t_{0})\sin\left(\sqrt{k}\chi\right)}{\sqrt{k}\sqrt{\cos^{2}\left(\sqrt{k}\chi\right)}}=\frac{a(t_{0})}{\sqrt{k}}\tan\left(\sqrt{k}\chi\right).\qquad(\text{for
}k>0)$

Writing $D_{\mathrm{par}}$ in terms of $\chi$ takes care of the sign ambiguity
in equations 36 and 40. For values $0\leq\chi<\frac{1}{\sqrt{k}}\frac{\pi}{2}$
(i.e. the northern hemisphere) the parallax distance is positive, for
$\chi=\frac{1}{\sqrt{k}}\frac{\pi}{2}$ (i.e. the equator) the parallax
distance diverges, and for values
$\frac{1}{\sqrt{k}}\frac{\pi}{2}<\chi\leq\frac{1}{\sqrt{k}}\pi$ (i.e. the
southern hemisphere) the parallax distance is negative.

##  Appendix B The Retarded Time Calculation

The goal is to calculate $\left|\vec{x}^{\prime}-\vec{x}\right|$ in equation
3, where $\vec{x}^{\prime}$ is the source location and $\vec{x}$ is the field
point. For our problem, we will primarily be interested in setting our field
point at the pulsar. With the center of our coordinate system at the Earth,
let the source be at the coordinate $r=R$ (i.e. $\chi=\chi_{S}$) with unit
vector $\hat{r}$, and let the pulsar be at the coordinate $r=L$ (i.e.
$\chi=\chi_{P}$) with unit vector $\hat{p}$.

Next we write the law of cosines for Euclidean, spherical, and hyperbolic
geometries:

$\begin{cases}\begin{tabular}[]{r l
r}$\cos\left(\sqrt{k}\left|\vec{x}^{\prime}-\vec{x}\right|\right)$&$=\cos\left(\sqrt{k}\chi_{S}\right)\cos\left(\sqrt{k}\chi_{P}\right)+\sin\left(\sqrt{k}\chi_{S}\right)\sin\left(\sqrt{k}\chi_{P}\right)\left(\hat{r}\cdot\hat{p}\right),$&$k>0$\\\\[4.0pt]
$\left|\vec{x}^{\prime}-\vec{x}\right|$&$=\chi_{S}^{2}+\chi_{P}^{2}-2\chi_{S}\chi_{P}\left(\hat{r}\cdot\hat{p}\right),$&$k=0$\\\\[4.0pt]
$\cosh\left(\sqrt{|k|}\left|\vec{x}^{\prime}-\vec{x}\right|\right)$&$=\cosh\left(\sqrt{|k|}\chi_{S}\right)\cosh\left(\sqrt{|k|}\chi_{P}\right)-\sinh\left(\sqrt{|k|}\chi_{S}\right)\sinh\left(\sqrt{|k|}\chi_{P}\right)\left(\hat{r}\cdot\hat{p}\right).$&$k<0$\end{tabular}\end{cases}$
(41)

Recall from Section 2.2 that at the present time the Gaussian curvature is
$k$. One approach to solving these equations is to Taylor expand in the small
parameter $\sqrt{|k|}\chi_{P}$ in the closed and open universe cases, and in
the small parameter $\frac{\chi_{P}}{\chi_{S}}$ the flat universe case. The
result is:

$\left|\vec{x}^{\prime}-\vec{x}\right|\approx\chi_{S}-\left(\hat{r}\cdot\hat{p}\right)\chi_{P}+\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\begin{cases}\frac{\chi_{P}^{2}}{\frac{1}{\sqrt{k}}\tan\left(\sqrt{k}\chi_{S}\right)},\quad&k>0\\\\[12.0pt]
\frac{\chi_{P}^{2}}{\chi_{S}},\quad&k=0\\\\[4.0pt]
\frac{\chi_{P}^{2}}{{\frac{1}{\sqrt{|k|}}\tanh\left(\sqrt{|k|}\chi_{S}\right)}}.\quad&k<0\end{cases}$
(42)

Additionally, using equation 2 we see that $\chi_{P}\approx L$ for each
universe geometry case (Taylor expanding once again for the closed and open
cases in $\sqrt{|k|}\chi_{P}$).

Consider for a moment a gravitational wave propagating directly from the
source to the Earth. The wave leaves it’s position at redshift $z$ and time
$t_{\mathrm{ret}}$, and arrives at Earth at redshift $z=0$ at time $t$ (recall
the boundary conditions from Appendix A). In this case our field position
$\vec{x}=0$, so $\eta_{\mathrm{ret}}=\eta-\frac{\chi_{S}}{c}$ from equation 3.
Therefore using equation 35 we can write:

$\displaystyle-\frac{\chi_{S}}{c}=\eta_{\mathrm{ret}}-\eta$
$\displaystyle=\int^{t_{\mathrm{ret}}}_{t}\frac{\mathrm{d}t^{\prime}}{a}$
$\displaystyle=-\frac{D_{c}}{a(t_{0})c}.$

This gives us a connection between $\chi_{S}$ and the comoving coordinate
distance $D_{c}$, namely that $\chi_{S}=\frac{D_{c}}{a(t_{0})}$, which we can
now substitute into equation 42. The result is that we can now see that the
relevant distance which appears in the Fresnel term in our Taylor expansion
(for all three universe cases) is the parallax distance $D_{\mathrm{par}}$
(substitute in equation 40). Therefore the final form of equation 42 which is
then used in equation 4 is:

$\left|\vec{x}^{\prime}-\vec{x}\right|\approx\frac{D_{c}}{a(t_{0})}-\left(\hat{r}\cdot\hat{p}\right)L+\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)a(t_{0})\frac{L^{2}}{D_{\mathrm{par}}}.$
(43)

Note that in the flat universe case $D_{\mathrm{par}}=D_{c}$, which is
consistent with the result of DL21 given that this was their working
assumption.

##  Appendix C Assumptions

Assumptions (ii)–(xvi) of MC21 carry over to this work. Additionally, we work
under the following assumptions:

  1. 1.

The universe is described by the FLRW metric with comoving coordinates
$(t,r,\theta,\phi)$ and curvature constant $k$.

  2. 2.

The FLRW scale factor $a$ does not change appreciably over the time-scale of
the travel time of the photon from the pulsar to the Earth:
$\dot{a}(t_{0})\Delta t_{\mathrm{p}\rightarrow\mathrm{E}}$.

  3. 3.

All “local” distances and times between the pulsar and the Earth are on a
Minkowski background,
$g_{\mu\nu}=\eta_{\mu\nu}=\text{diag}\left(-c^{2},1,1,1\right)$. The
transition from the global cosmological FLRW background to the local Minkowski
background requires that all points of interest have comoving coordinate
distance separations much smaller than the background curvature
$\left(r^{2}\ll\frac{1}{|k|}\right)$, and that the observation time of our
experiment is much smaller than the present day age of the universe
$\left(\frac{t_{\mathrm{obs}}}{t_{0}}\ll 1\right)$.

  4. 4.

For the results presented in this work we assume a flat universe ($k=0$).

##  Appendix D Pulsar Populations and Example Pulsar Distance Posteriors

Two pulsar populations are used for different studies within this paper, a
“fiducial PTA” and an “SKA-era PTA”. The general characteristics of these PTAs
are shown in Figure 9.

Figure 9: A side-by-side comparison of the fiducial PTA (left-hand column) and
SKA-era PTA (right-hand column). These figures show the general distribution
properties, in terms of the pulsar distances and locations throughout the
Milky Way Galaxy. In the middle row, Earth is located at the origin and
galactic center is located at $8$ kpc (at $\phi=0^{\circ}$). Left-hand panel:
This fiducial PTA contains 1000 pulsars generated from the distribution in
equation 22. Right-hand panel: This SKA-era PTA contains 1194 pulsars.

To create the SKA-era PTAs, we use the PSRPOPPY package (Bates et al., 2014)
to generate a realistic population of MSPs detectable in the SKA era. To
generate the underlying population of pulsars we assume the default model in
PSRPOPPY (except that we place the galactic center at 8 instead of 8.5 kpc),
with a spatial pulsar distribution following Lorimer et al. (2006), but with
the pulsar period distribution from Cordes & Chernoff (1997), chosen
specifically for MSPs. We then follow (Smits et al., 2009) in assuming that
there are $3\times 10^{4}$ MSPs in the galaxy to generate the discoverable
population.

We simulate an SKA-like survey by choosing survey parameters following Smits
et al. (2009) (model A) and summarized in Table 1. We find that our SKA survey
detects $8893$ of these MSPs with SNR greater than a threshold value of 9, at
distances ranging from $0.2$ kpc out to $18$ kpc. We assume that only some
percentage of these MSPs will be suitable for high precision timing, and
enforce this as described by the selection process in Section 3.3.2. We then
compute a timing uncertainty for each MSP as per equation 24. Timing
uncertainty distributions for the different choices of uncertainty scaling
with SNR are plotted in Figure 10.

Figure 10: Timing uncertainty distributions of all of the MSPs detected in a
mock SKA-era population, with $\sigma_{\min}=1$ ns. The blue circles and
orange crosses denote optimistic and pessimistic choices for $\sigma_{a}$, as
written in equation 24. Mock SKA-era PTAs are generated by choosing MSPs from
this set, as outlined in the text.

Finally, Figure 11 shows the complete set of 1D posteriors for the 40 pulsar
PTA shown in the right-hand panel of Figure 2, which motivates the modal
overlap criterion defined in equation 28.

Figure 11: The complete set of 1D posteriors for the same results shown in the right-hand panel of Figure 2 using the IB model. Here we have sorted the pulsars in ascending value by the modal overlap criterion ratio defined in equation 28 (given above each figure). As we can see, this particular criterion successfully shows the progression of the error envelope from “overlapping” to “resolved,” with a total of 12 overlapping, 18 intermediate, and 10 resolved pulsars. For ratios less than about $4$, we see that uncertainties about each mode are both sufficiently large and sufficiently close to the primary mode that they blend together around the true value. However, for ratios between about $4$ and $8$ these secondary modes begin to separate away from the primary mode far enough that they start to become distinguishable from the true mode. Finally, for ratios greater than about $8$, we see that the secondary modes are distinguishable from the true mode. In practice this therefore motivates us to search our results to see what the value of this ratio is for each of our pulsars in our desired source-PTA simulation. If it is greater than $8$, then we would expect that the true distance to that pulsar would be identifiable. Note the samples in this figure have been cleaned, such that all samples with log-posterior $\leq-100$ have been removed from the data set (see the discussion of cleaning in Section 4.3). Table 1: SKA survey parameters motivated by Smits et al. (2009) (their Model A), used to generate mock PTAs via the PSRPOPPY code (Bates et al., 2014). Name |  Value
---|---
Survey degradation factor |  1.0
Antenna gain (K $\mathrm{Jy}^{-1}$) |  130
Integration time (s) |  1800
Sampling time (ms) |  0.1
System temperature (K) |  30
Centre frequency (MHz) |  1400
Bandwidth (MHz) |  500
Channel bandwidth (MHz) |  0.009
No. polarizations |  2
FWHM (arcmin) |  65.5
Min RA (deg) |  0
Max RA (deg) |  360
Min DEC (deg) |  -90
Max DEC (deg) |  30
Frac. survey coverage. |  1
SNR threshold |  9