# Fresnel Models for Gravitational Wave Effects on Pulsar Timing

Casey McGrath,1 and Jolien Creighton,1
1Center for Gravitation, Cosmology & Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee, P.O. Box 413, Milwaukee, WI 53201, USA
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Merging supermassive black hole binaries produce low frequency gravitational
waves which pulsar timing experiments are searching for. Much of the current
theory is developed within the plane-wave formalism, and here we develop the
more general Fresnel formalism. We show that Fresnel corrections to
gravitational wave timing residual models allow novel new measurements to be
made, such as direct measurements of the source distance from the timing
residual phase and frequency, as well as direct measurements of chirp mass
from a monochromatic source. Probing the Fresnel corrections in these models
will require that future pulsar timing arrays include more distant pulsars
across our galaxy, and greatly decrease the uncertainty on these pulsar
distances (to within the order of the gravitational wavelength). However, we
find that with these conditions met the measured source distance uncertainty
can be made less than 10 per cent of the distance to the source for sources
out to $\sim 100$ Mpc, source sky localization can be reduced to sub-arcminute
precision, and source volume localization can be made to less than $1\
\text{Mpc}^{3}$ for sources out to $1$ Gpc distances.

###### keywords:

gravitational waves – quasars: supermassive black holes – pulsars: general –
methods: observational

††pubyear: 2020††pagerange: Fresnel Models for Gravitational Wave Effects on
Pulsar Timing–B

##  1 Introduction

Supermassive black hole binaries (SMBHBs) at the centers of coalescing
galaxies are expected to produce gravitational waves on the order
$\mathcal{O}(1-100\ \text{nHz})$. Currently pulsar timing collaborations such
as the North American Nanohertz Observatory for Gravitational Waves
(NANOGrav), the Parkes Pulsar Timing Array (PPTA), and the European Pulsar
Timing Array (EPTA) have been regularly collecting timing data on numerous
pulsars for over a decade in order to search for the gravitational wave signal
of these sources (Aggarwal et al., 2019; Kerr et al., 2020). Many pulsars have
rotational periods which behave like extremely regular astrophysical clocks,
but in the presence of a continuous gravitational wave these observed periods
will slowly oscillate in time. The gravitational wave induced period shifts in
every rotation of the pulsar then stack constructively over time, causing the
arrival times of the pulsar clock signals to drift in and out of
synchronization with a reference clock. Therefore by timing pulsars across our
galaxy in a pulsar timing array (PTA) and comparing the timing residual data
with our theoretical models for these residuals we can search for the presence
of a gravitational wave and if found, infer the model parameters that describe
the SMBHB source. While the first detection may be the cumulative effect of
gravitational waves coming from many sources all across the sky as a
stochastic background, in this paper we focus on the possibility of detecting
individually loud continuous wave signals coming from SMBHBs.

An important assumption that goes into deriving the gravitational wave timing
residual model is the description of the shape of the gravitational wavefront.
Many studies use models that assume the pulsar timing array receives plane-
waves coming from distant sources (Corbin & Cornish, 2010; Aggarwal et al.,
2019), with the notable exception of Deng & Finn (2011) which considered a
curved wavefront. A second important assumption that is built into these
models is the evolution of the gravitational wave in time. Models can either
treat the SMBHB as a purely monochromatic gravitational wave source, or
include the frequency evolution of the gravitational waves over the thousands
to tens of thousands of years it takes light to travel from the pulsar to the
Earth (Corbin & Cornish, 2010; Aggarwal et al., 2019).

In this paper we organize and classify four regimes of interest which we show
in the left panel of Figure 1, which increase in generality left-to-right and
top-to-bottom and govern the fundamental assumptions built into the timing
residual models. Of these four regimes the plane-wave models are well
established in previous literature. We add a new regime which we label
“Fresnel,” as we will show it becomes important for significant Fresnel
numbers describing the curvature of wavefronts. To this we derive a formula
for the timing residuals in the Fresnel monochromatic regime, and use it to
propose a physically motivated conjecture as to what the most general
gravitational wave residual model should be, the Fresnel frequency evolution
regime.

The inspiration of this work comes primarily from Deng & Finn (2011) and
Corbin & Cornish (2010). Our primary goal is to investigate how well we can
measure the source distance parameter $R$ purely from the Fresnel corrections
due to the curvature of the gravitational wave in our models, like in Deng &
Finn. We adopt a Fisher matrix approach as in Corbin & Cornish to perform this
analysis, and we also include the pulsar distances as free parameters in our
studies following their approach.

For this first study we do emphasize that in order to focus on the measurement
of the source distance and to try and determine what experimental PTA
requirements we will need for this measurement to be feasible, we adopt a
Euclidean cosmology and assume a static universe. This means there is only one
type of distance measurement of the physical distance between the Earth and
source in the work presented here. In a more generalized cosmology, the source
distance parameter that is searched for in timing residual models when working
in the plane-wave regime is the luminosity distance (Aggarwal et al., 2019),
which comes from the amplitude of the timing residual. However, the curvature
corrections from the Fresnel regime introduce a second distance parameter into
the timing residual models, the comoving distance. In this paper we do not
attempt to measure the source distance from the amplitude of the timing
residual, only from the wavefront curvature corrections. Therefore the
distance we recover would represent the comoving distance (not the luminosity
distance) in a more general cosmological framework. The idea of measuring both
distances in a generalized cosmological model was recently explored in
D’Orazio & Loeb (2020) as a means of measuring the Hubble constant using PTA
experiments, and it will also be the follow-up study to this work, using the
formalism and results that we establish in this paper.

It is expected that measuring the comoving distance from the Fresnel
corrections will be more difficult than measuring the luminosity distance from
the amplitude, which will therefore limit the ability to measure cosmological
parameters solely with the gravitational wave. In this study we show that the
source sky and volume localization from measurements of the gravitational wave
source parameters are greatly improved by highly accurate measurements of the
pulsar distances in our PTA. Localizing the source could be invaluable in the
identification of an optical counterpart. If a counterpart is found and the
redshift of its host galaxy determined, then the gravitational waves would
provide a standard siren that would determine the distance of the source and
hence a measurement of the Hubble constant (Holz & Hughes, 2005; Schutz,
1986).

The paper is outlined as follows. In Section 2 we present a new notational
description of gravitational wave timing residuals in pulsars, which we think
helps clarify and emphasize the underlying physics and assumptions that go
into building the timing residual models. Then we develop the actual residual
models in Section 3, derive the Fresnel monochromatic regime formula, and
motivate the full Fresnel frequency evolution regime. In Section 4 we explain
our Fisher matrix method for parameter estimation, and in Section 5 we provide
a fresh look at the pulsar distance wrapping problem and its importance in
parameter estimation. The results of our Fisher analysis are then given in
Section 6, which focus on the novel new insights that the Fresnel regime
offers to pulsar timing. Our final conclusions and future directions of this
work are then given in Section 7.

##  2 Gravitational Wave Timing Residual Theory

For this study we assume a flat and static universe. Additional assumptions
are listed in Appendix A.1. There are eight source parameters for a circular
binary system in the zeroth order Newtonian regime,
$\\{R,\theta,\phi,\iota,\psi,\theta_{0},\mathcal{M},\omega_{0}\\}$. These are
the Earth-source distance (recall we are assuming flat static space), 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}$
(a combination of the individual black hole masses), and orbital frequency
parameter. The pulsar parameters are $\\{L,\theta_{p},\phi_{p}\\}$, which are
the Earth-pulsar distance, and pulsar sky angles.

We choose to use the following notation and conventions in this work. With the
Earth at the center of the coordinate system, the source frame orientation
vectors and the pulsar vector are defined by (see also the right panel of
Figure 1):

$\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}$ (1)

which we use to 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}$
(2)

Writing the polarization tensors in this way groups all of the geometrical
orientation and location angles into the definition of the tensors and keeps
them from mixing the plus and cross metric perturbations, which we find to be
conceptually and mathematically convenient. Another important quantity we will
define here are the following antenna patterns:

$\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}$ (3)

These will be used when writing the residual models in Section 3, and help
decide the detector sensitivity to the gravitational wave based on the Earth-
pulsar-source geometrical alignment.

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

$\displaystyle h^{TT\hat{r}}_{ij}$
$\displaystyle=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],$
(4)
$\displaystyle\mathrm{where}\quad\begin{cases}h_{+}(t)&\equiv-h(t)\cos\big{(}2\Theta(t)\big{)},\\\
h_{\times}(t)&\equiv-h(t)\sin\big{(}2\Theta(t)\big{)},\\\
h(t)&\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega(t)^{2/3},\end{cases}$
(5)

with the angular phase and frequency functions $\Theta(t)$ and $\omega(t)$
defined in Sections 2.1.

A more detailed description of the connection between the gravitational wave
metric perturbation, induced fractional shift of a pulsar’s period $T$, and
its timing residual is offered in Appendix A. The gravitational wave induced
fractional shift of the pulsar’s period $T$ is:

$\frac{\Delta
T}{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{ret}}(t,\vec{x})\big{)}}{\partial
t}\Bigg{\rvert}_{\vec{x}=\vec{x}_{0}(t)}dt,$ (6)

where $t_{\mathrm{obs}}$ is the time a pulsar’s photon is observed arriving at
Earth, $t_{\mathrm{ret}}$ is the retarded time of the gravitational wave
(equation 11) , and $\vec{x}_{0}(t)$ 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; Maggiore,
2018):

$\mathrm{Res}(t)=\int\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\int\frac{T_{\mathrm{obs}}(t_{\mathrm{obs}})-T}{T}dt_{\mathrm{obs}}=\
\mathrm{Obs}(t)-\mathrm{Exp}(t),\\\ $ (7)

###  2.1 Regimes: Monochromatism (A) vs. Frequency Evolution (B)

In the simplest model of the gravitational waves emitted by a binary system,
one neglects the energy loss by emission of gravitational radiation (Maggiore,
2008). Under this assumption, the binary’s energy will remain conserved, the
orbit will not coalesce (it has “infinite” coalescence time), and hence the
orbital period/frequency will remain constant. This assumption means that the
gravitational waves emitted by the binary system will be monochromatic (their
amplitudes and frequencies remain constant for all time).

For the monochromatic gravitational wave model, the orbital frequency and
phase are given by:

$\displaystyle\begin{cases}\omega(t)&=\omega_{0},\\\
\Theta(t)&=\theta_{0}+\omega_{0}(t-t_{0}),\end{cases}$ (8)

where $t_{0}$ here (and below) denotes the fiducial time for the model. In
this work we will choose the fiducial time to be $t_{0}=-\frac{R}{c}$ (note
that $t=0$ would correspond to the present day start of our experiment on
Earth).

More realistically, however, energy in the binary system is lost due to the
emission of gravitational wave radiation. This means that over time the binary
will coalesce. As the binary coalesces, the orbital frequency will increase,
meaning that the gravitational wave frequency will also increase. This
evolution towards higher frequencies is called “frequency chirping.” Due to
the thousands of years it takes light to travel from the pulsar to the Earth,
it is possible that the gravitational waves emitted from our observed source
could evolve appreciably during this time.

The orbital frequency and phase of the frequency evolving gravitational wave
model are given by (Arzoumanian et al., 2014; Maggiore, 2008; Creighton &
Anderson, 2011):

$\displaystyle\begin{cases}\omega(t)&=\omega_{0}\left[1-\frac{t-t_{0}}{\Delta\tau_{c}}\right]^{-3/8},\\\
\Theta(t)&=\theta_{0}+\theta_{c}\left[1-\left(\frac{\omega(t)}{\omega_{0}}\right)^{-5/3}\right].\end{cases}$
(9)

Because a frequency evolving source will eventually coalesce, the model is
governed by the physically significant quantities $\Delta\tau_{c}$ which is
the “coalescence time,” and $\theta_{c}$ which is the “coalescence angle” (the
total angle swept out before the system coalesces). Both of these quantities
are measured from the chosen fiducial time, and are given by:

$\displaystyle\begin{cases}\begin{tabular}[]{c l
l}$\Delta\tau_{c}$&$\equiv\frac{5}{256}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\frac{1}{\omega_{0}^{8/3}}$&$\approx(430\text{
Myr})\left(\frac{10^{9}M_{\odot}}{\mathcal{M}}\right)^{5/3}\left(\frac{1\text{
nHz}}{\omega_{0}}\right)^{8/3}$ ,\\\
$\theta_{c}$&$\equiv\frac{8}{5}\Delta\tau_{c}\omega_{0}\ \ =\ \
\frac{1}{32}\left(\frac{c^{3}}{G\mathcal{M}\omega_{0}}\right)^{5/3}$&$\approx(3.5\times
10^{6}\text{
cycles})\left(\frac{10^{9}M_{\odot}}{\mathcal{M}}\right)^{5/3}\left(\frac{1\text{
nHz}}{\omega_{0}}\right)^{5/3}$ .\end{tabular}\end{cases}$ (10)

Note that in the “large coalescence time” limit
($\Delta\tau_{c}\rightarrow\infty$), the frequency evolution model (equation
9) reduces to the monochromatic model (equation 8) as we would expect.

###  2.2 Regimes: Plane-Wave (I) vs. Fresnel (II)

The wavefront of a traveling wave is traced by the retarded time, as shown in
the right panel of Figure 1. For a static flat universe, the retarded time is:

$t_{\mathrm{ret}}\quad=\quad
t-\frac{|\vec{x}-\vec{x}^{{}^{\prime}}|}{c}\quad=\quad
t\quad-\underbrace{\frac{|\vec{x}^{{}^{\prime}}|}{c}}_{\mathrm{``Far}\
\mathrm{Field"}}\quad+\quad\underbrace{\left(\hat{x}^{{}^{\prime}}\cdot\hat{x}\right)\frac{|\vec{x}|}{c}}_{\mathrm{``Plane-
Wave"}}\quad-\quad\underbrace{\frac{1}{2}\left(1-\left(\hat{x}^{{}^{\prime}}\cdot\hat{x}\right)^{2}\right)\frac{|\vec{x}|}{c}\frac{|\vec{x}|}{|\vec{x}^{{}^{\prime}}|}}_{\mathrm{``Fresnel"}}\quad+\quad\ldots$
(11)

where $\vec{x}^{{}^{\prime}}$ is the wave source’s position and $\vec{x}$ is
the field point of interest, and Earth is the origin. In this Taylor expansion
for $|\vec{x}^{{}^{\prime}}|\gg|\vec{x}|$ we label the first three terms which
are common in wave physics and optics. This gives us a way of categorizing
different wave solution approximation regimes, which will ultimately be the
crux of our work in this paper.

Figure 1: (Left) Our classification of the gravitational wave regimes. They
increase in generality from left-to-right and top-to-bottom. The models reduce
from frequency evolution to monochromatism in the large coalescence time limit
($\Delta\tau_{c}\rightarrow\infty$), and from Fresnel to plane-wave in either
the small Fresnel number limit ($F\rightarrow 0$) or the natural plane-wave
limit ($\hat{r}\cdot\hat{p}\rightarrow\pm 1$). The plane-wave regime has been
well studied in pulsar timing literature, and in this paper we add a new
“Fresnel” regime. (Right) An example contour diagram of the retarded time, in
the source-Earth-pulsar plane at fixed time $t=0$. The approximations indicate
up to what term is included in the expansion in equation 11. Since the lines
show constant $t_{\mathrm{ret}}$, this traces the curvature of the wavefront.
For illustrative purposes that exaggerate and show the differences in the
approximation regimes, we set the pulsar at a distance $L=5$ kpc and the
source at $R=30$ kpc from the Earth. A source with frequency $\omega_{0}=1$
nHz along with this pulsar would have a Fresnel number $F=27.3$ in this
example. As we can see here, the curvature of the wavefront means the wave
arrives at the pulsar’s location later than what is predicted by the plane-
wave approximation.

The Fresnel term gives the first order curvature of the physical wavefront.
Conceptually, it introduces an additional time delay to the arrival time of
the wavefront predicted by the plane-wave approximation (see the right panel
of Figure 1). We want to understand when this curvature term/time delay will
become significant in our timing residual model. For an order-of-magnitude
estimate consider a monochromatic wave. With our source at
$\vec{x}^{{}^{\prime}}=R\hat{r}$ and the pulsar at $\vec{x}=L\hat{p}$, the
timing residual solution (see Section 3) will contain terms that behave like
$\sin\left(2\omega_{0}t_{\mathrm{ret}}\right)$. These terms cycle on the
interval $0\rightarrow 2\pi$, so the Fresnel term in equation 11 will become
significant when it is on the order of a cycle, that is roughly when
$2\pi\sim\omega_{0}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\left(\frac{L}{c}\right)\left(\frac{L}{R}\right)\equiv\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\pi
F$, where:

$\displaystyle F$
$\displaystyle\equiv\frac{L^{2}}{\lambda_{\mathrm{gw}}R}\approx
0.0003\left(\frac{\omega_{0}}{1\text{ nHz}}\right)\left(\frac{L}{1\text{
kpc}}\right)^{2}\left(\frac{100\text{ Mpc}}{R}\right),$ (12)
$\displaystyle\mathrm{and}\quad\lambda_{\mathrm{gw}}\equiv\frac{\pi
c}{\omega_{0}}\approx(30.5\text{ pc})\left(\frac{1\text{
nHz}}{\omega_{0}}\right),$ (13)

which is the standard definition of the Fresnel number. If we ignore the
geometric terms here, then when $F\sim\mathcal{O}(1)$ the Fresnel term in
equation 11 will contribute significantly, suggesting that the simple plane-
wave approximation will be less accurate and that the Fresnel approximation
will be required. Deng & Finn (2011) reached this same conclusion in their
work. Even when there exists significant frequency evolution in the timing
residual model, the definition of the Fresnel number still offers us a
convenient proxy of the importance of wavefront curvature within the model
itself.

There are two important limits to keep in mind here. The “small Fresnel
number” limit ($F\rightarrow 0$) and the “natural plane-wave limit”
($\hat{r}\cdot\hat{p}\rightarrow\pm 1$). In either limit, we expect to recover
the plane-wave result. In the natural plane-wave limit the source, Earth, and
pulsar become perfectly aligned ($\hat{r}\cdot\hat{p}=1$) or anti-aligned
($\hat{r}\cdot\hat{p}=-1$), which means no amount of curvature in the
wavefront would affect the timing any differently than a plane-wave.

##  3 Four Model Regimes

Solving equations 6 and 7 under the various assumptions and approximations in
Sections 2.1 and 2.2 gives four different gravitational wave timing residual
models, which we compare here. Some derivation detail and references are
provided in Appendix A.

Numerous sources use alternative notations and conventions when expressing
these models. For example see Aggarwal et al. (2019); Arzoumanian et al.
(2014); Ellis (2013) for alternative ways of writing the residual model in the
plane-wave frequency evolution regime. One goal of this paper is simply to
present the familiar results in a new light, to make them more transparent and
streamlined. Furthermore, direct comparison of the models between the regimes
in this notation helps provide insight into the physical difference between
each model, and what is actually changing mathematically between them.

###  3.1 Plane-Wave, Monochromatic (IA)

In the plane-wave monochromatic gravitational wave regime, the timing residual
can be expressed compactly as:

$\displaystyle\mathrm{Res}(t)$
$\displaystyle=\frac{F^{\textsc{A}}}{4\omega_{0}}\bigg{[}h_{\textsc{A}}\Big{(}\Theta_{E}-\frac{\pi}{4}\Big{)}-h_{\textsc{A}}\Big{(}\Theta_{P}-\frac{\pi}{4}\Big{)}\bigg{]}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(14) $\displaystyle\underset{\left(t_{0}\ =\
-\frac{R}{c}\right)}{\mathrm{where}}\quad\begin{cases}\Theta_{E}\equiv\Theta\Big{(}t-\frac{R}{c}\Big{)}&\equiv\theta_{0}+\omega_{0}t,\\\
\Theta_{P}\equiv\Theta\Big{(}t-\frac{R}{c}-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\Big{)}&\equiv\theta_{0}+\omega_{0}\Big{(}t-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\Big{)},\end{cases}$
(15) $\displaystyle\qquad\qquad\ \ \ \begin{cases}h_{+}\equiv-
h_{0}\cos\big{(}2\Theta\big{)},\\\ h_{\times}\equiv-
h_{0}\sin\big{(}2\Theta\big{)},\\\
h_{0}\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega_{0}^{2/3},\end{cases}$
(16)

along with equations 1, 2, and 3. The term in brackets is the difference in
the metric perturbation at the Earth and pulsar locations when the pulsar’s
light first left the pulsar, and when it finally arrived at the Earth (denoted
with “E” and “P” subscripts here and below). The residual in the plane-wave
regime only depends on the endpoints of the photon’s motion, an “Earth term”
and a “pulsar term.” Note that in this regime the source chirp mass
$\mathcal{M}$ and distance $R$ parameters are fully degenerate in the
residual. Each parameter only appears in the combination
$\frac{\mathcal{M}^{5/3}}{R}$ with the other (see equation 16), therefore we
would not be able to measure these parameters independently through actual
observations (only this specific combination of the two).

An alternative but equivalent way of writing equation 14 is:

$\displaystyle\mathrm{Res}(t)$ $\displaystyle=A_{c(\mathrm{IA})}\
F^{\textsc{A}}\
s_{E,\textsc{A}}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$ (17)
$\displaystyle\mathrm{where}\quad\begin{cases}A_{c(\mathrm{IA})}\equiv\frac{h_{c(\mathrm{IA})}}{4\omega_{0}},\\\
h_{c(\mathrm{IA})}\equiv
2h_{0}\sin\left(\omega_{0}\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\right),\end{cases}$
(18) $\displaystyle\qquad\quad\ \
\begin{cases}s_{E,+}\equiv\sin\left(2\Theta_{E}-\frac{\pi}{2}-\omega_{0}\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\right),\\\
s_{E,\times}\equiv-\cos\left(2\Theta_{E}-\frac{\pi}{2}-\omega_{0}\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\right).\end{cases}$
(19)

The benefit of this notation is that it allows us to define in equation 18 a
“characteristic strain” $h_{c(\mathrm{IA})}$ and “characteristic timing
residual amplitude” $A_{c(\mathrm{IA})}$ for the IA regime. These
characteristic terms will be useful when comparing this model to the Fresnel
monochromatic model IIA, in Section 3.3.2.

Note that the plane-wave monochromatic model converges to zero in both of the
natural plane-wave limits (when $\hat{r}\cdot\hat{p}=\pm 1$). As the system
becomes anti-aligned the antenna patterns $F^{\textsc{A}}\rightarrow 0$. As
the system becomes aligned the pulsar phase $\Theta_{P}\rightarrow\Theta_{E}$
in equation 15, and the difference in the metric perturbations at the Earth
and pulsar go to zero fast enough to kill the entire timing residual
(conceptually the photons “surf” the gravitational waves in this second case).
In equation 17 this is reflected in the characteristic strain and amplitudes
themselves, which both go to zero in the alignment case.

###  3.2 Plane-Wave, Frequency Evolution (IB)

In the plane-wave frequency evolution gravitational wave regime, the timing
residual can be expressed compactly as:

$\displaystyle\mathrm{Res}(t)$
$\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{(}\omega_{0P},\Theta_{P}-\frac{\pi}{4}\Big{)}}{\omega_{0P}}\right]\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(20) $\displaystyle\underset{\left(t_{0}\ =\
-\frac{R}{c}\right)}{\mathrm{where}}\quad\begin{cases}\omega_{0E}&\equiv\omega_{0},\\\
\omega_{0P}&\equiv\omega_{0}\left[1+\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}}{\Delta\tau_{c}}\right]^{-3/8},\\\
\theta_{0E}&=\theta_{0},\\\
\theta_{0P}&=\theta_{0}+\theta_{c}\left(1-\left[1+\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}}{\Delta\tau_{c}}\right]^{5/8}\right),\\\
\Theta_{E}&=\theta_{0E}+\omega_{0E}t,\\\
\Theta_{P}&=\theta_{0P}+\omega_{0P}t,\end{cases}$ (21)

along with equations 1, 2, 3, 5, and 10. Like in the monochromatic case, this
model only depends on the endpoints of the photon’s motion. However, in the
case that frequency evolution is non-negligible there is a slight difference
in the frequency of the gravitational wave at the pulsar when the pulsar’s
light first left, and the frequency of the gravitational wave at the Earth
when the pulsar’s light finally arrived. Because the Earth and pulsar terms
combine, in the cases where frequency evolution is non-negligible we expect
this to produce “beats” in the residual data that we observe.

Equation 14 is recovered from equation 20 in the large coalescence time limit.
While this regime allows for frequency evolution over the thousands of years
between emission and observation of the pulsar’s light pulses, it still
assumes frequency evolution is negligible over the experiment’s observation
time scale. Notice that in this regime the degeneracy of the source chirp mass
$\mathcal{M}$ and distance $R$ parameters is now broken since $\mathcal{M}$
now appears in the phase and frequency of the pulsar term (within
$\Delta\tau_{c}$ and $\theta_{c}$ of equation 21) as well as the amplitude (in
equation 5). Unlike in the plane-wave monochromatic regime, $\mathcal{M}$ is
not in the same combination with $R$ everywhere. So in principle a source well
described by this model with significant frequency evolution should allow for
independent measurements of both of these source parameters.

###  3.3 Fresnel, Monochromatic (IIA)

Here we provide a new model to the literature which probes a regime we call
the “Fresnel monochromatic” regime. An outline of a derivation for this model
is given in Appendix A. The key to deriving this model is that we keep out to
the Fresnel term in the expansion of the retarded time in equation 11.
Conceptually, we account for the additional time delay of a curved
gravitational wavefront (as compared to a plane wavefront) along the path of a
photon traveling from the pulsar to the Earth, as the example in the right
panel of Figure 1 shows. This will slightly alter the phase and frequency of
the residual in the pulsar term, which we will show later can produce
measurable effects. It is important to realize that this derivation and
formalism does not correct for wavefront curvature effects in the geometrical
orientation terms like the antenna patterns. We expect those corrections to
introduce curvature effects of order $\mathcal{O}\left(\frac{L}{R}\right)$ in
the amplitude of the pulsar term, which we do not expect to be measurable.
This is unlike the order $\mathcal{O}\left(\frac{L}{R}\right)$ in the retarded
time which can lead to $F\sim\mathcal{O}(1)$ as discussed in Section 2.2. Due
to this higher sensitivity to small corrections in the phase and frequency of
the pulsar term, and thanks to how they combine with the Earth term to
potentially further produce beats in the timing residual signal, these
corrections can influence the timing residual model in a significant way.

The model we have derived can be expressed as:

$\displaystyle\mathrm{Res}(t)$
$\displaystyle=\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\left(\frac{1}{8}\frac{\lambda_{\mathrm{gw}}}{c}\frac{R}{c}\frac{1}{\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)}\right)^{1/2}\Bigg{[}\Big{\\{}C\left(\eta_{2}\right)-C\left(\eta_{1}\right)\Big{\\}}h_{\textsc{A}}\left(\Theta^{{}^{\prime}}\right)\
+\
\Big{\\{}S\left(\eta_{2}\right)-S\left(\eta_{1}\right)\Big{\\}}h_{\textsc{A}}\left(\Theta^{{}^{\prime}}-\frac{\pi}{4}\right)\Bigg{]}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(22) $\displaystyle\underset{\left(t_{0}\ =\
-\frac{R}{c}\right)}{\mathrm{where}}\quad\Theta^{{}^{\prime}}\equiv\Theta\left(t+\frac{1}{2}\frac{R}{c}\frac{\left(1-\hat{r}\cdot\hat{p}\right)}{\left(1+\hat{r}\cdot\hat{p}\right)}\right)\equiv\theta_{0}+\omega_{0}\left(t+\frac{1}{2}\frac{R}{c}\frac{\left(1-\hat{r}\cdot\hat{p}\right)}{\left(1+\hat{r}\cdot\hat{p}\right)}\right),$
(23) $\displaystyle\qquad\qquad\ \
\begin{cases}\eta_{1}\equiv\left(2\frac{R}{\lambda_{\mathrm{gw}}}\frac{\left(1-\hat{r}\cdot\hat{p}\right)}{\left(1+\hat{r}\cdot\hat{p}\right)}\right)^{1/2},\\\
\eta_{2}\equiv\eta_{1}\Big{[}1+\left(1+\hat{r}\cdot\hat{p}\right)\left(\frac{L}{R}\right)\Big{]},\end{cases}$
(24)

along with equations 1, 2, 3, 13, and 16. Here $S$ and $C$ are the Fresnel
integrals (see Appendix A). The Fresnel monochromatic model equation 22
approaches the plane-wave monochromatic model equation 14 as we would expect
in the small Fresnel number limit ($F\rightarrow 0$) and the natural plane-
wave limit ($\hat{r}\cdot\hat{p}\rightarrow\pm 1$). As discussed in Section
3.1, both of these effects result in zero residual for these two precise
alignments hence the natural plane-wave limit will not be of much interest to
us going forward in this study.

A very important conclusion here is to notice that in this regime the
degeneracy of the source chirp mass $\mathcal{M}$ and distance $R$ parameters
is also broken, like it was in the plane-wave frequency evolution regime. Here
(and in the models in Sections 3.3.1 and 3.3.2), the distance $R$ now appears
independently from a combination with the chirp mass $\mathcal{M}$. So in
principle a source well described by this model with significant Fresnel
numbers among the pulsars in the PTA should allow for independent measurements
of both of these source parameters.

####  3.3.1 Asymptotic Model

Both $\eta_{2},\ \eta_{1}\gg 1$ unless they are very near the natural plane-
wave limit. When using this model we will normally not be interested in the
natural plane-wave limit, in which case we can replace the Fresnel integrals
in equation 22 with their asymptotic expansions (see equation 64). Doing this
and dropping terms of $\mathcal{O}\left(\frac{L}{R}\right)$ in the amplitude
(not in any phase terms) let’s us replace equation 22 with:

$\mathrm{Res}(t)\approx\frac{F^{\textsc{A}}}{4\omega_{0}}\Bigg{[}\Bigg{\\{}\sin\bigg{(}\frac{\pi}{2}\eta_{2}^{2}\bigg{)}-\sin\bigg{(}\frac{\pi}{2}\eta_{1}^{2}\bigg{)}\Bigg{\\}}h_{\textsc{A}}\left(\Theta^{{}^{\prime}}\right)+\Bigg{\\{}-\cos\bigg{(}\frac{\pi}{2}\eta_{2}^{2}\bigg{)}+\cos\bigg{(}\frac{\pi}{2}\eta_{1}^{2}\bigg{)}\Bigg{\\}}h_{\textsc{A}}\left(\Theta^{{}^{\prime}}-\frac{\pi}{4}\right)\Bigg{]},$
(25)

where now we see the Fresnel formalism solution more closely resembles that of
the plane-wave solution equation 14, as they both share the same amplitude
factor.

####  3.3.2 Heuristic Model

Our study and derivation of the gravitational wave timing residual in both the
plane-wave and Fresnel regimes gives us an insight into the underlying physics
in these models. In the plane-wave regime the timing residual depends only on
an Earth term and a pulsar term. We know that a curved wavefront introduces
additional time delays in the arrival of the wave at these endpoints (see the
right panel of Figure 1). Therefore we begin with the plane-wave formalism,
but modify the retarded time expression at the pulsar endpoint to include the
higher order correction terms in the retarded time which come from a
considering the wavefront curvature. A priori, we define the following
quantity:

$\displaystyle\bar{t}_{\mathrm{ret}}$ $\displaystyle\equiv
t-\frac{L}{c}-\frac{1}{c}\sqrt{R^{2}-2\left(\hat{r}\cdot\hat{p}\right)RL+L^{2}},$
$\displaystyle=t-\frac{R}{c}-\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}{R}+\ldots.$
(26)

This expression is nearly identical to the retarded time equation 11 evaluated
at the source position ($\vec{x}^{{}^{\prime}}=R\hat{r}$) and pulsar position
($\vec{x}=L\hat{p}$), except we are explicitly adding in an additional
$-\frac{L}{c}$ term which does not appear in that equation. However, if we
define this quantity then first notice that the first three terms in the above
expression recover the same retarded time factor we had for the pulsar term in
the plane-wave monochromatic regime equation 15 (which is what we want if this
proposed model is to reduce to the plane-wave monochromatic formalism in the
appropriate limits). Second, notice that the fourth term then adds back in the
Fresnel correction we had previously established in equation 11.

On the basis of this physical interpretation, we therefore write the timing
residual in the Fresnel monochromatic regime as:

$\displaystyle\overline{\mathrm{Res}}(t)$
$\displaystyle\equiv\frac{F^{\textsc{A}}}{4\omega_{0}}\bigg{[}h_{\textsc{A}}\Big{(}\Theta_{E}-\frac{\pi}{4}\Big{)}-h_{\textsc{A}}\Big{(}\bar{\Theta}_{P}-\frac{\pi}{4}\Big{)}\bigg{]}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(27) $\displaystyle\underset{\left(t_{0}\ =\
-\frac{R}{c}\right)}{\mathrm{where}}\quad\begin{cases}\Theta_{E}\equiv\Theta\Big{(}t-\frac{R}{c}\Big{)}&\equiv\theta_{0}+\omega_{0}t,\\\
\bar{\Theta}_{P}\equiv\Theta\big{(}\bar{t}_{\mathrm{ret}}\big{)}&\approx\Theta\Bigg{(}t-\frac{R}{c}-\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}{R}\Bigg{)},\\\
&\equiv\theta_{0}+\omega_{0}\Bigg{(}t-\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}{R}\Bigg{)},\end{cases}$
(28)

along with equations 1, 2, 3, and 16. The overbar notation used here and below
is simply to distinguish that unlike the previous residual models, this has
not been analytically derived but rather proposed.

This proposed model has the required properties that it reduces to the plane-
wave formalism in the small Fresnel number and natural plane-wave limits. It
takes the general form of the plane-wave formalism (equation 14) but with the
corrections motivated by the Fresnel regime study. We tested it numerically
and found it to match the results given by the analytic model equation 22 and
the asymptotic model equation 25 extremely well. Moreover, we found from a
practical standpoint that this formula is easier to study numerically in our
Fisher analysis in Section 4.

Alternatively, we can decompose equation 27 into the sum of a plane-wave part
and a correction part due to the curvature of the gravitational wavefront,
similar to what was done in Deng & Finn (2011):

$\displaystyle\overline{\mathrm{Res}}(t)$
$\displaystyle\equiv\mathrm{Res}_{(\mathrm{IA})}(t)+\mathrm{Res}_{(\mathrm{IIA,correction})}(t),$
(29)
$\displaystyle\mathrm{where}\quad\mathrm{Res}_{(\mathrm{IIA,correction})}(t)\equiv\frac{F^{\textsc{A}}}{4\omega_{0}}\bigg{[}h_{\textsc{A}}\Big{(}\Theta_{P}-\frac{\pi}{4}\Big{)}-h_{\textsc{A}}\Big{(}\bar{\Theta}_{P}-\frac{\pi}{4}\Big{)}\bigg{]},$
$\displaystyle\hskip 98.16191pt=A_{c(\mathrm{IIA,correction})}\
F^{\textsc{A}}\ s_{P,\textsc{A}}$ (30) $\displaystyle\qquad\quad\
\begin{cases}A_{c(\mathrm{IIA,correction})}\equiv\frac{h_{c(\mathrm{IIA,correction})}}{4\omega_{0}},\\\
h_{c(\mathrm{IIA,correction})}\equiv
2h_{0}\sin\left(\frac{\omega_{0}}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{R}\right),\end{cases}$
(31) $\displaystyle\qquad\quad\
\begin{cases}s_{P,+}\equiv\sin\left(2\Theta_{P}-\frac{\pi}{2}-\frac{\omega_{0}}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{R}\right),\\\
s_{P,\times}\equiv-\cos\left(2\Theta_{P}-\frac{\pi}{2}-\frac{\omega_{0}}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{R}\right).\end{cases}$
(32)

Here $\mathrm{Res}_{(\mathrm{IA})}(t)$ is the plane-wave monochromatic model
from Section 3.1, and equation 30 is the correction to the plane-wave regime
due to the wavefront curvature from the Fresnel terms. Note that $\Theta_{P}$
in equations 30 and 32 comes from the IA model equation 15. Just like in
Section 3.1 we can express this correction in terms of its own
“characteristic” strain and timing residual amplitudes
$h_{c(\mathrm{IIA,correction})}$ and $A_{c(\mathrm{IIA,correction})}$. This
provides a very convenient interpretation and comparison of the IA and IIA
regimes. Since the IIA model is identical to the IA model with an additional
correction due to the curvature of the gravitational wavefront, and since we
now have expressions for the characteristic amplitudes of both the IA timing
residual as well as this correction term, we can define the ratio of the
amplitude of the correction to the amplitude of the plane-wave part as:

$\rho\equiv\left|\frac{A_{c(\mathrm{IIA,correction})}}{A_{c(\mathrm{IA})}}\right|=\left|\frac{h_{c(\mathrm{IIA,correction})}}{h_{c(\mathrm{IA})}}\right|=\left|\frac{\sin\left(\frac{\omega_{0}}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{R}\right)}{\sin\left(\omega_{0}\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\right)}\right|.$
(33)

This ratio compares the relative size of the Fresnel correction in the IIA
model to the plane-wave contribution, and is a similar quantity again to what
Deng & Finn calculated in their work. If the value of $\rho\sim\mathcal{O}(1)$
or larger, then the Fresnel terms contribute a very significant correction to
the otherwise plane-wave monochromatic model. Together with the Fresnel number
itself (equation 12), $F$ and $\rho$ are two useful metrics in evaluating the
significance of the curvature corrections to our models.

###  3.4 Fresnel, Frequency Evolution (Conjecture) (IIB)

Given the understanding of the previous timing regimes we now propose an even
more generalized gravitational wave timing residual model, one which includes
the Fresnel corrections, allows for frequency evolution, and still recovers
all of the previous three regimes in the appropriate limits. A full
mathematical derivation along the lines of what was given in Appendix A would
be a useful future project to validate our proposal here. Our conjecture,
however, builds upon the reasoning explained from our insights in Section
3.3.2. Again, this model only corrects the effects of wavefront curvature in
the frequency and phase terms, not in the geometrical orientation terms (as
discussed in Section 3.3).

In the monochromatic formalism we generalized the plane-wave solution by
modifying the retarded time in the pulsar term to include the Fresnel term,
which encoded information about the curvature of the wave. Now we propose
doing the same thing with the frequency evolution formalism. We begin with the
plane-wave frequency evolution regime equation 20, and we modify the pulsar
term’s retarded time using equation 26. This affects both the orbital
frequency of the pulsar term $\bar{\omega}_{0P}$ as well as the phase
$\bar{\theta}_{0P}$. The model can be expressed compactly as:

$\displaystyle\overline{\mathrm{Res}}(t)$
$\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{(}\bar{\omega}_{0P},\bar{\Theta}_{P}-\frac{\pi}{4}\Big{)}}{\bar{\omega}_{0P}}\right]\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(34) $\displaystyle\underset{\left(t_{0}\ =\
-\frac{R}{c}\right)}{\mathrm{where}}\quad\begin{cases}\omega_{0E}&\equiv\omega_{0},\\\
\bar{\omega}_{0P}&\equiv\omega_{0}\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}{R}}{\Delta\tau_{c}}\right]^{-3/8},\\\
\theta_{0E}&=\theta_{0},\\\
\bar{\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}{R}}{\Delta\tau_{c}}\right]^{5/8}\right),\\\
\Theta_{E}&=\theta_{0E}+\omega_{0E}t,\\\
\bar{\Theta}_{P}&=\bar{\theta}_{0P}+\bar{\omega}_{0P}t,\end{cases}$ (35)

along with equations 1, 2, 3, and 16. This model behaves as expected in the
previously established limits, recovering all previous regimes.

And just like we did with the Fresnel monochromatic model in Section 3.3.2, we
can also decompose this formula into a plane-wave part and a correction part
due to the curvature of the gravitational wavefront:

$\displaystyle\overline{\mathrm{Res}}(t)$
$\displaystyle\equiv\mathrm{Res}_{(\mathrm{IB})}(t)+\mathrm{Res}_{(\mathrm{IIB,correction})}(t),$
(36)
$\displaystyle\mathrm{where}\quad\mathrm{Res}_{(\mathrm{IIB,correction})}(t)\equiv\frac{F^{\textsc{A}}}{4}\left[\frac{h_{\textsc{A}}\Big{(}\omega_{0P},\Theta_{P}-\frac{\pi}{4}\Big{)}}{\omega_{0P}}-\frac{h_{\textsc{A}}\Big{(}\bar{\omega}_{0P},\bar{\Theta}_{P}-\frac{\pi}{4}\Big{)}}{\bar{\omega}_{0P}}\right].$
(37)

Here $\mathrm{Res}_{(\mathrm{IB})}(t)$ is the plane-wave frequency evolution
model from Section 3.2, and equation 37 is the correction to the plane-wave
regime due to the wavefront curvature from the Fresnel terms. Note that
$\omega_{0P}$ and $\Theta_{P}$ in equation 37 come from the IB model equation
21.

Finally, like in the Fresnel monochromatic and plane-wave frequency evolution
regimes, there are no degeneracies between any of the source parameters. In
the specific case of the parameters $\mathcal{M}$ and $R$ both appear in
combination in the amplitude terms (equation 16), and both appear
independently in the pulsar term phase and frequency functions (frequency
evolution effects bring in the chirp mass, while Fresnel corrections bring in
the source distance). So in principle a source well described by this model
with significant frequency evolution and Fresnel numbers among the pulsars in
the PTA should allow for independent measurements of all source parameters.

###  3.5 Regime & Model Comparisons

Figure 2 shows four examples that highlight the differences in these regimes.
It is important to keep in mind the limitations and assumptions under which
each of these models are derived when making direct comparisons. The models
can give significantly different and unreliable results if they are being
applied to a regime in which they cannot accurately account for that regime’s
physics. For example, all four models should give the same predictions for a
source with high coalescence time and low Fresnel number since they all can
describe the physics of a plane-wave monochromatic source. However, they
should all predict different timing residuals for a source with low
coalescence time and high Fresnel number since each of the models (apart from
model IIB) have assumptions built into them which prevent them from describing
all of the physics in this scenario (that is IA and IB can’t describe the
effects of a curved wavefront, and IIA can’t describe the effects of frequency
evolution).

The four model regimes can be further subdivided into different categories
given the parameter dependencies of the coalesence time and the Fresnel number
(equations 10 and 12), as shown in Table 1. The parameters
$\left\\{\mathcal{M},\omega_{0}\right\\}$ control the transition between the
monochromatic and frequency evolution regimes, and the parameters
$\left\\{R,\omega_{0},L\right\\}$ control the transition between the plane-
wave and Fresnel regimes (see the left panel of Figure 1). Consequently, the
orbital frequency parameter has the most direct influence on these regimes -
increasing the orbital frequency pushes the models towards both the Fresnel
and frequency evolution regimes, and vice versa. Special consideration should
be taken of this when studying these models, because systems at higher
frequencies will likely be more influenced by both the effects of frequency
evolution and wavefront curvature. This means understanding the the full
Fresnel frequency evolution model is very important, because it has the
potential to predict very different timing residuals than the currently used
plane-wave frequency evolution model.

Table 1: The gravitational wave timing residual model regimes are divided by the coalescence time and Fresnel number limits. These transitions can be expressed in terms of the model parameter limits indicated here. | large coalescence time limit
---
$\left(\Delta\tau_{c}\ \rightarrow\ \infty\right)$
| small Fresnel number limit
---
$\left(F\ \rightarrow\ 0\right)$
| low mass limit | $(\mathcal{M}\rightarrow 0)$
---|---
low orbital frequency limit | $(\omega_{0}\rightarrow 0)$
| far source limit | $(R\rightarrow\infty)$
---|---
low orbital frequency limit | $(\omega_{0}\rightarrow 0)$
small detector limit | $(L\rightarrow 0)$

We also see that the Fresnel number scales as $L^{2}$, which means increasing
the pulsar distance is the easiest way to push the residual model into the
Fresnel regime. This gives us motivation to look for and time pulsars at
greater distances in our PTAs, since it is the one factor that we have more
direct control over in enabling us to explore the Fresnel regime in pulsar
timing experiments. We explore this idea in much greater detail Section 6.

Figure 2: To emphasize regime differences we plot the gravitational wave
timing residuals for a single set of parameters as predicted by each of the
four models. Each panel is labelled by the approximate regime the source is in
based on its physical properties (coalescence time and Fresnel number). In all
cases:
$\left\\{\theta,\phi,\iota,\psi,\theta_{0},\theta_{p},\phi_{p}\right\\}=\left\\{\frac{\pi}{3},0,\frac{\pi}{4},\frac{\pi}{4},0,0,0\right\\}\mathrm{rad}$.
(IA) Low Fresnel number, high coalescence time: As a sanity check, we find
that all four models do correctly converge on the plane-wave monochromatic
regime prediction as expected. Here:
$\left\\{R,\mathcal{M},\omega_{0},L\right\\}=\left\\{50\ \mathrm{Mpc},1\times
10^{8}\ M_{\odot},8\ \mathrm{nHz},0.8\ \mathrm{kpc}\right\\}$. (IB) Low
Fresnel number, low coalescence time: The monochromatic models separate from
the frequency evolution models since they cannot capture the significant
change in the orbital evolution of the source. Additionally, since the Fresnel
number is not significant, the Fresnel models converge on the plane-wave
models. Here: $\left\\{R,\mathcal{M},\omega_{0},L\right\\}=\left\\{100\
\mathrm{Mpc},5\times 10^{8}\ M_{\odot},50\ \mathrm{nHz},0.5\
\mathrm{kpc}\right\\}$. (IIA) High Fresnel number, high coalescence time: The
plane-wave models separate from the Fresnel models since they cannot capture
the significant effect of the wavefront’s curvature. Additionally, since
frequency evolution is not significant, those models converge on the
monochromatic model predictions. Here:
$\left\\{R,\mathcal{M},\omega_{0},L\right\\}=\left\\{50\ \mathrm{Mpc},6\times
10^{7}\ M_{\odot},10\ \mathrm{nHz},7\ \mathrm{kpc}\right\\}$. (IIB) High
Fresnel number, low coalesence time: Here the monochromatic models can no
longer capture the interesting physics such as the beat frequency due to the
frequency evolution of the source, and the plane-wave models can’t reliably
account for the additional effects that wavefront curvature introduces. All
four models begin to predict very different results for the same set of
parameters due to their respective limitations and underlying assumptions, and
the most reliable model becomes the full Fresnel frequency evolution model.
Here: $\left\\{R,\mathcal{M},\omega_{0},L\right\\}=\left\\{100\
\mathrm{Mpc},6\times 10^{8}\ M_{\odot},50\ \mathrm{nHz},10\
\mathrm{kpc}\right\\}$.

###  3.6 Regime Search Parameters & Degeneracies

Each of the four model regimes allow for the measurement of different source
parameters. However certain parameters are degenerate or highly covariant with
others in certain regimes, which has motivated us to parametrize the model in
specific ways. Table 2 indicates the default parameters we chose in each of
our model regimes.

Table 2: Our default source parameters $\vec{s}$ for the timing residual models in each regime, based on the model degeneracies and covariances. The parameter $A_{\mathrm{E,res}}$ is the “Earth term timing residual amplitude” (equation 38). In general, frequency evolution allows the direct measurement of chirp mass $\mathcal{M}$ in the frequency $\omega(t)$ and phase $\Theta(t)$ terms in the models, and Fresnel corrections allow for direct measurement of source distance $R$ from the frequency and phase terms. Fisher matrices with alternative parametrization choices were computed using 49. | Plane-Wave | Fresnel
---|---|---
Monochrome | ${\vec{s}_{\mathrm{IA}}=\left\\{A_{\mathrm{E,res}},\theta,\phi,\iota,\psi,\theta_{0},\omega_{0}\right\\}}$ | $\vec{s}_{\mathrm{IIA}}=\vec{s}_{\mathrm{IA}}\cup\left\\{R\right\\}$
Frequency Evolution | $\vec{s}_{\mathrm{IB}}=\vec{s}_{\mathrm{IA}}\cup\left\\{\mathcal{M}\right\\}$ | $\vec{s}_{\mathrm{IIB}}=\vec{s}_{\mathrm{IA}}\cup\left\\{R,\mathcal{M}\right\\}$

Starting with the most basic model, the plane-wave monochromatic regime, the
source distance $R$ and chirp mass $\mathcal{M}$ are fully degenerate, and
appear only in the amplitude of the residual (equation 14). We use equations
14 and 16 to introduce the following parameter:

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

which is the “Earth term timing residual amplitude,” a similar quantity to
characteristic amplitudes (equations 18 and 31). Therefore in the plane-wave
monochromatic regime, this Earth term timing residual amplitude parameter
replaces the chirp mass and source distance parameters. This new parameter
also removes the orbital frequency $\omega_{0}$ from the timing residual
amplitude terms in our model. Therefore $\omega_{0}$ is only being measured
from its contribution to the time evolution of the phase of the wave (equation
15), not from its contribution to the timing residual amplitude.

As we discussed in the previous sections the $R$-$\mathcal{M}$ degeneracy is
broken in the other three regimes due to either frequency evolution or Fresnel
curvature effects. However, we found that if we removed the
$A_{E,\mathrm{res}}$ parameter from these model and included $R$,
$\mathcal{M}$, and $\omega_{0}$ each separately in the amplitude of the timing
residual, it introduced strong covariances between the three parameters which
made it very difficult to measure each of them independently. This is because
the combination of $\frac{\mathcal{M}}{R\omega_{0}^{1/3}}$ appears in both the
Earth and pulsar terms of every timing residual measured, independent of the
pulsar being timed, which makes these parameters covariant. We need the
differences caused by frequency evolution and Fresnel corrections that are due
to different pulsar distances and locations across the sky to help reduce the
covariances between these three parameters. Therefore, in every regime we
continue to include the Earth term timing residual amplitude
$A_{E,\mathrm{res}}$ as its own separate parameter.111In the two frequency
evolution regimes IB and IIB, we wrote the amplitude of the pulsar term as the
combination of parameters
$A_{\mathrm{P,res}}=A_{\mathrm{E,res}}\left(\frac{\omega_{0}}{\omega_{0P}}\right)^{1/3}$
for the plane-wave regime, and
$A_{\mathrm{P,res}}=A_{\mathrm{E,res}}\left(\frac{\omega_{0}}{\bar{\omega}_{0P}}\right)^{1/3}$
for the Fresnel regime. The other parameters $R$, $\mathcal{M}$, and
$\omega_{0}$ are then measured from the other components of the model where
they appear.

##  4 Fisher Matrix Analysis

To study these models we perform a Fisher matrix analysis to predict how well
pulsar timing experiments will be able to measure and constrain the model
parameters (Wittman, ND; Ly et al., 2017), similar to Corbin & Cornish (2010).
This type of analysis is useful in that it allows us to quickly forecast the
results of different PTA experiments, and look for ways of improving model
parameter measurement through the experimental design of the timing experiment
itself. The inverse of the Fisher matrix gives an estimate of the parameter
covariance matrix $\textbf{{F}}=\textbf{{C}}^{-1}$. Therefore computing a
model’s Fisher matrix then inverting it gives us the covariances of all the
model parameters, which then tells us how well we should be able to measure
these parameters given the experimental design set-up.

Consider our observed gravitational wave timing residual data for every pulsar
at every observation time (indexed by $a$) arranged in a $d$-dimensional
random variable,
$\overrightarrow{\mathrm{Res}}=\left\\{\mathrm{Res}_{a}:a=1,2,\ldots
d\right\\}$, and the model parameters in Table 2 (indexed by $i,j$) arranged
in a $k$-dimensional vector $\vec{X}$. The vector containing the timing
residual random variable data is the sum of some true timing residual due to
the gravitational wave with added random noise, that is
$\overrightarrow{\mathrm{Res}}=\overrightarrow{\mathrm{Res}}_{\mathrm{true}}+\vec{N}$.
The timing data for every pulsar at every observation can then be modeled as a
function of $\vec{X}$ by choosing a timing residual model from Section 3,
organized into the vector $\overrightarrow{\mathrm{Res}}\left(\vec{X}\right)$.
In this notation,
$\overrightarrow{\mathrm{Res}}_{\mathrm{true}}=\overrightarrow{\mathrm{Res}}\left(\vec{X}_{\mathrm{true}}\right)$,
where $\vec{X}_{\mathrm{true}}$ are ultimately the true parameters of the
source that we are trying to measure from the data. For our timing data we
assume that timing observations have no covariances and that all timing
variances are the same, that is
$\mn@boldsymbol{\Sigma}=\mathrm{diag}\left[\sigma^{2}\right]$. Here $\sigma$
is the uncertainty in a timing residual measurement, which in our studies we
typically set to $\sigma=100$ ns for the order of magnitude of present
capabilities (Arzoumanian et al., 2018, 2014; Cordes & Shannon, 2010; Liu et
al., 2011). With this we model the likelihood function as a multivariate
Gaussian in the data $\overrightarrow{\mathrm{Res}}$ as:

$\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}}\left(\vec{X}\right)\bigg{)}^{T}\mn@boldsymbol{\Sigma}^{-1}\bigg{(}\overrightarrow{\mathrm{Res}}-\overrightarrow{\mathrm{Res}}\left(\vec{X}\right)\bigg{)}\right]},$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{d}}\sigma^{d}}\prod_{a}\exp{\left[-\frac{1}{2}\frac{1}{\sigma^{2}}\left(\mathrm{Res}_{a}-\mathrm{Res}_{a}\left(\vec{X}\right)\right)^{2}\right]}.$
(39)

Next we use the definition of the Fisher matrix (Ly et al.):

$\textsf{F}_{ij}\equiv-\Bigg{\langle}\frac{\partial^{2}}{\partial
X_{i}\partial
X_{j}}\ln{\mathcal{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{X}_{\mathrm{true}}\right)}\Bigg{\rangle}=\Bigg{\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)\Bigg{\rangle},$
(40)

where the expectation values are taken over the random variable
$\overrightarrow{\mathrm{Res}}$, to find the matrix elements for our
likelihood model in equation 39:

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

With this an estimate of the posterior distribution of our parameters
$\vec{X}$ is written as a multivariate Gaussian whose covariance matrix is the
inverse of our model’s Fisher matrix:

$p\left(\vec{X}\mid\overrightarrow{\mathrm{Res}}\right)=\sqrt{\frac{\mathrm{det}(\textbf{{F}})}{(2\pi)^{k}}}\exp{\left[-\frac{1}{2}\left(\vec{X}-\vec{X}_{\mathrm{true}}\right)^{T}\textbf{{F}}\left(\vec{X}-\vec{X}_{\mathrm{true}}\right)\right]}.$
(42)

Therefore knowing the Fisher matrix and computing its inverse will give us an
estimate of the parameter covariance matrix, which is our goal.

The main parameters of interest in this model are the source parameters.
However, distances to most pulsars are not known to a high degree of accuracy.
For instance numerous pulsars, including ones timed by NANOGrav, have parallax
uncertainties on the order of 100pc (Arzoumanian et al., 2018; Deller et al.,
2019), which is much larger than the typical gravitational wavelength
$\lambda_{\mathrm{gw}}$ for our sources of interest (see equation 13). The
reason why the size of the uncertainty on the pulsar distances is critical for
our investigation and results is the topic of Section 5. Because of these
uncertainties previous studies have thought to include the pulsar distances as
free parameters in addition to the source parameters in their analyses (Corbin
& Cornish, 2010). This allows us to use the gravitational wave data to also
help measure the distances to pulsars in our PTA.

When including pulsar distances as model parameters, the Fisher matrix takes a
symmetric block-matrix form, separating into a source parameter only symmetric
matrix $\textbf{{F}}^{\textbf{{S}}}$, a pulsar distance parameter only
symmetric matrix $\textbf{{F}}^{\textbf{{L}}}$, and a matrix with cross terms
$\textbf{{F}}^{\textbf{{SL}}}$. Dividing the model parameters into “source”
parameters and “pulsar distance” parameters,
$\vec{X}=\left[\vec{s},\vec{L}\right]$, equation 41 then becomes:

$\displaystyle\begin{cases}\textsf{F}^{\textsf{S}}_{ij}&=\sum\limits_{a}\frac{1}{\sigma^{2}}\left(\frac{\partial\mathrm{Res}_{a}\left(\vec{X}_{\mathrm{true}}\right)}{\partial
s_{i}}\right)\left(\frac{\partial\mathrm{Res}_{a}\left(\vec{X}_{\mathrm{true}}\right)}{\partial
s_{j}}\right)\\\\[8.0pt]
\textsf{F}^{\textsf{SL}}_{ij}&=\sum\limits_{a}\frac{1}{\sigma^{2}}\left(\frac{\partial\mathrm{Res}_{a}\left(\vec{X}_{\mathrm{true}}\right)}{\partial
s_{i}}\right)\left(\frac{\partial\mathrm{Res}_{a}\left(\vec{X}_{\mathrm{true}}\right)}{\partial
L_{j}}\right)\\\\[8.0pt]
\textsf{F}^{\textsf{L}}_{ij}&=\sum\limits_{a}\frac{1}{\sigma^{2}}\left(\frac{\partial\mathrm{Res}_{a}\left(\vec{X}_{\mathrm{true}}\right)}{\partial
L_{i}}\right)^{2}\delta_{ij}\end{cases}\qquad\longleftrightarrow\qquad\textbf{{F}}=\left[\begin{array}[]{c|ccc}\\\\[-5.0pt]
\textbf{{F}}^{\textbf{{S}}}&\cdots&\textbf{{F}}^{\textbf{{SL}}}&\cdots\\\\[8.0pt]
\hline\cr\vdots&\ddots&&0\\\
\left(\textbf{{F}}^{\textbf{{SL}}}\right)^{T}&&\textbf{{F}}^{\textbf{{L}}}&\\\
\vdots&0&&\ddots\end{array}\right].$ (48)

In our work we computed the Fisher matrix and resulting covariance matrix in
terms of the parameters $\vec{s}$ indicated within each regime in Table 2, and
the parameters $\vec{L}=\left[L_{1},L_{2},\ldots\right]$ for each pulsar in
our array. However, sometimes it was useful to compute the Fisher matrix for a
different parametrization choice, $\vec{X}^{{}^{\prime}}$ (for example,
parametrizing the model in the log of a parameter). We computed the new
transformed Fisher matrix through the Jacobian J of our transformation:

$\displaystyle\textsf{F}^{{}^{\prime}}_{ij}=\textsf{J}_{mi}\textsf{J}_{nj}\textsf{F}_{mn}\quad\longleftrightarrow\quad\textbf{{F}}^{\mn@boldsymbol{{}^{\prime}}}=\textbf{{J}}^{T}\textbf{{F}}\textbf{{J}},$
(49) $\displaystyle\mathrm{where}\quad\textsf{J}_{ij}\equiv\frac{\partial
X_{i}}{\partial X^{{}^{\prime}}_{j}}.$

For our Fisher analysis we used the Python SymPy package to first symbolically
write and manipulate our models from Section 3. Inverting the Fisher matrix
was accomplished through singular value decomposition (SVD). We checked that
the condition number of each Fisher matrix was below $10^{14}$ before
performing the inversion (Belsley et al., 1980). This requirement helped
ensure that the SVD procedure accurately calculated the inverse matrix - if
the condition number exceeded this value, then we did not calculate the
covariance matrix. As a final check for accuracy the computed covariance
matrix was multiplied against its original Fisher matrix, and the result
subtracted from the identity matrix. This was checked against an “error
threshold” matrix defined by:

$\textbf{{I}}-\textbf{{F}}\textbf{{C}}<\epsilon\left[\begin{array}[]{cc}1&\cdots\\\
\vdots&\ddots\\\ \end{array}\right],$ (50)

where $\epsilon$ was our “error threshold.” For all results presented in this
paper $\epsilon\leq 0.01$ (and in almost all cases $\leq 10^{-4}$).

##  5 Pulsar Distance Wrapping Problem

A crucial problem in parameter estimation from continuous gravitational wave
pulsar timing residuals is the “pulsar distance wrapping problem” (Ellis,
2013; Corbin & Cornish, 2010). The pulsar distance $L$ affects the phase of
the timing residual terms (see equations 15, 21, 28, and 35). Since the phase
wraps around the interval $[0,2\pi)$, if the pulsar distance is an
unrestricted free parameter, it can become difficult or even impossible to
recover in parameter estimation.222Interestingly, it is also because the phase
wraps around $2\pi$ that we can even try to measure $R$ from the Fresnel
formalism in the first place! Corrections of
$\mathcal{O}\left(\frac{L}{R}\right)$ to the amplitude of the timing residual
model would likely be impossible to measure. However, corrections of
$\mathcal{O}\left(\frac{L}{R}\right)$ that appear in the phase they can become
appreciable over the $2\pi$-interval through the Fresnel number. If the phase
did not cycle on the interval $[0,2\pi)$, then the Fresnel term in the
retarded time equation 11 would remain a much smaller correction to the plane-
wave approximation and not contribute in a significant way to the overall
timing residual. This idea was briefly discussed in Sections 2.2 and 3.3.
Consider for example the two monochromatic regimes IA and IIA. Given their
pulsar term phase dependence $\Theta_{p}$ on $L$ (equations 15 and 28), any
$L^{{}^{\prime}}=L+\Delta L_{n}$ where:

$\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{R}{\left(1+\hat{r}\cdot\hat{p}\right)}+L\right)+\sqrt{\left(\frac{R}{\left(1+\hat{r}\cdot\hat{p}\right)}+L\right)^{2}+n\frac{2\lambda_{\mathrm{gw}}R}{\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)}},&\quad\text{({Heuristic}
Fresnel, monochromatic)}\end{cases}$ (51)

for $n\in\mathbb{Z}$ will give the same timing residual. We will refer to
$\Delta L_{1}$ for the plane-wave monochromatic case as “one wrapping cycle”
in this paper and in the calculations in our results section. Note, assuming
$\hat{r}\cdot\hat{p}\neq\pm 1$, the Fresnel $\Delta L_{n}$ reduces to the
plane-wave $\Delta L_{n}$ in the limit $R\rightarrow\infty$ as expected. The
frequency evolution regimes don’t have an easily defined quantity like this
because their frequencies are time-dependent.

These relations give us a helpful proxy of how small our observational
uncertainty on the pulsar distances $\sigma_{L}$ need to be if we are to
include pulsar distances as free parameters in our models. Namely,
$\sigma_{L}\sim\Delta L_{1}$, a single “wrapping cycle,” which means that the
uncertainty on our pulsar distance measurements needs to be on the order of
the wavelength of the gravitational wave $\lambda_{\mathrm{gw}}$ that we are
trying to measure. For typical sources of interest ($\omega_{0}\sim$ 1 nHz to
100 nHz), the $\lambda_{\mathrm{gw}}$ can range between tens of parsecs to
subparsec distances, which imposes a very strong experimental challenge by
today’s standards.

Note additionally that the geometric scaling factors in equation 51 can either
help or hurt us, depending on if the pulsar happens to be more aligned or
anti-aligned with the source. Considering the plane-wave monochromatic pulsar
wrapping distance $\Delta L_{n}$, if our pulsar happens to be more aligned
with our source (within 90°, $0<\hat{r}\cdot\hat{p}<1$), then the condition
relaxes and we don’t need as tight of an experimental uncertainty on the
pulsar’s distance measurement a priori. But if our pulsar happens to be more
anti-aligned with our source (more than 90°, $-1<\hat{r}\cdot\hat{p}<0$), then
we require even more precise distance measurements (the smallest $\Delta
L_{n}$ being half a gravitational wavelength). In Table 4 we provide for
reference the values of the wrapping cycle $\Delta L_{1}$ proxy (from the
plane-wave monochromatic formula) for each pulsar in our PTA with Source 1
from Table 3.

In practice we found that this problem didn’t pose a significant issue in our
Fisher matrix analyses when using the plane-wave regime models IA and IB, but
did pose a significant issue when using with the Fresnel regime models IIA and
IIB. In the case of the Fresnel regime models, if no restrictions were placed
on the pulsar distance measurements at all then parameter estimation was not
possible because the Fisher matrices could not be accurately inverted (their
condition number exceeding $10^{14}$). Therefore we addressed this problem in
our simulations by adding an uncorrelated Gaussian prior to our knowledge of
the pulsar distances in our Fisher matrix (Wittman, ND). Adding
$\textbf{{F}}^{\textbf{{L}}\mn@boldsymbol{{}^{\prime}}}=\textbf{{F}}^{\textbf{{L}}}+\left(\textbf{{C}}^{\textbf{{L}}}\right)^{-1}$
to the pulsar distance sub-matrix in equation 48, where:

$\textbf{{C}}^{\textbf{{L}}}=\text{diag}\left(\sigma_{L}^{2}\right),$ (52)

represents pulsar distance measurement constraints placed on our experiments,
each pulsar with its own distance uncertainty $\sigma_{L}$. In practice these
uncertainties would likely come from electromagnetic observations of the
pulsars.

In the case of the plane-wave regime models, we did find it possible to
accurately invert the Fisher matrix and recover the source parameters without
the need for this additional pulsar distance prior. Therefore we conclude that
the cause for this is the introduction of the distance parameter $R$ in the
Fresnel regimes which is not in the plane-wave regimes. As we discuss in
detail in Section 6 the parameter $R$ is the most difficult parameter to
measure from our models. This parameter enters the Fresnel models IIA and IIB
in the combination $L/R$, therefore it is highly sensitive to the pulsar
distance wrapping problem and requires the additional pulsar distance prior in
order to constrain the measurement.

##  6 Results

For this study we used the PTA in Table 4, an array of 40 of the pulsars from
NANOGrav’s PTA in Arzoumanian et al. (2018). In our simulations we timed each
of these pulsars for an observation time of 10 years, with a timing cadence of
30 observations/year (unless specified otherwise). Each of these pulsars were
timed at the exact same times, and all of the timings were evenly spaced over
the observation period. Additionally we simulated timing uncertainty $\sigma$
uniformly in all of the timing residuals (see equation 41) and distance
uncertainties $\sigma_{L}$ for each of the pulsars (see equation 52). Finally,
in order to allow use to control the baseline distances to all of the pulsars
in our PTA we introduced a PTA distance “scale factor” term. Therefore a scale
factor of 1 indicates that all of the pulsars in the PTA have their standard
distances given by the $L$ column in Table 4, and at most we increased the
scale factor to 7 (because that would place our farthest pulsar # 40 at
roughly the distance of the Large Magellanic Cloud).

In a real pulsar timing experiment, there are other physical effects that
create timing residuals which must also be included in a complete timing
residual model (Arzoumanian et al., 2018). For example, underestimating the
pulsar’s own period derivative when modeling the time-of-arrivals of the
pulsar’s signal will result in a parabolic timing residual in the fitted data
(Lorimer & Kramer, 2004). In our studies we wanted to focus solely on the
effects of the gravitational waves, without simultaneously modeling other
timing residual sources, to help ensure that the gravitational wave parameters
would not be degenerate with other model parameters.

However, if for example we considered a low frequency gravitational wave
source whose timing residual over our simulated 10 year observation time
didn’t complete a full wave cycle and appeared parabolic (for example, a half
wave cycle), then our parameter estimation results may be overly optimistic.
This is because a real timing experiment would simultaneously fit a quadratic
polynomial and the gravitational wave signal to the timing data (as well as
other models which we are not considering for simplicity here). Therefore we
would be less sensitive to a lower frequency source whose gravitational wave
timing residual predicted by one of our models from Section 3 resembled a
parabola.

So in order to help decide what source frequencies we should not consider in
our Fisher matrix analyses in this paper we first performed the following
calculation. Using following proxy for the SNR of our experiment:

$\mathrm{SNR}=\sqrt{\sum_{a}\left(\frac{\mathrm{Res}_{a}}{\sigma}\right)^{2}},$
(53)

we calculated the SNR of just the gravitational wave timing residuals for
every pulsar at every time, using our IIB model (Section 3.4). Then we
recalculated the SNR, but first subtracted off a best-fit quadratic polynomial
to the original gravitational wave timing residual:

$\mathrm{Res}_{\mathrm{polyfit}}=\mathrm{Res}-\left(a+bt+ct^{2}\right),$ (54)

This polynomial fitting procedure was done for every pulsar in the array,
given its timing data, meaning in a PTA of 40 pulsars there were 120
polynomial best-fit parameters. The result of this is given in Figure 3, which
shows the fractional difference between these two SNR values. If the
fractional difference is low, then we can safely perform our Fisher matrix
analysis on that source given that PTA without needing to model the additional
parabolic timing residual. Therefore we considered only sources in this study
where the fractional difference in the SNR was small, that is with
$\omega_{0}\gtrsim 10$ nHz.

Figure 3: The fractional difference in the gravitational wave timing residual
SNR defined in equation 53, with vs. without a quadratic polynomial best-fit
subtracted from the data, as a function of the source orbital frequency
$\omega_{0}$ and chirp mass $\mathcal{M}$. The fractional difference in SNR is
independent of timing uncertainty $\sigma$, and we also found that it remains
effectively unchanged for different values of $R$ (this specific result was
computed for $R=100$ Mpc). For this study we want to consider sources where
the fractional difference in SNR is negligible, since we are only modelling
gravitational wave effects in our Fisher matrix analyses. Additionally
contours of coalescence time $\Delta\tau_{c}$ (equation 10) are indicated,
with a cut along $\Delta\tau_{c}=1$ kyr due to assumption 14. Contours of
$A_{E,\mathrm{res}}\left(\frac{R}{100\ \mathrm{Mpc}}\right)$ are also
indicated in red, which serve as a useful proxy for how strong the timing
residual signal is (especially when compared to the timing noise $\sigma$).
The fiducial sources from Table 3 are also indicated by their number 1 and 2
here. This was computed using the IIB model (Section 3.4).

Further restrictions on the searchable region of parameter space come from
consideration of the coalescence time (our fundamental assumption 14), as well
as the strength of the timing residual, which as a proxy can be estimated by
the parameter $A_{E,\mathrm{res}}$ (equation 38). These are also indicated in
Figure 3.

In our work we use the coefficient of variation (CV) of a given model
parameter $x$ as a proxy for that parameter’s measureability:

$\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}$ (55)

where $\mu$ and $\sigma$ are the distributions’ parameters (the mean and
standard deviation in the case of a normally distributed $x$). Once we have
the covariance matrix of our parameters from our Fisher analysis, we take the
$\sigma$ values from that matrix and compute the CVs based on how each
respective parameter is distributed (in the case of the regular normal
distribution, $\mu$ is the parameter value we inject into the simulation). All
parameters within a model are normally distributed in a Fisher matrix
approximation (see again equation 42), so we only need the second line of
equation 55 if we parametrize a model parameter $x$ as $\log_{10}(x)$.
Therefore a parameter with a small CV suggests to us that it would be
measurable from the experimental set-up, while a CV of order unity or larger
would suggest an unmeasurable parameter.

###  6.1 Measuring Source Distance

As a point of reference we used Source 1 in Table 3 as a fiducial reference.
This is a nice source because it has a low $\Delta\tau_{c}$ and is therefore
strongly chirping, has a high $A_{E,\mathrm{res}}$ even when the source is at
$R=100$ Mpc (which is good if our timing uncertainty $\sigma\sim 100$ ns), and
a $\lambda_{\mathrm{gw}}\sim 1$ pc which sets the length scale of the typical
wrapping cycle.

####  6.1.1 Figure-of-Merit (F.O.M)

One of the primary goals of this work was to determine how well we could
measure the source distance $R$ from the Fresnel corrections in our pulsar
timing models. As a starting question we asked, “What is the most important
factor in governing the recoverability of the source distance parameter?” The
main take-away we found is that fundamentally, in order to exploit the Fresnel
corrections to our model to measure $R$, our precision is governed by the
tightest we can constrain the pulsar distance uncertainties $\sigma_{L}$ on
the farthest pulsars in our PTA.

To arrive at this conclusion we first simulated two PTA experiments, each with
a selection of 39 “near” pulsars ($L<1$ kpc) and a single “far” pulsar (using
the F.O.M. $L$ values indicated in Table 4). Source 1 was placed at a very
nearby distance of $R=10$ Mpc. In the first case we simulated a prior distance
constraint only on the single far pulsar to within its own wrapping cycle
$\sigma_{L}=\Delta L_{1}$, while all near pulsars were unconstrained. In the
second case we did the opposite, placing no prior constraint on the position
of the single far pulsar and giving all nearby pulsars priors equal to their
wrapping cycles. For both scenarios we built up the size of the PTA by first
simulating the timing array with only the single far pulsar # 40. Then for
every subsequent simulation we added in one additional nearby pulsar, starting
at pulsar # 1. For every simulation we recorded the PTA’s ability to measure
the source distance parameter $R$ by calculating the $\mathrm{CV}_{R}$ value.

Figure 4: Two figure-of-merit studies investigating what is needed in a PTA to
measure the source distance $R$ from Fresnel corrections. See Table 4 for the
pulsars and their F.O.M $L$ values which were used in these calculations.
Source 1 from Table 3 is simulated, at a distance $R=10$ Mpc (therefore
$A_{E,\mathrm{res}}=21\ \mu$s). Note that because this is simply meant to be a
figure-of-merit study to understand the general behavior of changing the
pulsar and PTA parameters, we purposefully simulated a very near and loud
source. In both panels we first simulated the PTA with only the single “far”
pulsar (pulsar # 40), placed at the distances indicated, with a distance
uncertainty constraint of a wrapping cycle $\sigma_{L}=\Delta L_{1}$. For
every subsequent simulation we added in one additional “nearby” pulsar ($L<1$
kpc) with no distance prior constraints, starting with pulsar # 1 through
pulsar # 39. In each simulation we calculated the $\mathrm{CV}_{R}$ value for
that experiment. In all cases, we found that we needed a minimum PTA size of 3
in order to accurately invert our Fisher matrix using SVD. Interestingly
beyond 18 pulsars in the PTA, adding additional pulsars didn’t seem to improve
the recoverability of $R$. These results were computed using the IIB model
(Section 3.4).

Only the first scenario succeeded in producing a reliable measurement of the
source distance - the PTA with one well constrained far pulsar, and many
unconstrained nearby pulsars. Results for this scenario are shown in the top
panel of Figure 4. This experimental set-up shows a marked improvement in the
$\mathrm{CV}_{R}$ value as the distance to the furthest pulsar is increased
(while still holding its $\sigma_{L}$ fixed). However, in the opposite
scenario when all nearby sources were known to within their wrapping cycles
but the single far pulsar was unconstrained, all measurements of
$\mathrm{CV}_{R}$ were well above unity. Therefore even when increasing the
distance of the far pulsar in the array, we found that prior knowledge only on
the nearby pulsars was not enough to allow us to measure the value of the
source distance from the Fresnel corrections.

In summary, even if the distances of many nearby pulsars are known to a great
deal of precision, we cannot recover $R$ unless we can strongly constrain the
distance to the furthest pulsars in the array. In fact that is all we need in
principle - we don’t even need to have distance measurement constraints on the
nearest pulsars in the array a priori. Our conclusion for the reason that
these two F.O.M. studies gave us these results is because nearby pulsars will
likely have negligible Fresnel numbers, meaning that we can’t probe their
Fresnel corrections to measure $R$. Even if these nearby pulsars have very
well constrained distance measurements it is not enough since the actual
magnitude of the Fresnel corrections are too small for these pulsars.
Therefore, in order to actually exploit knowledge coming from the Fresnel
corrections in our models we need to have distant pulsars in the array, which
have sufficiently high Fresnel numbers. But due to the pulsar distance
wrapping problem, if these distant pulsars are not well constrained then we
once again lose the ability to measure the source distance $R$.

We also explored what the effects of increased observation cadence and
improved timing uncertainty $\sigma$ might have on measuring $R$. Using again
a single constrained far pulsar and 39 nearby pulsars, the bottom panel of
Figure 4 shows that if the PTA is small in size a higher timing cadence does
help, but timing resolution makes the largest difference. However, what is
interesting is that regardless of the experimental timing uncertainty or
cadence, all three schemes shown produce the same level of accuracy at
recovering $R$ once our PTA is beyond about 17 pulsars in size. This means a
small PTA with very well timed pulsars is good, but if high timing precision
can’t be achieved then adding more pulsars to the PTA will help just as much.
Similarly cutting a PTA in half and doubling its observation cadence will
produce roughly the same results.

As a final study, since we found that nearby pulsars don’t help with the
measurement of $R$ we therefore investigated the effect of adding additional
constrained far pulsars into the PTA. In Figure 5 we started the initial PTA
with pulsars # 1-30 from Table 4 (again using their F.O.M. $L$ distances), and
then consecutively added the final 10 (each of these with the distances
indicated in Figure 5). Again for this study, nearby pulsars were
unconstrained, while each added far pulsar was simulated with its own wrapping
cycle prior $\sigma_{L}=\Delta L_{1}$. As suspected, we can see that adding
additional well constrained far pulsars improves the measurement of $R$. In a
separate study we also found that if only the first of the far pulsars was
given a wrapping cycle prior constraint, then adding additional unconstrained
far pulsars did not improve the measurement of $R$ notably. So it appears that
the only way to improve and actually probe the source distance from Fresnel
corrections really is to have a PTA with numerous distance pulsars, whose
distances are very well measured.

Figure 5: A figure-of-merit study that shows the improvement of measuring $R$
as additional constrained distant pulsars are added to a PTA. Here the base
PTA consists of nearby pulsars # 1-30 from Table 4 using their F.O.M $L$
values, with no prior knowledge distance constraints. Pulsars # 31-40 are then
added consecutively into the PTA, each with the $L$ distance indicated in the
figure, and each with an uncertainty prior of its own wrapping cycle
$\sigma_{L}=\Delta L_{1}$. We used the same source and source distance $R$ as
in Figure 4. These results were computed using the IIB model (Section 3.4).

####  6.1.2 Distant Sources

Based on our findings, a direct measurement of $R$ given current and near
future standards would be difficult, unless perhaps the source happened to be
very nearby, such as in the Virgo Cluster. But for consideration of the
direction of future technologies, distant sources on the order of $100$ Mpc to
$1$ Gpc would require significant improvements in our PTA. From our figure-of-
merit studies we know that nearby pulsars don’t contribute their Fresnel
corrections towards improving the measurement of $R$ \- for this we need many
distant pulsars in the array with $\sigma_{L}$ uncertainties on the order of
their wrapping distances. Therefore we focused our investigation on what
scenarios would give us good measurements of $R$, and how this might be useful
in the future.

The left panel of Figure 6 shows how increasing the distances of all of the
pulsars with the distance scale factor improves the measurement of $R$. In
effect we will need pulsars across the entire span of our galaxy, and ideally
out to the Large and Small Magellanic Clouds, with distance uncertainties on
the order of $\Delta L_{1}$ or $\lambda_{\mathrm{gw}}$ in order to probe
source distances beyond $100$ Mpc. Since the Fresnel number scales as $L^{2}$,
increasing the scale factor of our PTA yields a marked effect towards
recovering $R$. If, however, the pulsar distance uncertainties also scaled as
$L^{2}$, which we would expect for measurements of pulsar distances made using
parallax, then we found that these effects effectively cancelled each other
out, resulting in no improvement in the recovery of $R$ even with larger PTAs.
Therefore it is crucial that as the pulsar distances in our PTA increases,
that their uncertainties continue to remain of the order of the wrapping cycle
(a quantity which is independent of the pulsar’s distance).

Sources with certain physical properties will also be favored when measuring
Fresnel corrections, as seen in the right panel of Figure 6. In terms of
intrinsic parameters, strongly chirping sources with coalescence times
$\Delta\tau_{c}\sim\mathrm{kyr}-\mathrm{Myr}$ will yield the best measurements
of $R$ through Fresnel corrections. The bias we see in this figure leans
towards lower frequency sources, largely because these sources will have
larger wrapping cycles and therefore require less pulsar distance precision
than higher frequency sources (see equations 51 and 13). Therefore keeping
$\sigma_{L}=1$ pc pinned means the lower frequency sources here benefit from a
higher degree of precision than the higher frequency sources.

Figure 6: Source distance measurement as a function of distance and as a
function of source intrinsic parameters $\mathcal{M}$ and $\omega_{0}$. In
both panels the timing uncertainty $\sigma=100$ ns, the pulsar distance
uncertainty $\sigma_{L}=1$ pc. The red circle indicates corresponding points
between the two graphs, and the red square corresponds to Figure 7. Note, a
scale factor of $7$ puts the furthest pulsar in our PTA at roughly the
distance of the Large Magellanic Cloud. These results were computed using the
IIB model (Section 3.4). (Left) Here we measure $R$ from Source 1 in Table 3.
Note that $\sigma_{L}=1\ \mathrm{pc}\ \sim\lambda_{\mathrm{gw}}$ for this
source. (Right) These sources have the same angular parameters as our fiducial
Source 1 (indicated in Table 3) and are set at $R=100$ Mpc, and the PTA has a
scale factor of $7$. Bias leans towards lower frequencies, where using a fixed
$\sigma_{L}$ value will benefit sources with larger wrapping cycles as
compared to smaller ones.

While our focus so far has been on the recovery of the distance parameter, we
also looked at the measurement of all source parameters and their covariances.
An example of this is shown in Figure 7, where Source 1 was placed at $R=500$
Mpc. Again this is a loud source and $\mathrm{CV}_{A_{E,\mathrm{res}}}=0.0097$
for our simulated $\sigma=100$ ns timing uncertainty. In this case
$\mathrm{CV}_{R}=0.37$, while all remaining parameters had their respective
values of $\mathrm{CV}\leq 0.018$. Sky angles tend to be measured well for the
reasons discussed in Section 6.2. Plus as this is a highly chirping source,
chirp mass and frequency are well measured.

Figure 7: Fisher matrix measurements of the parameters for Source 1, which has
been placed at $R=500$ Mpc (therefore $A_{E,\mathrm{res}}=414$ ns). Here the
PTA has a scale factor of $7$ and all pulsars have been given $\sigma_{L}=1\
\mathrm{pc}\ \sim\lambda_{\mathrm{gw}}$, with timing uncertainty $\sigma=100$
ns. This simulation and experimental set-up corresponds to the red square
indicated in Figure 6. The Fisher matrix analysis predicts that this source
can be recovered with: $\left\\{\mathrm{CV}_{A_{E,\mathrm{res}}},\
\mathrm{CV}_{R},\ \mathrm{CV}_{\theta},\ \mathrm{CV}_{\phi},\
\mathrm{CV}_{\iota},\ \mathrm{CV}_{\psi},\ \mathrm{CV}_{\theta_{0}},\
\mathrm{CV}_{\mathcal{M}},\ \mathrm{CV}_{\omega_{0}}\right\\}=\left\\{0.0097,\
0.367,\ 3\times 10^{-6},\ 10^{-6},\ 0.0148,\ 0.0175,\ 0.0182,\ 3.6\times
10^{-5},\ 5.3\times 10^{-5}\right\\}$. These results were computed using the
IIB model (Section 3.4).

For this particular example, we found that improving the timing uncertainty by
an order of magnitude to $\sigma=10$ ns improved the indicated $\mathrm{CV}$
values by one order of magnitude for the parameters
$\left\\{A_{E,\mathrm{res}},\iota,\psi,\theta_{0}\right\\}$, while all
remaining parameter $\mathrm{CV}$ values remained the same. In general we
found that these four parameters were most sensitive to timing accuracy, while
mostly insensitive to pulsar distance accuracy. While the measurement of $R$
is the opposite of this as we have already discussed, we found that the
remaining four source parameters
$\left\\{\theta,\phi,\mathcal{M},\omega_{0}\right\\}$ could benefit from both
improvements to timing as well as pulsar distance uncertainties.

###  6.2 Source Localization

As we investigated in the previous section, high precision measurements of
pulsar distances (to within $\sigma_{L}\sim\Delta L_{1}$) allow for the direct
measurement of the source distance $R$ purely through Fresnel corrections. In
this section, we show that this level of precision can also improve
localization of the source on the sky. As we will explain, this improved sky
localization mostly comes from phase and frequency parallax terms which can be
exploited when the pulsar distance measurements are well constrained. When we
combine the sky localization with the recovered distance information, there
exists great potential for future PTAs to pinpoint the galaxy source of the
gravitational waves, which could be of great benefit to astronomers seeking
electromagnetic counterparts to these coalescing binary sources.

To measure the solid angle on the sky that our uncertainty in the measurement
of the source angles $\theta$ and $\phi$ sweeps out, we first compute the
confidence ellipse area for those two parameters $\Delta
A=\pi\chi^{2}\sigma_{\theta^{{}^{\prime}}}\sigma_{\phi^{{}^{\prime}}}$ as
outlined in (Coe, 2009)333Note that in Coe (2009), equation 7 has a typo in
that it is missing the factor of $\chi^{2}$ which we have corrected here in
equation 56.. Here $\sigma_{\theta^{{}^{\prime}}}$ and
$\sigma_{\phi^{{}^{\prime}}}$ are the measured uncertainties along the
principle axes of the ellipse itself, and this quantity represents an area in
the $\theta$-$\phi$ parameter space. On the sky, we can define a small solid
angle as $\Delta\Omega\approx\sin(\theta)\Delta\theta\Delta\phi$, so to
connect these two quantities we measure the solid angle by multiplying our
confidence ellipse by $\sin(\theta)$ to get:

$\displaystyle\Delta\Omega$
$\displaystyle=\pi\chi^{2}\sin(\theta)\sigma_{\theta^{{}^{\prime}}}\sigma_{\phi^{{}^{\prime}}},$
$\displaystyle=\pi\chi^{2}\sin(\theta)\sigma_{\theta}\sigma_{\phi}\sqrt{1-\rho^{2}}.$
(56)

The uncertainties $\sigma_{\theta}$ and $\sigma_{\phi}$ are measured from the
inverse of the Fisher matrix, and $\rho$ is the correlation coefficient of the
parameters $\theta$ and $\phi$. In all of our results here we set
$\chi^{2}=2.279$, which gives the approximate $68$ per cent likelihood area
for our source.

In a similar calculation we also compute the volume uncertainty region for our
source as:

$\displaystyle\Delta V$
$\displaystyle=\frac{4\pi}{3}\chi^{3}\sigma_{x^{{}^{\prime}}}\sigma_{y^{{}^{\prime}}}\sigma_{z^{{}^{\prime}}},$
$\displaystyle\approx\frac{4\pi}{3}\chi^{3}R^{2}\sin(\theta)\sigma_{R^{{}^{\prime}}}\sigma_{\theta^{{}^{\prime}}}\sigma_{\phi^{{}^{\prime}}},$
(57)

where again, primes on the uncertainty variables denote they are the
uncertainties measured along the principle axes of the ellipsoid. We found
that the Fisher matrix parametrized in cartesian coordinates often had a
condition number too large and therefore made the matrix inversion unreliable
and inaccurate. But spherical coordinates were in general a better choice of
parametrization and produced Fisher matrices which could be reliably inverted,
thus we use the approximation for small volumes in the second line of equation
57 to transform to spherical coordinates. In 3D space, the approximate $68$
per cent likelihood volume of our source corresponds to a $\chi^{2}=3.5059$,
which is what we used for all of our results shown here. To calculate the
actual values of the uncorrelated uncertainties along the principle axes
$\sigma_{R^{{}^{\prime}}}$, $\sigma_{\theta^{{}^{\prime}}}$, and
$\sigma_{\phi^{{}^{\prime}}}$ we first took just the 3D sub-matrix of our
inverted Fisher matrix (for parameters $R$, $\theta$, and $\phi$). Next we
applied SVD once again to that matrix since this mathematically decomposes the
matrix into three matrices, one of which is a rectangular diagonal matrix that
contains the singular values of our original matrix, that is the uncorrelated
uncertainties along the principle axes.

First we consider how the sky angles $\theta$ and $\phi$ are being measured
within our models. In all of our timing regimes, part of the measurement comes
from the antenna patterns contained in the amplitudes of the Earth and pulsars
terms of the timing residuals (see equation 3). Additionally, all timing
regimes also have the “phase parallax” term
$\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}$ which contains information of
the sky angles, and appears in the phase $\Theta(t)$ of the pulsar terms of
each model (see equations 15, 21, 28, and 35). Furthermore, in the frequency
evolution regimes IB and IIB this becomes a “phase-frequency parallax” term,
because not only does it appear in the pulsar term phase $\Theta(t)$ but also
in the pulsar term frequency $\omega(t)$. Lastly, the Fresnel regimes add an
additional small order correction
$\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{R}$
to these phase/phase-frequency parallax terms.

Therefore with this in mind we looked at examples of improvements to three PTA
qualities which can directly improve the measurement of the sky localization:
timing precision $\sigma$, pulsar distance precision $\sigma_{L}$, and PTA
distance scale factor. Timing precision helps to measure the sky angles from
the overall amplitude of the timing residual. Due to the $\frac{L}{c}$ factor
in the phase/phase-frequency parallax terms (both with and without the Fresnel
correction), in order to measure information of the sky angles from these
terms we need high precision pulsar distance measurements. Otherwise without
significantly constrained distances the pulsar distance wrapping problem makes
it impossible to probe this term directly. In fact it is partially for this
reason that studies such as Aggarwal et al. (2019) choose to group these
parallax terms into separate phase parameters for every individual pulsar in
their models (which is a similar but alternative approach to introducing the
pulsar distances as free parameters in the model, as discussed in Sections 4
and 5). This sacrifices the additional sky localization information contained
in these parallax terms, but also helps to avoid the wrapping problem.
Finally, the PTA scale factor can help improve the sky angle measurements in
two ways. Creating a longer baseline between the Earth and pulsar both
amplifies the chirping effects due to the greater disparity between the
Earth/pulsar frequencies experienced by the pulsar photons, and boosts the
size of the Fresnel corrections.

Examples of improvements in these three PTA qualities (timing precision,
pulsar distance precision, and scale factor) as applied to Source 1 are shown
in the three panels of Figure 8, respectively. It is very difficult to
disentangle the many different factors that are all competing to change the
measurement of the source’s sky localization (including the PTA qualities just
discussed, as well as intrinsic source parameters which we discuss below).
Therefore we emphasize here that Figure 8 is only meant to serve as one
specific illustration of the general comments that we have made so far, for
the PTA set-ups indicated in each panel and for one particular source.

Figure 8: Various effects of PTA improvements on source sky localization,
shown here for the example of Source 1. The bottom left corner of the first
panel is excluded because our SVD matrix inversion failed our condition number
requirement (see Section 4), therefore the Fisher matrix for these cases
couldn’t be stability inverted. Nevertheless, we can clearly see the trend is
that the overall sky accuracy improves as the source gets closer to us and the
value of $\sigma$ decreases, as we would expect. Compared to the other PTA
parameters, timing precision improvements follow the power law relationship
$\Delta\Omega\propto\sigma R$. Furthermore, timing precision controls the
largest range of magnitude differences in $\Delta\Omega$ as a function of
source distance, followed by pulsar distance precision, then finally PTA scale
size. These results were computed using the IIB model (Section 3.4).

We find that unlike pulsar distance precision and PTA scale factor, timing
precision follows a nice power law relationship, specifically:
$\Delta\Omega\propto\sigma R$. In the first panel the pulsar distance
precision $\sigma_{L}\approx 100\lambda_{\mathrm{gw}}$ and the scalefactor is
1, therefore we expect that most of the knowledge of our sky angles is coming
from the overall amplitude of the timing residual, and not from the phase-
frequency parallax or Fresnel corrections. The second and third panels show
the additional benefits that can be gained when probing the phase-frequency
parallax and Fresnel corrections. Recall that $\lambda_{\mathrm{gw}}\sim 1$ pc
for Source 1. In panel two the Fresnel corrections are likely not adding a
great deal of information towards the sky angles themselves since most of
these values are for $\sigma_{L}>\lambda_{\mathrm{gw}}$. However we still do
see that information can be gained through the phase-frequency parallax as the
precision of the distance measurements to our pulsars improves. As the value
of $\sigma_{L}$ decreases, we see for a given $\Delta\Omega$ resolution there
are turnover points where even small improvements to the pulsar distance
measurements can allow the same source to be localized with that precision at
much greater distances. Finally, with pulsar distance measurements on the
order of the source’s gravitational wavelength, panel three shows that
increasing the baseline distance between pulsars further improves our ability
to measure the source’s sky angles. All in all we see here that timing
precision controls the largest range of magnitude differences in
$\Delta\Omega$ as a function of source distance, followed by pulsar distance
precision, then finally PTA scale size.

The localization of the source also depends greatly on the source’s own
intrinsic parameters $\mathcal{M}$ and $\omega_{0}$. Figure 9 shows that for
the same distance (here $R=100$ Mpc) and PTA set-up, different sources are
localized to $\sim\mathcal{O}\left(0.1\ \mathrm{deg}^{2}\ -\ 10\
\mathrm{arcsec}^{2}\right)$. Notice that this PTA set-up matches that in the
right panel of Figure 6 so we know that the Fresnel corrections are being
probed for the sources where $\mathrm{CV}_{R}<1$, since the distance parameter
$R$ is measured entirely from these corrections in our model. (Also note that
once again, this plot does place a bias on the lower frequency sources, for
the same reason discussed earlier in Section 6.1.2).

Figure 9: Source sky localization, as a function of source frequency and chirp
mass. These are for sources at $R=100$ Mpc, with the same angular parameters
as our fiducial sources (indicated in Table 3), timing uncertainty
$\sigma=100$ ns, pulsar distance uncertainty $\sigma_{L}=1$ pc, and a PTA
scale factor of $7$. Bias leans towards lower frequencies, where using a fixed
$\sigma_{L}$ value will benefit sources with larger wrapping cycles as
compared to smaller ones. These results were computed using the IIB model
(Section 3.4). When we re-calculated this plot using the IB model (Section
3.2) and IA model (Section 3.1) we found that the relative difference in
$\Delta\Omega$ between these two models and the original fully general IIB
model was about $0.1$ and $0.7$, respectively.

In general we conclude that the Fresnel corrections themselves don’t
significantly improve sky localization of a source. Figure 9 shows
$\Delta\Omega$ predicted from the Fisher matrix analysis using the fully
generalized IIB timing residual model (Section 3.4). We repeated the same
calculations of the sky angles for these sources using the IB model which
removes the Fresnel corrections (Section 3.2), and the IA model which removes
both the Fresnel corrections and the frequency evolution effects (Section
3.1), in order to compare the predictions between models. It is important to
note that this isn’t exactly an apples-to-apples comparison, because
fundamentally these three models will compute different timing residuals for
the same source (see again Figure 2 \- the level of difference will depend on
what regime the source is actually in). And the Fisher matrix analysis only
asks how well can we recover the parameters from a given model. Nevertheless,
comparing the results of these sets of calculations still gives us an idea of
how different the sky localization is with vs. without the Fresnel
corrections. What we found was that the magnitude of the relative difference
between $\Delta\Omega$ computed using the models IIB vs. IB and IIB vs. IA was
at most about $0.1$ and $0.7$ for what is shown in Figure 9, respectively.
This suggests that we lose some sky localization precision when we leave out
Fresnel corrections, but we lose more precision when we don’t account for
frequency evolution. However, even with these relative differences we still
find that in the IB and IA models (which don’t include Fresnel corrections)
sky localization ranges from $\sim\mathcal{O}\left(0.1\ \mathrm{deg}^{2}\ -\
10\ \mathrm{arcsec}^{2}\right)$ depending on the intrinsic source parameters.
This along with the earlier observations and statements about Figure 8 is what
leads us to believe sky localization itself is not greatly improved by the
inclusion of Fresnel effects.

As an additional comment, in their separate studies Corbin & Cornish (2010)
reported sky localization could be measured with $\sim\mathcal{O}\left(1-10\
\mathrm{deg}^{2}\right)$ while Deng & Finn (2011) reported sky localization
could be measured with $\sim\mathcal{O}\left(100\ \mathrm{arcsec}^{2}\right)$,
which is a significant improvement. The results discussed here from our study
help to explain why. These studies use different physical models (with and
without Fresnel corrections), different example sources and PTAs, and simulate
very different pulsar distance uncertainties. Corbin & Cornish used values of
$\sigma_{L}\sim\mathcal{O}\left(1-100\ \mathrm{pc}\right)$, while Deng & Finn
used $\sigma_{L}\sim\mathcal{O}\left(0.001-0.01\ \mathrm{pc}\right)$. Here we
see in Figure 9 approximately five orders of magnitude difference in sky
localization depending on the type of source, and in the middle panel of
Figure 8 another five orders of magnitude difference between $\sigma_{L}=1\ -\
100$ pc. Therefore we easily found scenarios where $\Delta\Omega$ ranged from
$\mathrm{arcsec}^{2}$ to $>10\ \mathrm{deg}^{2}$ precision, simply by
accounting for these different effects, PTA qualities, and sources.

While Fresnel corrections in our models may not give much added precision to
the sky location of our source, we can combine the distance measured from
Fresnel corrections with this sky location, to pinpoint our source to within
an uncertainty volume of space. Figure 10 shows how well we could recover
Source 1 as a function of distance and sky position. The results have only
been given where all points have a corresponding $\mathrm{CV}<1$ for the
parameters $R$, $\theta$, and $\phi$. Since the source distance $R$ is the
hardest parameter to measure in general, our ability to localize the source to
within some uncertainty volume in space is most directly determined by how
well we can measure that parameter. Some parts of the sky are less sensitive
than others, and even if we lose the sensitivity within our PTA required to
measure the volume localization $\Delta V$, we can still typically measure a
sky location $\Delta\Omega$ with striking precision. Even at $1$ Gpc, Source 1
was measured here with sub-square arcminute precision, no matter where it was
located on the sky, shown in the right panel of Figure 10.

Figure 10: Source sky and volume and localization calculated for Source 1,
with timing uncertainty $\sigma=100$ ns, and pulsar distance uncertainty
$\sigma_{L}=1\ \mathrm{pc}\ \sim\lambda_{\mathrm{gw}}$. The areas and volumes
represent the $68$ per cent likelihood regions of our source. These results
were computed using the IIB model (Section 3.4). (Right) The top panel shows
volume localization of the source at $R=100$ Mpc (therefore
$A_{E,\mathrm{res}}=2\ \mu$s). In the bottom panel the source is at $R=1$ Gpc
(therefore $A_{E,\mathrm{res}}=207$ ns), and although we can no longer measure
its distance with this PTA (and hence calculate the volume localization), we
can still measure the sky angles from the chirping in the signal to localize
it on the sky. In both the top and bottom panel the PTA has a scalefactor of
$5$.

###  6.3 Measuring Chirp Mass from a Monochromatic Source

Another important novelty of the Fresnel regime is that it allows a
measurement of the source’s chirp mass even if the source is producing
monochromatic gravitational waves. For a monochromatic source, Fresnel
corrections break the $R$-$\mathcal{M}$ degeneracy in the IIA regime (see
Section 3.3). As long as we can measure $A_{E,\mathrm{res}}$, $R$, and
$\omega_{0}$, then from equation 38 we can infer the the chirp mass.

As pointed out in Section 6 and shown in Figure 3, within the limitations of
this study low frequency sources below $\omega_{0}<10$ nHz require a greater
analysis of the impact of other timing residual sources before confidently
measuring the gravitational wave parameters and their uncertainties. Therefore
what we show here is simply meant to be a proof-of-concept.

For this we use the fiducial Source 2 in Table 3. Again considering Figure 3,
an ideal monochromatic source would have infinite $\Delta\tau_{c}$, so we
chose a low chirp mass source on the edge of our SNR limitation. At $R=100$
Mpc, the amplitude of such as source would be very small, so we also simulate
a future timing experiment capable of $\sigma\sim 1$ ns timing uncertainties.
In a future study, it would be interesting to simulate the measurement of a
higher chirp mass system with an orbital frequency $\omega_{0}\sim 1$ nHz, as
this type of source would produce a stronger amplitude.

The left panel of Figure 11 shows an example of the uncertainty propagation to
our measurement of the system’s chirp mass. We found that measuring
$\mathcal{M}$ was most strongly correlated to measuring $A_{E,\mathrm{res}}$
and $R$. Since measuring $A_{E,\mathrm{res}}$ depends on the timing
uncertainty $\sigma$, and measuring $R$ depends on the pulsar distance
uncertainty $\sigma_{L}$, this means different experimental set-ups may depend
more on timing or parallax for the recovery of the system chirp mass. As an
example, this can be seen in the right panel of Figure 11. Initially,
increasing the PTA scale factor notably improves the recovery of
$\mathcal{M}$, since the higher Fresnel numbers and improved parallax
measurements allow for better recovery of $R$ (similar to what was shown
earlier in the left panel of Figure 6). However, around a scale factor of $5$
we reach the threshold where timing precision dominates over parallax in our
ability to recover $\mathcal{M}$. Moreover for this source we found that
increasing $\sigma>1$ ns more dramatically made the recovery of $\mathcal{M}$
dependent on the recovery of $A_{E,\mathrm{res}}$.

Figure 11: Error propagation to the measurement of the source’s chirp mass
$\mathcal{M}$ for a monochromatic source recovered from the measurements of
$A_{E,\mathrm{res}}$, $R$, and $\omega_{0}$, calculated here in both panels
for Source 2. Here all pulsars have been given timing uncertainty $\sigma=1$
ns, and pulsar distance uncertainty $\sigma_{L}=1.5\ \mathrm{pc}\
\sim\lambda_{\mathrm{gw}}$. Note that in this study we chose to parametrize
the model in the log of the parameters $\left\\{A_{E,\mathrm{res}},\ R,\
\omega_{0}\right\\}$. These results were computed using the IIA model (Section
3.3.2). (Left) Here we show the recovery and correlations of the log of the
parameters $\left\\{A_{E,\mathrm{res}},\ R,\ \omega_{0}\right\\}$ in the
triangle plot, and in the upper right corner we show the propagated
uncertainty on the measurement of the log of $\mathcal{M}$. Source 2 has been
placed at $R=100$ Mpc (therefore $A_{E,\mathrm{res}}=0.15$ ns), and here the
PTA has a scale factor of $5$. In this particular example, we find that the
uncertainty on the chirp mass is dominated by the pulsar distance uncertainty
$\sigma_{L}$, which most strongly affects the recovery of $R$. The Fisher
matrix analysis predicts that this source can be recovered with:
$\left\\{\mathrm{CV}_{A_{E,\mathrm{res}}},\ \mathrm{CV}_{R},\
\mathrm{CV}_{\theta},\ \mathrm{CV}_{\phi},\ \mathrm{CV}_{\iota},\
\mathrm{CV}_{\psi},\ \mathrm{CV}_{\theta_{0}},\
\mathrm{CV}_{\omega_{0}}\right\\}=\left\\{0.32,\ 0.23,\ 7\times 10^{-6},\
4\times 10^{-6},\ 0.48,\ 0.58,\ 0.59,5.4\times 10^{-5}\right\\}$. Uncertainty
propagation gives an inferred $\mathrm{CV}_{\mathcal{M}}=0.23$ . (Right) For
small PTA scale factors, recovery of $\mathcal{M}$ is dominated by pulsar
distance uncertainty $\sigma_{L}$. As the scale factor increases, the larger
Fresnel numbers make the measurement of $R$ more precise, and as a result for
scale factors above about $5$ the measurements of $\mathcal{M}$ don’t improve
because they become dominated by timing uncertainty $\sigma$. Note a scale
factor of $7$ puts the furthest pulsar in our PTA at roughly the distance of
the Large Magellanic Cloud.

##  7 Conclusions

In this work we motivate the importance of Fresnel corrections in the effects
of gravitational waves on pulsar timing residual models, and develop what we
call the “Fresnel” timing regime, separate from the previously established
plane-wave regimes which are currently used in pulsar timing searches for
continuous waves (Aggarwal et al., 2019). We derive the Fresnel monochromatic
model analytically, study it asymptotically, and provide a more simple and
heuristically motivated model which is easier to implement in timing
simulations and analyses. With this we build a well motivated conjecture of
the Fresnel frequency evolution timing residual model, which fully generalizes
the previously studied plane-wave frequency evolution models with Fresnel
order corrections to the phase and frequency of the timing residual.
Surrounding this we also offer a framework for understanding the relevant
limits of each model, namely by considering the coalescence time of the source
and the size of its Fresnel number with the pulsars in the timing array.

Then we perform a Fisher matrix analysis for parameter estimation within these
Fresnel regimes, using a representative NANOGrav-related PTA as the base for
the experimental design within our simulations. This was meant to build upon
and bridge the gap between the studies by Corbin & Cornish (2010) and Deng &
Finn (2011). The main challenge in measuring the source distance purely from
Fresnel corrections is the pulsar distance wrapping problem. We offer a new
look at this problem and show how well measured the pulsar distances will need
to be in order to accurately recover the source distance from our searches.
Figure-of-merit studies indicate that the source distance can only be measured
from the timing residuals of distant sources in our PTAs, whose Fresnel
corrections are large enough to be probed. In general distance constraints are
needed on the order of the wrapping cycle on the most distant pulsars in the
PTA, but not on the nearby pulsars whose Fresnel corrections are negligible.

One way to increase the Fresnel corrections in the timing residuals is by
including pulsars at farther distances, since the Fresnel number scales as
$F\propto L^{2}$ (see equation 12). Intrinsic source properties $\mathcal{M}$
and $\omega_{0}$ also play an important role as they control both the amount
of frequency evolution in the residuals and the characteristic wrapping cycle.
In our simulations we find that future PTAs with distant pulsars across the
Milky Way galaxy whose distances are constrained to the order of the wrapping
cycle can measure the source distance with a $\mathrm{CV}_{R}\la 0.1$ for
sources at out to $\mathcal{O}(100\ \text{Gpc})$ distances.

High precision pulsar measurements not only allow us to probe the Fresnel
corrections, but they also help localize the source through the combined
phase-frequency parallax of all of the pulsars in the PTA. The previous
studies of source localization by Deng & Finn and Corbin & Cornish had shown
many orders of magnitude difference in the sky angle ranging from
$\mathcal{O}(100\ \text{arcsec}^{2}\ -\ 10\text{deg}^{2})$. In this work we
explain that this difference is due to many factors including the timing
accuracy $\sigma$, the pulsar distance accuracy $\sigma_{L}$, the PTA scale
factor, and the intrinsic source parameters. We also show how knowledge of the
source distance through Fresnel corrections can be combined with the sky angle
measurements in order to localize the source to a volume of space. Even for
sources near $\sim 1$ Gpc in distance could be localized to within a volume
$\Delta V<1\ \text{Mpc}^{3}$. This would allow us to localize the galaxy
source of the gravitational waves which could then be targeted with an in-
depth multi-messenger follow-up, such as what was done recently in Arzoumanian
et al. (2020).

A novelty of the Fresnel formalism is that it offers a way of measuring the
chirp mass of a monochromatic source, which is something that cannot be done
using the plane-wave formalism. As long as the amplitude, distance, and
frequency parameters can be recovered from the source, then the chirp mass can
be calculated. Furthermore, the Fresnel formalism opens a path for future
cosmological studies using PTAs. As discussed recently in D’Orazio & Loeb
(2020), when the Fresnel models are generalized to a cosmologically expanding
universe, then the distance parameter $R$ which comes from the Fresnel
corrections becomes the comoving distance $D_{c}$. If the source is also
significantly chirping, then the system’s observed (redshifted) chirp mass can
be measured. Combining this with the measured observed frequency and amplitude
parameters produces a measurement of the system’s luminosity distance $D_{L}$,
which is related to the comoving distance through $D_{L}=(1+z)D_{c}$.
Therefore in the full Fresnel frequency evolution regime IIB, frequency
evolution effects allow for the measurement of the source $D_{L}$, and Fresnel
corrections allow for the measurement of the source $D_{c}$. If both of these
distances are recovered from the timing residual model then a measurement of
$H_{0}$ can be obtained (as shown in D’Orazio & Loeb). Or alternatively if
only one of these two distances can be measured accurately, but localization
of the gravitational wave source and multimessenger counterparts identify the
host galaxy, then the source host galaxy’s redshift can be measured. Combining
the redshift with either of the measured source distances once again provides
a measurement of $H_{0}$.

While we acknowledge that many of our “best-case-scenarios” required
simulating idealized PTAs well beyond our current experimental capabilities,
timing and pulsar measurements will continue to improve (Lam et al., 2018),
and we do see a lot of potential for future studies that could be made
possible from the Fresnel timing regimes. For example, seeing the direct
impact that PTA size has on parameter estimation hopefully offers strong
motivation to actively search for and include more distant pulsars with higher
precision distance measurements in our current PTAs. Finally, we reiterate
here that in this work we have explicitly assumed a flat static universe for
the sake of simplicity and to help isolate the measurement of the distance
parameter $R$ (see assumption 1). More generally, for an expanding universe
the fully generalized Fresnel frequency evolution model contains both a
luminosity and a comoving distance parameter. We will explore this generalized
model in much greater detail in a follow-up study to this paper, and the
intriguing cosmological implications that the Fresnel regime may offer to
future PTAs.

## Acknowledgements

We thank Luke Kelley for helpful comments during the preparation of this
paper. This work was supported by the National Science Foundation (NSF)
through the NANOGrav collaboration Physics Frontier Center award NSF
PHY-1430284, as well as the awards NSF PHY-1607585 and NSF PHY-1912649. This
work also received support from the Wisconsin Space Grant Consortium Graduate
and Professional Research Fellowship Program under the NASA Training Grant
#NNX15AJ12H, as well as from the University of Wisconsin-Milwaukee, for which
we are grateful for access to their computational resources supported by NSF
PHY-1626190.

## Data Availability

The 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, The Astrophysical Journal, 880, 116
  * Arzoumanian et al. (2014) Arzoumanian Z., et al., 2014, The Astrophysical Journal, 794, 141
  * Arzoumanian et al. (2018) Arzoumanian Z., et al., 2018, The Astrophysical Journal Supplement Series, 235, 37
  * Arzoumanian et al. (2020) Arzoumanian Z., et al., 2020, ApJ, 900, 102
  * Belsley et al. (1980) Belsley D. A., Kuh E., Welsch R. E., 1980, Regression Diagnostics: Identifying Influential Data and Sources of Collinearity. John Wiley & Sons, Ltd, doi:10.1002/0471725153
  * Coe (2009) Coe D., 2009, Fisher Matrices and Confidence Ellipses: A Quick-Start Guide and Software (arXiv:0906.4123)
  * Corbin & Cornish (2010) Corbin V., Cornish N. J., 2010, arXiv e-prints, p. arXiv:1008.1782
  * Cordes & Shannon (2010) Cordes J. M., Shannon R. M., 2010, arXiv e-prints, p. arXiv:1010.3785
  * 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 (2020) D’Orazio D. J., Loeb A., 2020, arXiv e-prints, p. arXiv:2009.06084
  * Deller et al. (2019) Deller A. T., et al., 2019, The Astrophysical Journal, 875, 100
  * Deng & Finn (2011) Deng X., Finn L. S., 2011, Monthly Notices of the Royal Astronomical Society, 414, 50
  * Ellis (2013) Ellis J. A., 2013, Classical and Quantum Gravity, 30, 224004
  * Finn (2009) Finn L. S., 2009, Phys. Rev. D, 79, 022002
  * Holz & Hughes (2005) Holz D. E., Hughes S. A., 2005, The Astrophysical Journal, 629, 15
  * Kerr et al. (2020) Kerr M., et al., 2020, Publ. Astron. Soc. Australia, 37, e020
  * Lam et al. (2018) Lam M. T., McLaughlin M. A., Cordes J. M., Chatterjee S., Lazio T. J. W., 2018, The Astrophysical Journal, 861, 12
  * Liu et al. (2011) Liu K., Verbiest J. P. W., Kramer M., Stappers B. W., van Straten W., Cordes J. M., 2011, Monthly Notices of the Royal Astronomical Society, 417, 2916
  * Lorimer & Kramer (2004) Lorimer D. R., Kramer M., 2004, Handbook of Pulsar Astronomy. Vol. 4
  * Ly et al. (2017) Ly A., Marsman M., Verhagen J., Grasman R., Wagenmakers E.-J., 2017, Journal of Mathematical Psychology, 80, 40–55
  * 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
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, The Astronomical Journal, 129, 1993
  * Schutz (1986) Schutz B. F., 1986, Nature, 323, 310
  * Wittman (ND) Wittman D., N.D., Fisher Matrix for Beginners, UC Davis, http://wittman.physics.ucdavis.edu/Fisher-matrix-guide.pdf

##  Appendix A Deriving the Gravitational Wave Timing Residual

In this section we give a brief overview of the calculations which go into
modeling the gravitational wave induced timing residual for a pulsar timing
(or one-way Doppler tracking) experiment. There exist various approaches to
the derivation of the residual in the literature (Creighton & Anderson, 2011;
Finn, 2009). The approach we take is to study the path of a photon by
integrating the spacetime metric, as shown in Maggiore (2018). Finn (2009)
notes that, under certain circumstances, the method used here will not be
valid; however, as we assume a flat Minkowski background (see below) and
express the metric in the transverse-traceless gauge, the conditions for
validity required in Finn (2009) are satisfied. Below we indicate the
assumptions which go into this work and our derivations.

###  A.1 Assumptions

Einstein Field Equations & Source Binary

  1. 1.

Cosmologically static and flat universe (hence the background metric is flat
Minkowski, ${\eta_{\mu\nu}=\mathrm{diag}(-c^{2},\ 1,\ 1,\ 1)}$).

  2. 2.

Weak field limit (gravitational waves are a metric perturbation on top of the
Minkowski background).

  3. 3.

Small source compared to the distance to the observer and the wavelength of
the wave.

  4. 4.

Slow-moving source (non-relativistic, no post-Newtonian analysis required).

  5. 5.

Transverse-traceless gauge.

  6. 6.

The binary source is circular.

  7. 7.

The source is located at a fixed angular sky position and distance from the
Earth.

  8. 8.

Far field approximation of the metric perturbation amplitude \- the binary
source is sufficiently far from the field point of interest that
$|\vec{x}_{\mathrm{field}}-\vec{x}_{\mathrm{source}}|\approx R$ in the
amplitude of the metric perturbation. This assumption won’t be made when
evaluating the retarded time.

Earth & Pulsar

  1. 10.

Earth is at the center of our coordinate system.

  2. 11.

The pulsar is located at a fixed angular sky position and distance from the
Earth.

  3. 12.

The antenna patterns are assumed to remain constant over the Earth-pulsar
baseline.

  4. 13.

The pulsar’s own rotation period $T$ timescale must be appropriately small in
comparison to the orbital period and coalescence time of the gravitational
wave source. Specifically we have two case requirements depending on the
regime:

     * •

Monochromatic Source: $\omega_{0}T\sim\frac{T}{T_{\mathrm{orbit}}}\ll 1$

     * •

Frequency Evolving Source: $\frac{T}{\Delta\tau_{c}}\ll 1$

Our pulsars are millisecond pulsars (Arzoumanian et al., 2018) and the orbital
period of our sources of interest are $\mathcal{O}\left(\text{weeks}\ -\
\text{decades}\right)$, so this assumption is safe.

Frequency Evolution Formalism

  1. 14.

The observation time scale is much, much less than the time to coalescence
(measured from the fiducial time): $t_{\mathrm{obs}}\ll\Delta\tau_{c}$

Fresnel Formalism

  1. 15.

When deriving the Fresnel formalism we assume we are not working in any of the
plane-wave limits, that is when $F=0$ or $\ \hat{r}\cdot\hat{p}=\pm 1$.

  2. 16.

Corrections in the geometrical orientation components which affect the
amplitude of the timing residual, namely the antenna patterns, which are due
to the curvature of a gravitational wave are negligible. Only corrections
which directly affect the phase and frequency of these quantities are
significant.

###  A.2 The Gravitational Wave Timing Residual

We solve for the model of an induced gravitational wave timing residual on a
pulsar by first considering the time of arrival of photons in successive
pulsar periods. Photons will travel the radial null geodesic between the
pulsar and Earth. We start with the spacetime metric:

$0=-(cdt)^{2}+\left[\eta_{ij}+h^{TT\hat{r}}_{ij}\right]\hat{p}^{i}\hat{p}^{j}dr^{2},$
(58)

where
$h^{TT\hat{r}}_{ij}=h^{TT\hat{r}}_{ij}\big{(}t_{\mathrm{ret}}(t,\vec{x})\big{)}$
is evaluated at the retarded time of the gravitational wave. The superscript
$TT\hat{r}$ denotes that we are in the transverse-traceless gauge along some
axis $\hat{r}$, which is the direction to our source, and $\hat{p}$ is the
direction to our pulsar. The photon’s path to zeroth order in the metric
perturbation is
$\vec{x}_{0}(t)\approx\left[ct_{\mathrm{obs}}-ct\right]\hat{p}$, and we
integrate along that path between the endpoints
${(t,r)=(t_{\mathrm{obs}}-\frac{L}{c},L)\rightarrow(t_{\mathrm{obs}},0)}$.
This gives us an expression for the time the first photons are observed on
Earth after they were emitted from the pulsar. The integral is then repeated
for the photons one pulsar period $T$ later. The difference in these two
expressions gives us the observed pulsar period on Earth,
$T_{\mathrm{obs}}\approx T+\Delta T$, where:

$\Delta
T=\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}}\bigg{[}h_{\textsc{A}}\Big{(}t_{\mathrm{ret}}\big{(}t+T,\vec{x}_{0}(t)\big{)}\Big{)}-h_{\textsc{A}}\Big{(}t_{\mathrm{ret}}\big{(}t,\vec{x}_{0}(t)\big{)}\Big{)}\bigg{]}dt.$
(59)

Dividing both sides of this expression by $T$ we can write $\frac{\Delta
T}{T}$, which is the gravitational wave induced fractional shift in the pulsar
period. Next we take the resulting integrand and recognize that as long as the
non-dimensional quantities indicated in assumption 13 are small, then we can
approximate that integrand as as a derivative. This allows us to now write:

$\frac{\Delta