literature, and in my work I add a new “Fresnel” regime.

The key to this derivation of the Fresnel formalism is that I keep out to the
Fresnel term in the expansion of the retarded time in equation 3.21.
Conceptually, we account for 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 Figure 3.4 shows).
This will slightly alter the frequency and phase of the timing residual in the
pulsar term, which later we will show can produce a measurable effect.

It is important to realize that this derivation and formalism does not correct
for wavefront curvature effects in the amplitude or the geometrical
orientation terms like the antenna patterns in Section 3.5.1. In reality, the
pulsar-source distance is different from the Earth-source distance $R$. This
means that the amplitude of the metric perturbation (equation 2.23) should be
different for the photon when it leaves the pulsar versus when it arrives at
the Earth. Furthermore, since the pulsar is at a different location, the
pulsar-source geometry will be different, and hence a photon leaving the
pulsar will be affected by a metric perturbation with different orientation
angles $\iota$ and $\psi$ (and $\theta_{0}$ since that is also an Euler angle
as discussed in Section 2.5). We could Taylor expand
$|\vec{x}-\vec{x}^{{}^{\prime}}|$ in the small parameter $\frac{L}{R}$ and
look for $\mathcal{O}\left(\frac{L}{R}\right)$ corrections in the orientation
geometry as well, but all of these corrections will directly affect the
amplitude of the metric perturbation and timing residual. For a typical
$L=1$kpc pulsar, even a source at $R=1$Mpc would only give a correction of
order $\mathcal{O}(0.001)$. In general we just don’t expect any of these
geometrical curvature effects to produce any measurable effect in the timing
residual.

However, corrections of order $\mathcal{O}\left(\frac{L}{R}\right)$ in the
frequency and phase of the metric perturbation and timing residual are
different. This is primarily because the phase of a sinusoidal function wraps
around the interval $[0,2\pi)$. Additionally, as we discussed in Section 3.3,
corrections of order $\mathcal{O}\left(\frac{L}{R}\right)$ in the retarded
time combine with factors $\omega_{\mathrm{gw}}\frac{L}{c}$ to produce net
corrections of the order of the Fresnel number. Overall, we find that these
corrections on a $2\pi$-interval can produce measurable effects in the overall
timing residual. It is one of the main goals of my work to study when and how
measureable these effects will be in this new Fresnel regime.

###  3.5 Plane-Wave Formalism (I)

We start by keeping only terms in the plane-wave regime in equation 3.27:

$\displaystyle t_{\mathrm{ret}}^{0}$ $\displaystyle\approx
t-\frac{R}{c}+\left(\hat{r}\cdot\hat{p}\right)(t_{\mathrm{obs}}-t),$
$\displaystyle=\left[-\frac{R}{c}+\left(\hat{r}\cdot\hat{p}\right)t_{\mathrm{obs}}\right]+\left(1-\hat{r}\cdot\hat{p}\right)t.$
(3.28)

The integral in equation 3.13 is over the dummy time variable $t$. The
explicit dependence on $t$ in that expression is in $t^{0}_{\mathrm{ret}}$,
which in this section is given in equation 3.28. Again remember, conceptually
our time variable $t$ is defined for the interval $t_{\mathrm{obs}}-L/c\leq
t\leq t_{\mathrm{obs}}$ which gives us the path of the photon as it traverses
between the pulsar and the Earth.

As we found in equation 3.13, there are a number of useful ways to express the
fractional change in the pulsar’s period. We could approach solving the
integral using any of these forms. However, in the plane-wave formalism, it
turns out the solution is rather remarkably simple and elegant once we make
the following critical observation. Let’s change the integration variable from
$t\rightarrow t^{0}_{\mathrm{ret}}$ using equation 3.28, which tells us
$\frac{dt_{\mathrm{ret}}^{0}}{dt}=\left(1-\hat{r}\cdot\hat{p}\right)$. This
one statement is actually the crux of the entire plane-wave formalism -
because this derivative is just a constant, it can immediately be pulled out
of the integral. The result is:

$\displaystyle\frac{\Delta T}{T}(t_{\mathrm{obs}})$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{obs}}}_{t_{\mathrm{obs}}-L/c}\frac{dh_{\textsc{A}}(t^{0}_{\mathrm{ret}})}{dt^{0}_{\mathrm{ret}}}dt,$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{obs}}-R/c}_{t_{\mathrm{obs}}-R/c-\left(1-\hat{r}\cdot\hat{p}\right)L/c}\frac{dh_{\textsc{A}}(t^{0}_{\mathrm{ret}})}{dt^{0}_{\mathrm{ret}}}\frac{dt^{0}_{\mathrm{ret}}}{\left(1-\hat{r}\cdot\hat{p}\right)},$
$\displaystyle=\frac{1}{2}\frac{\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}\int\limits^{h_{\textsc{A}}\left(t_{\mathrm{obs}}-R/c\right)}_{h_{\textsc{A}}\left(t_{\mathrm{obs}}-R/c-\left(1-\hat{r}\cdot\hat{p}\right)L/c\right)}dh_{\textsc{A}},$
$\displaystyle\equiv\frac{1}{2}F^{\textsc{A}}\int\limits^{h_{\textsc{A}}\left(t_{\mathrm{obs}}-R/c\right)}_{h_{\textsc{A}}\left(t_{\mathrm{obs}}-R/c-\left(1-\hat{r}\cdot\hat{p}\right)L/c\right)}dh_{\textsc{A}}.$
(3.29)

The function $F^{\textsc{A}}$ is named the antenna response function and is
discussed in detail in the next section 3.5.1. In this final expression we see
that all dependence on the actual functional form of the metric perturbation
$h_{\textsc{A}}$ actually drops out of the integration. The importance of this
statement is that it means the fractional change in the pulsar’s period does
not depend on the intermediate motion of the photon between the pulsar and the
Earth, only on those two endpoints of its journey - when the pulsar’s light
first left the pulsar, and when it finally arrived at the Earth. These
endpoints we denote with “E” and “P” subscripts. So the timing residual in the
plane-wave regime breaks apart into an “Earth term” and a “pulsar term.”
Furthermore, this result is independent of whether or not the source is a
monochromatic or a frequency evolving gravitational wave. Therefore,
completing the integration we have have the general result:

$\displaystyle\frac{\Delta T}{T}(t_{\mathrm{obs}})$
$\displaystyle=\frac{1}{2}F^{\textsc{A}}\Bigg{[}h_{\textsc{A}}\left(t_{\mathrm{obs}}-\frac{R}{c}\right)-h_{\textsc{A}}\left(t_{\mathrm{obs}}-\frac{R}{c}-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\right)\Bigg{]},$
$\displaystyle\equiv\frac{1}{2}F^{\textsc{A}}\bigg{[}h_{\textsc{A}}\Big{(}\omega^{0}_{E},\Theta^{0}_{E}\Big{)}-h_{\textsc{A}}\Big{(}\omega^{0}_{P},\Theta^{0}_{P}\Big{)}\bigg{]},$
(3.30)
$\displaystyle\mathrm{where}\quad\begin{cases}\omega^{0}_{E}&\equiv\omega\Big{(}t_{\mathrm{obs}}-\frac{R}{c}\Big{)},\\\
\omega^{0}_{P}&\equiv\omega\Big{(}t_{\mathrm{obs}}-\frac{R}{c}-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\Big{)},\\\
\Theta^{0}_{E}&\equiv\Theta\Big{(}t_{\mathrm{obs}}-\frac{R}{c}\Big{)},\\\
\Theta^{0}_{P}&\equiv\Theta\Big{(}t_{\mathrm{obs}}-\frac{R}{c}-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\Big{)}.\end{cases}$
(3.31)

Now we want to solve for the timing residual equation 3.19 by integrating the
fractional change in the pulsar’s period equation 3.30:

$\mathrm{Res}(t)=\frac{1}{2}F^{\textsc{A}}\int\bigg{[}h_{\textsc{A}}\Big{(}\omega^{0}_{E},\Theta^{0}_{E}\Big{)}-h_{\textsc{A}}\Big{(}\omega^{0}_{P},\Theta^{0}_{P}\Big{)}\bigg{]}dt_{\mathrm{obs}}.$
(3.32)

At this point depending on which model we are investigating, there will be a
couple of important differences (primarily conceptual) in the solution, so we
will treat the two cases separately.

####  3.5.1 The Plane-Wave Antenna Response

The quantity
$F^{\textsc{A}}\equiv\frac{\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$
has an important physical interpretation. Using the definitions
$E^{\hat{r}\textsc{A}}_{ij}$ and $e^{\hat{r}\textsc{A}}_{ij}$ from equation
2.31 and 2.29, and changing from Einstein notation to dot product notation we
can write:

$\displaystyle\begin{cases}\begin{tabular}[]{l 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)}$&$=\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)}$
,\\\
$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)}$&$=\frac{2\left(\hat{p}\cdot\hat{\theta}\right)\left(\hat{p}\cdot\hat{\phi}\right)}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\end{tabular}\end{cases}$ (3.33)
$\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{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)f^{+}+\sin(2\psi)f^{\times}\right]$
,\\\
$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.34)

These are purely geometrical functions which depend on the source location
($\theta$ and $\phi$) and orientation ($\iota$ and $\psi$) in relation to the
pulsar location ($\theta_{p}$ and $\phi_{p}$). We can think of the Earth-
pulsar baseline as an “antenna” for our experiment. These functions are
therefore an antenna response that then determines how sensitive $\frac{\Delta
T}{T}$ in equation 3.30 and $\mathrm{Res}(t)$ in equation 3.32 will be to the
geometrical alignment of our pulsar and source. We derive these relations
specifically from the plane-wave regime. Later in Section 3.6 we will find it
is not as easy to write the antenna response cleanly in this way in the
Fresnel formalism, but that asymptotically the Fresnel formalism also has the
same antenna response as the plane-wave formalism.

An example of this antenna response as a function of source sky position is
shown in Figure 3.6. It is important to note that the definitions of these
response functions in equations 3.33 and 3.34 are mathematically undefined
when the source and pulsar are perfectly aligned, i.e. when
$\hat{r}\cdot\hat{p}=+1$. (Note that if $\hat{r}\cdot\hat{p}=\pm 1$, then
$\hat{p}\cdot\hat{\theta}\equiv\hat{p}\cdot\hat{\phi}\equiv 0$ because by
definition of the basis vectors, $\hat{\theta}$ and $\hat{\phi}$ are
orthogonal to $\hat{r}$). However, the expressions for the fractional period
shift and timing residual will remain well defined. Looking first at the
fractional period shift equation 3.30, we see that in the limit
$\hat{r}\cdot\hat{p}\rightarrow 1$ the term in square brackets goes to zero.
It turns out that this goes to zero fast enough to kill the diverging antenna
response, so the overall fractional period shift approaches zero as the source
and pulsar become perfectly aligned. We will see in the next section 3.5.2
that the same thing happens to the timing residual in this limit. And in the
limit that the source and pulsar become anti-aligned,
$\hat{r}\cdot\hat{p}=-1$, all of these expressions, including the antenna
response, are well defined and also go to zero.

Finally, here are a few important case examples:

  1. 1.

If the source is both face-on $\iota=0$ and $\psi=0$, then $F^{+}=f^{+}$ and
$F^{\times}=f^{\times}$ (see Figure 3.6).

  2. 2.

If the source is both face-off $\iota=\pi$ and $\psi=0$, then $F^{+}=f^{+}$
and $F^{\times}=-f^{\times}$. Physically, the only difference in this case as
compared to case 1 is that the binary source now appears to orbit in the
reverse direction. The plus pattern is invariant to this orbital direction
reversal, but the cross pattern becomes “mirrored” (considering Figure 3.6
again, $F^{+}$ remains the same, but the blue and red areas in $F^{\times}$
are swapped).

  3. 3.

If $\iota=\pi/2$ (edge-on) and $\psi=0$, then $F^{+}=\frac{1}{2}f^{+}$ and
$F^{\times}=0$, meaning only the plus polarization of the gravitational wave
is affecting the change in the pulsar’s period.

  4. 4.

If $\iota=\pi/2$ (edge-on) and $\psi=\pi/4$, then
$F^{+}=\frac{1}{2}f^{\times}$ and $F^{\times}=0$, meaning only the cross
polarization of the gravitational wave is affecting the change in the pulsar’s
period.

Figure 3.6: An example of the antenna patterns as a function of sky position
angles $(\theta,\phi)$ of a source with $\iota=\psi=0$ (face-on, not
polarized). Note this corresponds to example 1 in Section 3.5.1, so really
these plots are also of $f^{+}$ and $f^{\times}$. The pulsar has been fixed at
$\theta_{p}=\frac{\pi}{2}$ and $\phi_{p}=0$, as indicated in the figure by the
star. The quadrupolar pattern of the gravitational wave’s effect can clearly
be seen in the antenna patterns. Here red and blue simply indicate regions of
the sky that would apply opposite signs to the fractional period shift and
timing residual quantities. So for example, a source in a blue region would
appear “blueshifted” with respect to the same source at the same time which
appears in the red region (which would appear “redshifted,” respectively). Of
course, in reality the gravitational waves cause this change to vary
sinusoidally in time in the fractional period shift and the timing residual,
so really here the blue and red colors indicate oppositely induced changes in
these quantities at any given time. The white regions show where the antenna
would not be sensitive to detecting a source.

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

For a monochromatic gravitational wave, the orbital frequency remains constant
and hence the amplitude of the metric perturbation loses its time-dependence
(equation 2.23). Explicitly in this model (see equations 2.33 and 2.34)
equation 3.31 becomes:

$\underset{t_{0}\ =\
-\frac{R}{c}}{\mathrm{choosing}}\quad\begin{cases}\omega^{0}_{E}&=\omega^{0}_{P}\quad\equiv\omega_{0},\\\
\Theta^{0}_{E}&=\theta_{0}+\omega_{0}t_{\mathrm{obs}},\\\
\Theta^{0}_{P}&=\theta_{0}+\omega_{0}\Big{(}t_{\mathrm{obs}}-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\Big{)}.\end{cases}$
(3.35)

Here we have also chosen our model’s fiducial time, primarily motivated by
convenience here. As described earlier in Section 2.7, let’s pick
$t_{0}=-\frac{R}{c}$. This will cancel out some factors of $\frac{R}{c}$ in
these expressions, but also it will keep our formula here consistent with our
later formulae, wherein we will make this same choice.

We use one last change of variable furnished by equation 3.35 to write the
timing residual equation 3.32 as:

$\displaystyle\mathrm{Res}(t)$
$\displaystyle=\frac{1}{2}F^{\textsc{A}}\left[\int
h_{\textsc{A}}\left(\omega_{0},\Theta_{E}^{0}\right)dt_{\mathrm{obs}}\quad-\int
h_{\textsc{A}}\left(\omega_{0},\Theta_{P}^{0}\right)dt_{\mathrm{obs}}\right],$
$\displaystyle=\frac{1}{2}F^{\textsc{A}}\left[\int
h_{\textsc{A}}\left(\omega_{0},\Theta_{E}^{0}\right)\frac{d\mathrm{\Theta^{0}_{E}}}{\omega_{0}}\quad-\int
h_{\textsc{A}}\left(\omega_{0},\Theta_{P}^{0}\right)\frac{d\mathrm{\Theta^{0}_{P}}}{\omega_{0}}\right].$
(3.36)

Because the orbital frequency is constant it can be pulled out of the
integral, and the remaining integration in equation 3.36 then becomes the
simple integration of the sine and cosine terms in eqn 2.22. Using the helpful
identities in equation 3.14, the final result 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{]},$
(3.37) $\displaystyle\mathrm{for}\quad\textsc{A}\in[+,\times],$
$\displaystyle\underset{\left(t_{0}\ =\
-\frac{R}{c}\right)}{\mathrm{where}}\quad\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}$
(3.38) $\displaystyle\hskip 56.9055pt\begin{cases}h_{+}\equiv-
h_{0}\cos\big{(}2\Theta\big{)},\\\ h_{\times}\equiv-
h_{0}\sin\big{(}2\Theta\big{)},\end{cases}$ (3.39) $\displaystyle\hskip
56.9055pth_{0}\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega_{0}^{2/3},$
(3.40) $\displaystyle\hskip 56.9055pt\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{\left(\hat{p}\cdot\hat{\theta}\right)^{2}-\left(\hat{p}\cdot\hat{\phi}\right)^{2}}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\\\
$f^{\times}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}e^{\hat{r}\times}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\frac{2\left(\hat{p}\cdot\hat{\theta}\right)\left(\hat{p}\cdot\hat{\phi}\right)}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\end{tabular}\end{cases}$ (3.33 r) $\displaystyle\hskip
56.9055pt\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{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)f^{+}+\sin(2\psi)f^{\times}\right]$
,\\\
$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.34 r) $\displaystyle\hskip
56.9055pt\begin{cases}\begin{tabular}[]{r l l
r}$\hat{r}$&$=\big{[}\sin(\theta)\cos(\phi),$&$\sin(\theta)\sin(\phi),$&$\cos(\theta)\big{]}$,\\\
$\hat{\theta}$&$=\big{[}\cos(\theta)\cos(\phi),$&$\cos(\theta)\sin(\phi),$&$\sin(\theta)\big{]}$,\\\
$\hat{\phi}$&$=\big{[}-\sin(\phi),$&$\cos(\phi),$&$0\big{]}$,\\\
$\hat{p}$&$=\big{[}\sin(\theta_{p})\cos(\phi_{p}),$&$\sin(\theta_{p})\sin(\phi_{p}),$&$\cos(\theta_{p})\big{]}$.\end{tabular}\end{cases}$
(2.28, 3.3 r)

In this model regime, looking at eqn 3.38 we see that the Earth and pulsar
terms have the same frequencies. However, their difference is in the phase
offset between the two. Comparing both of these expressions, we see that the
pulsar term $\Theta_{p}$ differs in phase by
$-\omega_{0}\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}$ compared to the
Earth term $\Theta_{E}$.

As mentioned in Section 3.5.1, the timing residual expression equation 3.37
remains well defined in both limits $\hat{r}\cdot\hat{p}\rightarrow\pm 1$
(which in Section 3.6 we will refer to as the “natural plane-wave” limits),
despite the fact that the antenna response is not well defined in the limit
that $\hat{r}\cdot\hat{p}\rightarrow+1$. As the system becomes anti-aligned
$\hat{r}\cdot\hat{p}\rightarrow-1$ the antenna patterns
$F^{\textsc{A}}\rightarrow 0$ and this kills the entire timing residual. And
as the system becomes aligned $\hat{r}\cdot\hat{p}\rightarrow+1$ the pulsar
phase $\Theta_{P}\rightarrow\Theta_{E}$ and the difference in the metric
perturbations at the Earth and pulsar goes to zero fast enough to also kill
the entire timing residual. Conceptually the photons “surf” the gravitational
waves in this second case, and thus the photons is never “feel” the
gravitational wave itself alter their own motion between the pulsar and Earth.

Note that in terms of our notation, we are writing the residual as a function
$t$, so we are changing notation in equation 3.35 from
$t_{\mathrm{obs}}\rightarrow t$ through the integration. It is also important
to realize that in this regime the source chirp mass $\mathcal{M}$ and
distance $R$ parameters are fully degenerate in the timing residual - they
only appear in the combination $\frac{\mathcal{M}^{5/3}}{R}$ with each other
in the amplitude terms here (see equation 3.40). Therefore we would not be
able to measure these parameters independently through actual observations
(only this specific combination of the two).

After applying a couple of trigonometric identities (namely,
$\cos(x)-\cos(x+a)\equiv
2\sin\left(\frac{a}{2}\right)\sin\left(\frac{a}{2}+x\right)$, and similar for
the sine identity), an alternative but equivalent way of writing equation 3.37
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],$ (3.41)
$\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}$
(3.42) $\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}$
(3.43)

The benefit of this notation is that it allows us to define in equation 3.42 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 heuristic
Fresnel monochromatic model IIA, in Section 3.6.1.

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

Here we will make the following assumptions:

Assumptions: Frequency Evolving Timing Residuals 1. The observation time scale
is much, much less than the time of coalescence (measured from the fiducial
time): $t_{\mathrm{obs}}\ll\Delta\tau_{c}$

The frequency evolution model of a gravitational wave described in Section
2.7.2 is more realistic as it captures the changing frequency of the orbit of
the binary. Conceptually, as the photon first leaves the pulsar it will feel
gravitational waves of some initial frequency which were emitted by the
source. As the pulsar-Earth baseline is on the order of kiloparsecs in
distance, the photon will then take thousands of years before it arrives at
the Earth. During this time, the binary will lose energy due to the emission
of gravitational waves and its orbit will collapse, increasing the overall
frequency of the emitted gravitational waves with time. So as the photon
travels toward the Earth, it will continuously experience gravitational waves
of higher frequencies. However, crucially, as we found in equation 3.30 and
discussed, in the plane-wave formalism the result of the fractional change in
the pulsar’s period (and subsequently, in the observed timing residual), only
depends on the frequency the photon experiences at the two endpoints, when it
initially leaves the pulsar and when it finally arrives at Earth. In the
monochromatic model, these frequencies are the same since the orbit does not
coalesce, but in this model the thousands of years travel can potentially
alter the frequency felt by the photon at the endpoints in a significant way,
as well will see here. Therefore in this model, the Earth and pulsar terms
will not only have different phases as we saw in the plane-wave monochromatic
regime, but also different frequencies.

Explicitly in this model (see equations 2.37 and 2.38) equation 3.31 becomes:

$\underset{t_{0}\ =\
-\frac{R}{c}}{\mathrm{choosing}}\quad\begin{cases}\omega^{0}_{E}&=\omega_{0}\left[1-\frac{t_{\mathrm{obs}}}{\Delta\tau_{c}}\right]^{-3/8},\\\
\omega^{0}_{P}&=\omega_{0}\left[1-\frac{t_{\mathrm{obs}}-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}}{\Delta\tau_{c}}\right]^{-3/8},\\\
\Theta^{0}_{E}&=\theta_{0}+\theta_{c}\left[1-\left(1-\frac{t_{\mathrm{obs}}}{\Delta\tau_{c}}\right)^{5/8}\right],\\\
\Theta^{0}_{P}&=\theta_{0}+\theta_{c}\left[1-\left(1-\frac{t_{\mathrm{obs}}-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}}{\Delta\tau_{c}}\right)^{5/8}\right].\end{cases}$
(3.44)

Again here we have made the choice to set the fiducial time
$t_{0}=-\frac{R}{c}$. This means that below in equation 3.45 when we specify
the value of $\omega_{0}$, we are directly specifying the value of the
frequency of the gravitational wave affecting the Earth end of the photon’s
journey, which is arguably a conceptual convenience.

Integrating the timing residual in equation 3.32 would be difficult using the
frequency evolution model since the orbital frequency is no longer a constant,
but a function of $t_{\mathrm{obs}}$. However, this is where we make our final
important assumption. Fortunately, while the frequency evolution model allows
the frequency of the gravitational wave to evolve appreciably in the thousands
of years it takes the photon to travel from the pulsar to the Earth, we still
make the assumption that on the observation time scales of our experiment
itself (which lasts on the order of years to decades) the frequency remains
constant. This means we can Taylor expand the functions in the small parameter
$\frac{t_{\mathrm{obs}}}{\Delta\tau_{c}}$ equation 3.44, the result being:

$\underset{t_{0}\ =\
-\frac{R}{c}}{\mathrm{choosing}}\quad\begin{cases}\begin{tabular}[]{ l l l
}$\omega^{0}_{E}$&$\approx\omega_{0}$&$\equiv\omega_{0E}$ ,\\\
$\omega^{0}_{P}$&$\approx\omega_{0}\left[1+\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}}{\Delta\tau_{c}}\right]^{-3/8}$&$\equiv\omega_{0P}$
,\\\
$\Theta^{0}_{E}$&$\approx\theta_{0}+\omega_{0}t_{\mathrm{obs}}$&$\equiv\theta_{0E}+\omega_{0E}t_{\mathrm{obs}}$
,\\\
$\Theta^{0}_{P}$&$\approx\left[\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)\right]+\omega_{0P}t_{\mathrm{obs}}$&$\equiv\theta_{0P}+\omega_{0P}t_{\mathrm{obs}}$
,\end{tabular}\end{cases}$ (3.45)

In the orbital frequency we only keep to the zeroth order term because again
we are assuming it doesn’t evolve in time on our observation time scales, and
in the phase we keep out to linear order in
$\frac{t_{\mathrm{obs}}}{\Delta\tau_{c}}$ (note the $\Delta\tau_{c}$ cancels
in the algebra leaving it to linear order in $t_{\mathrm{obs}}$ as we see
here).

Therefore, on observation time scales, the gravitational wave affecting the
endpoints of the pulsar-Earth baseline is approximately monochromatic
(comparing equations 3.35 and 3.45 we see they are of the same basic form, a
constant frequency and a phase linearly increasing in time). This allows us to
employ the same change of variables that we did when solving equation 3.36 in
the monochromatic formalism. The only difference here is that now when we
consider equation 3.45, we see that there are two frequencies, an initial
“pulsar” frequency and an initial “Earth” frequency, corresponding to the
respective frequencies of the gravitational waves which arrived at those two
endpoints as the photon made its journey from the pulsar to the Earth (when we
initially begin our experiment). In the monochromatic model those two
frequencies were of course the same, because in that model the frequency never
changes, but here the difference between the two frequencies is the result of
the slowly changing orbital frequency over the thousands of years travel for
the photon between the pulsar and Earth.

Therefore making the appropriate change of variables through equation 3.45
(mirroring what we did earlier in the monochromatic case with equation 3.36),
now we can write equation 3.32 as:

$\displaystyle\mathrm{Res}(t)$
$\displaystyle\approx\frac{1}{2}F^{\textsc{A}}\left[\int
h_{\textsc{A}}\left(\omega_{0E},\Theta_{E}^{0}\right)dt_{\mathrm{obs}}\quad-\int
h_{\textsc{A}}\left(\omega_{0P},\Theta_{P}^{0}\right)dt_{\mathrm{obs}}\right],$
$\displaystyle\approx\frac{1}{2}F^{\textsc{A}}\left[\int
h_{\textsc{A}}\left(\omega_{0E},\Theta_{E}^{0}\right)\frac{d\mathrm{\Theta^{0}_{E}}}{\omega_{0E}}\quad-\int
h_{\textsc{A}}\left(\omega_{0P},\Theta_{P}^{0}\right)\frac{d\mathrm{\Theta^{0}_{P}}}{\omega_{0P}}\right].$
(3.46)

Just like we found before, the orbital frequencies in each integral are
constants so they can be pulled out of the integration. And like before, the
remaining integration in equation 3.46 then becomes the simple integration of
the sine and cosine terms in equation 2.22. Using the helpful identities in
equation 3.14, the final result 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],$
(3.47) $\displaystyle\mathrm{for}\quad\textsc{A}\in[+,\times],$
$\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}$ (3.48) $\displaystyle\hskip
45.52458pt\begin{cases}\Delta\tau_{c}\equiv\frac{5}{256}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\frac{1}{\omega_{0}^{8/3}},\\\
\theta_{c}\equiv\frac{8}{5}\Delta\tau_{c}\omega_{0},\end{cases}$ (2.36, 2.39
r) $\displaystyle\hskip
45.52458pt\begin{cases}h_{+}\equiv-h(\omega)\cos\big{(}2\Theta\big{)},\\\
h_{\times}\equiv-h(\omega)\sin\big{(}2\Theta\big{)},\end{cases}$ (2.22 r)
$\displaystyle\hskip
45.52458pth(\omega)\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega^{2/3},$
(2.23 r) $\displaystyle\hskip 45.52458pt\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{\left(\hat{p}\cdot\hat{\theta}\right)^{2}-\left(\hat{p}\cdot\hat{\phi}\right)^{2}}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\\\
$f^{\times}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}e^{\hat{r}\times}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\frac{2\left(\hat{p}\cdot\hat{\theta}\right)\left(\hat{p}\cdot\hat{\phi}\right)}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\end{tabular}\end{cases}$ (3.33 r) $\displaystyle\hskip
45.52458pt\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{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)f^{+}+\sin(2\psi)f^{\times}\right]$
,\\\
$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.34 r) $\displaystyle\hskip
45.52458pt\begin{cases}\begin{tabular}[]{r l l
r}$\hat{r}$&$=\big{[}\sin(\theta)\cos(\phi),$&$\sin(\theta)\sin(\phi),$&$\cos(\theta)\big{]}$,\\\
$\hat{\theta}$&$=\big{[}\cos(\theta)\cos(\phi),$&$\cos(\theta)\sin(\phi),$&$\sin(\theta)\big{]}$,\\\
$\hat{\phi}$&$=\big{[}-\sin(\phi),$&$\cos(\phi),$&$0\big{]}$,\\\
$\hat{p}$&$=\big{[}\sin(\theta_{p})\cos(\phi_{p}),$&$\sin(\theta_{p})\sin(\phi_{p}),$&$\cos(\theta_{p})\big{]}$.\end{tabular}\end{cases}$
(2.28, 3.3 r)

Once again note that in terms of our notation, we are writing the residual as
a function $t$, so we change notation in equation 3.35 from
$t_{\mathrm{obs}}\rightarrow t$ through the integration. As explained earlier,
from equation 3.48 we can now see explicitly how the Earth and pulsar terms
differ in this regime in both phase and frequency. Because the Earth and
pulsar terms combine in equation 3.47, in these cases where frequency
evolution is non-negligible we expect this to physically produce “beats” in
the timing residual data that we observe.

Note, 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 as well as the amplitude, and not in
the same combination with $R$ everywhere. So in principle a loud enough
observed source with significant frequency evolution should allow for
independent measurements of both of these source parameters.

Also, as a sanity check, we can see that this result equation 3.47 reduces to
the monochromatic result equation 3.37 in the appropriate limit, as we would
expect. For the source to produce monochromatic gravitational waves, the time
of coalescence must be infinite. Therefore, we take the limit that
$\Delta\tau_{c}\rightarrow\infty$, and see that this immediately implies
$\omega_{0P}\rightarrow\omega_{0}$ as we would expect (now both endpoints,
pulsar and Earth, feel the same frequency). But it also implies
$\theta_{c}\rightarrow\infty$, so we must be careful when taking the limit of
$\theta_{0P}$ in order to do it correctly. In the large $\Delta\tau_{c}$ limit
we can Taylor expand $\theta_{0P}$ and note the following cancellations:

$\displaystyle\theta_{0P}$
$\displaystyle=\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),$
$\displaystyle=\theta_{0}+\frac{8}{5}\Delta\tau_{c}\omega_{0}\left(1-\left[1+\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}}{\Delta\tau_{c}}\right]^{5/8}\right),$
$\displaystyle\approx\theta_{0}+\frac{8}{5}\Delta\tau_{c}\omega_{0}\left(1-\left[1+\frac{5}{8}\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}}{\Delta\tau_{c}}\right]\right),$
$\displaystyle=\theta_{0}-\omega_{0}\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}.$

Therefore in this limit we have
$\Theta_{P}\rightarrow\theta_{0}+\omega_{0}\Big{(}t-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\Big{)}$.
This along with the expression for $\Theta_{E}$ now matches equation 3.38.

A chronological summary of events and times is shown in Figure 3.7. While this
figure is specifically for the plane-wave formalism, it also holds in general
for the Fresnel formalism, which would only slightly alter the time the
gravitational waves are produced that hit the photon as it leaves the pulsar.

Figure 3.7: Here we give a conceptual chronological summary of the important
times and scales involved in the calculation of our timing residuals. This
also shows the terminology we are using for the various measured quantities of
the source’s orbital phase and frequency, with the time they correspond to.
The fiducial time can be chosen at any time, and in this work we choose it to
be at $t_{0}=-\frac{R}{c}$, which means that $\theta_{0E}=\theta_{0}$ and
$\omega_{0E}=\omega_{0}$. The coalescence time $\tau_{c}$ is measured from
present day, whereas $\Delta\tau_{c}$ is measured from the fiducial time. We
begin the experiment at $t=0$, and at this time our measurement of the pulsar
timing residual depends on the initial observation time $t_{\mathrm{obs}}$ and
the time the photon first left the pulsar, $t_{\mathrm{obs}}-\frac{L}{c}$.
Furthermore we show how these times are mapped to the retarded times
corresponding to when the gravitational waves were produced,
$t^{0}_{\mathrm{ret}}\left(t_{\mathrm{obs}}\right)$ and
$t^{0}_{\mathrm{ret}}\left(t_{\mathrm{obs}}-\frac{L}{c}\right)$. These time
scales are continually integrated as we observe the pulsar for the duration of
our experiment, so we can think of them as “sliders” of fixed width that move
forward in time. The time at which the gravitational waves are produced that
hit the photon as it leaves the pulsar depends on the Earth-pulsar-source
geometry (indicated by the $\hat{r}\cdot\hat{p}$ term). It will fall somewhere
within the range $t=-\frac{R}{c}$ (if aligned) to
$t=-\frac{R}{c}-2\frac{L}{c}$ (if anti-aligned). Finally we see the
retardation delay caused by a plane-wave gravitational wavefront. In the case
of the Fresnel formalism, additional curvature terms would factor in here.
Namely we see in our heuristic model in Section 3.6.1 that the time the
gravitational waves produced which hit the photon as it leaves the pulsar
would become approximately
$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}$.

###  3.6 Fresnel Formalism (II)

For the work below we will make the following assumptions:

Assumptions: Fresnel Formalism 1. We are not working in a plane-wave limit.
The Fresnel formalism will reduce to the plane-wave formalism in the following
limits: (a) $F\rightarrow 0$. As the Fresnel number between the pulsar and
source decreases, the Fresnel effects become more and more negligible. (b)
$\hat{r}\cdot\hat{p}\rightarrow\pm 1$. These are the natural plane-wave
limits. In this case the source becomes perfectly aligned/anti-aligned with
the pulsar, which means there are geometrically no curvature effects possible
in this precise alignment. 2. 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.

The key to deriving this model is that we keep out to the Fresnel term in the
expansion of the retarded time in equation 3.21. Conceptually, we account for
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 Figure 3.4 shows). This will slightly alter the
frequency and phase of the timing 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, as stated in
assumption 2. We expect those corrections to introduce curvature effects of
order $\mathcal{O}\left(\frac{L}{R}\right)$ in the amplitudes of the terms
that appear in the timing residuals, which we do not expect to be measurable.
Curvature corrections to the phase and frequency terms we will show later can
be measured.

Now we return to equation 3.27 and this time we keep out to the Fresnel regime
and drop all higher order terms after that. The result can be written in the
form:

$\displaystyle t^{0}_{\mathrm{ret}}$
$\displaystyle\approx\left(t_{\mathrm{obs}}-\frac{R}{c}\right)-\frac{R}{c}\left(1-\hat{r}\cdot\hat{p}\right)u-\frac{1}{2}\frac{R}{c}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)u^{2},$
$\displaystyle\mathrm{where}\quad u\equiv\frac{ct_{\mathrm{obs}}-ct}{R}.$
(3.49)

The key difference here this time is that if we try to perform the same change
of variables as we did before in equation 3.29, the $\frac{\partial
t^{0}_{\mathrm{ret}}}{\partial t}$ term will no longer be a constant (it will
be a function of $t$), and so we won’t be able to pull it out of the integral.
Instead let’s go back to equation 3.13 and look at the third line, using the
definition of $u$ to change variables and write the integral as the following:

$\displaystyle\frac{\Delta T}{T}$
$\displaystyle=\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{[}\frac{2}{3}h_{\textsc{A}}(t^{0}_{\mathrm{ret}})\frac{\dot{\omega}(t^{0}_{\mathrm{ret}})}{\omega^{0}}+2h_{\textsc{A}}\bigg{(}\omega^{0},\Theta^{0}+\frac{\pi}{4}\bigg{)}\omega^{0}\Bigg{]}dt,$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{0}_{L/R}\Bigg{[}\frac{2}{3}h_{\textsc{A}}(t^{0}_{\mathrm{ret}})\frac{\dot{\omega}(t^{0}_{\mathrm{ret}})}{\omega^{0}}+2h_{\textsc{A}}\bigg{(}\omega^{0},\Theta^{0}+\frac{\pi}{4}\bigg{)}\omega^{0}\Bigg{]}\frac{-R}{c}du,$
$\displaystyle=+\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\frac{R}{c}\int\limits^{L/R}_{0}\Bigg{[}\frac{2}{3}h_{\textsc{A}}(t^{0}_{\mathrm{ret}})\frac{\dot{\omega}(t^{0}_{\mathrm{ret}})}{\omega^{0}}+2h_{\textsc{A}}\bigg{(}\omega^{0},\Theta^{0}+\frac{\pi}{4}\bigg{)}\omega^{0}\Bigg{]}du.$
(3.50)

####  3.6.1 Fresnel, Monochromatic (IIA)

For a monochromatic gravitational wave the orbital frequency is a constant
$\omega^{0}\equiv\omega\left(t^{0}_{\mathrm{ret}}\right)\equiv\omega_{0}$
(from equation 2.33), which means $\dot{\omega}\equiv 0$, $h(t)=h_{0}$
(equation 3.40), and the first term in square brackets in equation 3.50
vanishes. For the second term, let’s take a look at the explicit expressions
for both the plus and cross terms. Using equation 3.15 we can write them as:

$\displaystyle 2h_{\textsc{A}}\bigg{(}\Theta^{0}+\frac{\pi}{4}\bigg{)}$
$\displaystyle=\left\\{\begin{tabular}[]{r
l}$2h_{0}\sin\left(2\Theta^{0}\right)$&$=2h_{0}\sin\left(2\Theta\left(t^{0}_{\mathrm{ret}}(u)\right)\right)$
,\\\\[8.0pt]
$-2h_{0}\cos\left(2\Theta^{0}\right)$&$=-2h_{0}\cos\left(2\Theta\left(t^{0}_{\mathrm{ret}}(u)\right)\right)$
,\end{tabular}\right.$ (3.53) $\displaystyle=\left\\{\begin{tabular}[]{l
r}$=2h_{0}\sin\left(\textsf{A}-\textsf{B}u-\textsf{C}u^{2}\right)$
,&$(+)$\\\\[8.0pt]
$=-2h_{0}\cos\left(\textsf{A}-\textsf{B}u-\textsf{C}u^{2}\right)$
,&$(\times)$\end{tabular}\right.$ (3.56) $\displaystyle\underset{\left(t_{0}\
=\ -\frac{R}{c}\right)}{\mathrm{where}}\quad\left\\{\begin{tabular}[]{l l
l}$\textsf{A}$&$\equiv 2\theta_{0}+2\omega_{0}t_{\mathrm{obs}}$ ,&\\\\[8.0pt]
$\textsf{B}$&$\equiv
2\omega_{0}\frac{R}{c}\left(1-\hat{r}\cdot\hat{p}\right)$&$=\frac{2\textsf{C}}{\left(1+\hat{r}\cdot\hat{p}\right)}$
,\\\\[8.0pt]
$\textsf{C}$&$\equiv\omega_{0}\frac{R}{c}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)$&$=\omega_{0}\frac{R}{c}\left(1-\hat{r}\cdot\hat{p}\right)\left(1+\hat{r}\cdot\hat{p}\right)$
,\\\ &&$=\frac{\left(1+\hat{r}\cdot\hat{p}\right)\textsf{B}}{2}$
.\end{tabular}\right.$ (3.61)

Once again we make the familiar choice that $t_{0}=-\frac{R}{c}$. Now we can
write equation 3.50 in the following compact form:

$\displaystyle\frac{\Delta T}{T}$
$\displaystyle=\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\omega_{0}\frac{R}{c}h_{0}g_{\textsc{A}},$
(3.62)
$\displaystyle\mathrm{where}\quad\begin{cases}g_{+}\equiv\int\limits^{L/R}_{0}\sin\left(\textsf{A}-\textsf{B}u-\textsf{C}u^{2}\right)du,\\\
g_{\times}\equiv-\int\limits^{L/R}_{0}\cos\left(\textsf{A}-\textsf{B}u-\textsf{C}u^{2}\right)du.\end{cases}$
(3.63)

As a sanity check, note that if we drop the $u^{2}$ term and carry out the
integration in equation 3.63, we arrive back at the expression in equation
3.36 and hence we do indeed regain the plane-wave monochromatic approximation
equation 3.37 as we would expect.

By appealing to complex methods, both the integrals in equation 3.63 can be
solved simultaneously by taking the real and imaginary parts of the solution
of the following expression:

$\displaystyle\int\limits_{0}^{L/R}\exp{\bigg{[}i\Big{(}\textsf{A}-\textsf{B}u-\textsf{C}u^{2}\Big{)}\Bigg{]}}du$
$\displaystyle=\int\limits_{0}^{L/R}\exp{\left[i\left(\left(\textsf{A}+\frac{\textsf{B}^{2}}{4\textsf{C}}\right)-\frac{\left(2\textsf{C}u+\textsf{B}\right)^{2}}{4\textsf{C}}\right)\right]}du,$
$\displaystyle\equiv\int\limits_{0}^{L/R}\exp{\Bigg{[}i\bigg{(}\Phi-\frac{\pi}{2}x^{2}\bigg{)}\Bigg{]}}du,$
$\displaystyle=\sqrt{\frac{\pi}{2\textsf{C}}}e^{i\Phi}\int\limits_{\eta_{1}}^{\eta_{2}}e^{-i\frac{\pi}{2}x^{2}}dx,$
$\displaystyle=\sqrt{\frac{\pi}{2\textsf{C}}}e^{i\Phi}\left[\int\limits_{0}^{\eta_{2}}e^{-i\frac{\pi}{2}x^{2}}dx-\int\limits_{0}^{\eta_{1}}e^{-i\frac{\pi}{2}x^{2}}dx\right],$
$\displaystyle=\sqrt{\frac{\pi}{2\textsf{C}}}\Big{(}\cos(\Phi)+i\sin(\Phi)\Big{)}\Big{[}C(\eta_{2})-iS(\eta_{2})-C(\eta_{1})+iS(\eta_{1})\Big{]},$
(3.64) $\displaystyle\mathrm{where}\quad$
$\displaystyle\begin{cases}\Phi\equiv\textsf{A}+\frac{\textsf{B}^{2}}{4\textsf{C}},\\\\[8.0pt]
\eta_{1}\equiv\frac{\textsf{B}}{\sqrt{2\pi\textsf{C}}},\\\\[8.0pt]
\eta_{2}\equiv\eta_{1}\Big{[}1+\frac{2\textsf{C}}{B}\left(\frac{L}{R}\right)\Big{]},\end{cases}$
(3.65)

where we make the substitution
$\sqrt{\frac{\pi}{2}}x=\frac{2\textsf{C}u+\textsf{B}}{2\sqrt{\textsf{C}}}$.
Note that the only required assumption here is that $\textsf{C}>0$. Given the
form of C in equation 3.61 this will only fail if the gravitational wave
source and the pulsar are perfectly aligned or anti-aligned, such that
$\hat{r}\cdot\hat{p}=\pm 1$. However, if this is the case and the source and
pulsar are perfectly aligned or anti-aligned, then the quadratic term in $u$
would vanish, and the problem would once again be a plane-wave problem (hence
assumption 1b).

Here we note that we can write our solution using the following definitions of
the well known Fresnel integrals:

$\displaystyle\begin{cases}S(\eta)\equiv\int\limits^{\eta}_{0}\sin\left(\frac{\pi}{2}x^{2}\right)dx\quad\xrightarrow[|\eta|\gg
1]{\mathrm{asymptotic}}\quad\frac{\mathrm{sgn}(\eta)}{2}-\frac{1}{\pi\eta}\cos\left(\frac{\pi}{2}\eta^{2}\right),\\\
C(\eta)\equiv\int\limits^{\eta}_{0}\cos\left(\frac{\pi}{2}x^{2}\right)dx\quad\xrightarrow[|\eta|\gg
1]{\mathrm{asymptotic}}\quad\frac{\mathrm{sgn}(\eta)}{2}+\frac{1}{\pi\eta}\sin\left(\frac{\pi}{2}\eta^{2}\right),\end{cases}$
(3.66)

Taking the real and imaginary parts of equation 3.64 gives us the solutions to
equation 3.63:

$\displaystyle\begin{cases}g_{+}=\sqrt{\frac{\pi}{2\textsf{C}}}\bigg{[}\Big{\\{}C\left(\eta_{2}\right)-C\left(\eta_{1}\right)\Big{\\}}\sin(\Phi)-\Big{\\{}S\left(\eta_{2}\right)-S\left(\eta_{1}\right)\Big{\\}}\cos(\Phi)\bigg{]},\\\
g_{\times}=-\sqrt{\frac{\pi}{2\textsf{C}}}\bigg{[}\Big{\\{}C\left(\eta_{2}\right)-C\left(\eta_{1}\right)\Big{\\}}\cos(\Phi)+\Big{\\{}\left(\eta_{2}\right)-S\left(\eta_{1}\right)\Big{\\}}\sin(\Phi)\bigg{]},\end{cases}$
(3.67)

Now we can solve equation 3.19 to get the timing residual. The convenient
thing here is that the only term that depends on $t_{\mathrm{obs}}$ is
$\Phi=\Phi\big{(}\textsf{A}\left(t_{\mathrm{obs}}\right)\big{)}$, which we can
use to write
$\frac{d\Phi}{dt_{\mathrm{obs}}}=\frac{d\Phi}{d\textsf{A}}\frac{d\textsf{A}}{dt_{\mathrm{obs}}}=2\omega_{0}$
and perform the following change of variables in the integration:

$\mathrm{Res}(t)=\int\left[\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\omega_{0}\frac{R}{c}h_{0}g_{\textsc{A}}\right]dt_{\mathrm{obs}}=\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\omega_{0}\frac{R}{c}h_{0}\int
g_{\textsc{A}}\frac{d\mathrm{\Phi}}{2\omega_{0}}.$ (3.68)

Finishing out the integration of $g_{\textsc{A}}$ from equation 3.67 and
tidying up the algebra (and again, integration changes
$t_{\mathrm{obs}}\rightarrow t$), the entire result can be expressed as:

$\displaystyle\mathrm{Res}(t)$
$\displaystyle=\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\sqrt{\frac{\pi
R/c}{8\omega_{0}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)}}\Bigg{[}\Big{\\{}C\left(\eta_{2}\right)-C\left(\eta_{1}\right)\Big{\\}}h_{\textsc{A}}\left(\Theta^{{}^{\prime}}\right)\Bigg{.}$
$\displaystyle\hskip 184.9429pt+\
\Bigg{.}\Big{\\{}S\left(\eta_{2}\right)-S\left(\eta_{1}\right)\Big{\\}}h_{\textsc{A}}\left(\Theta^{{}^{\prime}}-\frac{\pi}{4}\right)\Bigg{]},$
(3.69) $\displaystyle\mathrm{for}\quad\textsc{A}\in[+,\times],$
$\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),$
(3.70) $\displaystyle\hskip
42.67912pt\begin{cases}\eta_{1}\equiv\sqrt{\frac{2\omega_{0}R/c\left(1-\hat{r}\cdot\hat{p}\right)}{\pi\left(1+\hat{r}\cdot\hat{p}\right)}},\\\
\eta_{2}\equiv\eta_{1}\Big{[}1+\left(1+\hat{r}\cdot\hat{p}\right)\left(\frac{L}{R}\right)\Big{]},\end{cases}$
(3.71) $\displaystyle\hskip 42.67912pt\begin{cases}h_{+}\equiv-
h_{0}\cos\big{(}2\Theta\big{)},\\\ h_{\times}\equiv-
h_{0}\sin\big{(}2\Theta\big{)},\end{cases}$ (3.39 r) $\displaystyle\hskip
42.67912pth_{0}\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega_{0}^{2/3},$
(3.40 r) $\displaystyle\hskip 42.67912pt\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{\left(\hat{p}\cdot\hat{\theta}\right)^{2}-\left(\hat{p}\cdot\hat{\phi}\right)^{2}}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\\\
$f^{\times}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}e^{\hat{r}\times}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\frac{2\left(\hat{p}\cdot\hat{\theta}\right)\left(\hat{p}\cdot\hat{\phi}\right)}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\end{tabular}\end{cases}$ (3.33 r) $\displaystyle\hskip
42.67912pt\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{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)f^{+}+\sin(2\psi)f^{\times}\right]$
,\\\
$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.34 r) $\displaystyle\hskip
42.67912pt\begin{cases}\begin{tabular}[]{r l l
r}$\hat{r}$&$=\big{[}\sin(\theta)\cos(\phi),$&$\sin(\theta)\sin(\phi),$&$\cos(\theta)\big{]}$,\\\
$\hat{\theta}$&$=\big{[}\cos(\theta)\cos(\phi),$&$\cos(\theta)\sin(\phi),$&$\sin(\theta)\big{]}$,\\\
$\hat{\phi}$&$=\big{[}-\sin(\phi),$&$\cos(\phi),$&$0\big{]}$,\\\
$\hat{p}$&$=\big{[}\sin(\theta_{p})\cos(\phi_{p}),$&$\sin(\theta_{p})\sin(\phi_{p}),$&$\cos(\theta_{p})\big{]}$.\end{tabular}\end{cases}$
(2.28, 3.3 r)

The Fresnel monochromatic model equation 3.69 approaches the plane-wave
monochromatic model equation 3.37 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$). However, it is also important to
remember that the plane-wave monochromatic model itself converges to zero in
both of the natural plane-wave limits, as discussed in Section 3.5.2. Thus
both of these effects result in zero timing residual for these two precise
alignments, and hence the natural plane-wave limit will not be of much
interest to us going forward in this study.

Additionally, in this regime the degeneracy of the source chirp mass
$\mathcal{M}$ and distance $R$ parameters is also broken, so in principle a
loud enough observed source with a significant Fresnel number with the pulsar
should allow for independent measurements of both of these source parameters.

##### Asymptotic Model

Also of interest here we note the result when considering the asymptotic
expansion of the Fresnel integrals. Both $\eta_{2},\ \eta_{1}\gg 1$, so long
as we are not looking at a source that is nearly perfectly aligned or anti-
aligned in the plane-wave limit (where $\hat{r}\cdot\hat{p}\approx\pm 1$). In
this case, after replacing the Fresnel integrals with their asymptotic
expansions given in equation 3.66, the result simplifies to:

$\displaystyle\mathrm{Res}(t)\approx\frac{F^{\textsc{A}}}{4\omega_{0}}\Bigg{[}$
$\displaystyle\Bigg{\\{}\frac{1}{1+\left(1+\hat{r}\cdot\hat{p}\right)\left(\frac{L}{R}\right)}\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)$
$\displaystyle+\Bigg{\\{}\frac{-1}{1+\left(1+\hat{r}\cdot\hat{p}\right)\left(\frac{L}{R}\right)}\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{]},$
$\displaystyle\approx\frac{F^{\textsc{A}}}{4\omega_{0}}\Bigg{[}$
$\displaystyle\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{.}$
$\displaystyle\hskip
14.22636pt+\Bigg{.}\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{]}.$
(3.72)

In the final line we dropped the $\mathcal{O}\left(\frac{L}{R}\right)$ factors
in the amplitude that came from the asymptotic expansion, and now we see the
Fresnel formalism solution more closely resembles that of the plane-wave
solution, as they both share the same dominant amplitude term.

##### Heuristic Model

Our careful study and derivation of the timing residual in both the plane-wave
and Fresnel regimes reveals some interesting insights into the underlying
physics in these models. In the previous Sections 3.5 and 3.6 we have solved
for the timing residual model including higher order terms in the expansion of
the retarded times.

In the limit that the source is further away, we know that the Fresnel
formalism should reduce to the plane-wave formalism. And looking at the plane-
wave solution of the timing residual in equation 3.37 we see that it has two
terms which are evaluated at the endpoints of our photon’s journey. Evaluating
equation 3.28 at the start of the photon’s motion $t_{\mathrm{obs}}-L/c$ and
the end of the photon’s motion $t_{\mathrm{obs}}$ ultimately leads us to the
two times that the phase terms are evaluated at in equation 3.38, namely $t$
and $t-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}$ (given we choose
$t_{0}=-\frac{R}{c}$).

Referring again to Figure 3.4, we know the biggest potential difference
between the Fresnel and plane-wave regimes will come when evaluating the
retarded time at the pulsar’s position, since that is at the largest spatial
separation along the Earth-pulsar baseline. So these two observations lead us
to make the following guess: perhaps we can model the timing residual by using
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.

To try this, first note that, a priori, if we write the following expression:

$\displaystyle\overline{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{L}{c}-\frac{R}{c}\left[1-\left(\hat{r}\cdot\hat{p}\right)\left(\frac{L}{R}\right)+\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\left(\frac{L}{R}\right)^{2}+\ldots\right],$
$\displaystyle=t\quad-\frac{R}{c}\quad-\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}\quad-\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{R}\quad+\ldots.$
(3.73)

that the first three terms in this expression, when later combined with our
model evaluated at $t_{0}=-\frac{R}{c}$, will recover the retarded time factor
for the pulsar term in the plane-wave expansion equation 3.38. And further
given what we have established previously it suggests that this fourth term is
the Fresnel regime correction for the pulsar term.

These ideas lead us on the basis of physical interpretation of this problem
only, to propose the following model for the monochromatic regime timing
residual of a pulsar which includes higher order corrections of the effect of
the gravitational wavefront’s curvature coming from the retarded time:

$\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{(}\overline{\Theta}_{P}-\frac{\pi}{4}\Big{)}\bigg{]},$
(3.74) $\displaystyle\mathrm{for}\quad\textsc{A}\in[+,\times],$
$\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,\\\
\overline{\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}$
(3.75)

The overbar notation used here and below is simply to distinguish that, unlike
the previous timing residual models, this has not been analytically derived,
but rather proposed.

This proposed model has the required property that it reduces to the plane-
wave formalism in the limits given in assumption 1 in Section 3.6.1. It takes
the general form of the plane-wave monochromatic formalism (equation 3.37) but
with the corrections motivated by the Fresnel regime study. It is critical to
understand that this formula is not a mathematically derived formula for the
timing residual. It is merely meant to be a model which attempts to accurate
encapsulate the basic physics principles behind the true timing residual model
motivated from our previous sections. As we showed, the real Fresnel regime
solution equation 3.69 is more complicated. But from a practical standpoint,
we found that this formula is much easier to study numerically in our Fisher
analysis in Chapter 4.

As a model we can still test it to see how accurately it matches the known
analytic solution equation 3.69, and it turns out, both the asymptotic and the
heuristic models do exceptionally well at predicting the behavior of the
timing residual. Figure 3.8 shows an example of the difference between the
analytic and asymptotic formulae, and the analytic and the heuristic formulae.
For this example, the difference in models is on the order of femtoseconds,
which is well below the precision with which we can currently measure timing
residuals. But perhaps even more importantly, we see that these two
approximation models, equations 3.72 and 3.74 match each other, which further
validates our heuristic understanding of the underlying physics affecting the
timing residual model.

Figure 3.8: An example of the difference between the monochromatic timing
residual models in the Fresnel formalism. Here we take the difference between
the analytically derived model equation 3.69, and the two approximations we
discussed, the asymptotic model (equation 3.72) and the heuristic model
(equation 3.74). Both approximations accurately recover the analytic formula,
with the difference here on the order of femtoseconds. This is well below
current pulsar timing precision, which is around $100$ns (Arzoumanian et al.,
2018). For this example the following sets of parameters were used:
$\left\\{R,\theta,\phi,\iota,\psi,\theta_{0},\mathcal{M},\omega_{0}\right\\}=\left\\{100\
\mathrm{Mpc},\frac{\pi}{2}\ \mathrm{rad},0\ \mathrm{rad},0\ \mathrm{rad},0\
\mathrm{rad},0\ \mathrm{rad},10^{7}\ M_{\odot},10\ \mathrm{nHz}\right\\}$ and
$\left\\{L,\theta_{p},\phi_{p}\right\\}=\left\\{20\ \mathrm{kpc},0\
\mathrm{rad},0\ \mathrm{rad}\right\\}$, which gives the source a coalescence
time $\Delta\tau_{c}=2.0$ Gyr, and a Fresnel number of $F=1.31$ with this
pulsar.

Alternatively, we can decompose equation 3.74 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),$
(3.76)
$\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{(}\overline{\Theta}_{P}-\frac{\pi}{4}\Big{)}\bigg{]},$
$\displaystyle\hskip 98.16191pt=A_{c(\mathrm{IIA,correction})}\
F^{\textsc{A}}\ s_{P,\textsc{A}}$ (3.77) $\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}$
(3.78) $\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}$
(3.79)

Here $\mathrm{Res}_{(\mathrm{IA})}(t)$ is the plane-wave monochromatic model
from Section 3.5.2, and equation 3.77 is the correction to the plane-wave
regime due to the wavefront curvature from the Fresnel terms. Note that
$\Theta_{P}$ in equations 3.77 and 3.79 comes from the IA model equation 3.38.
Using the same trick as in Section 3.5.2 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|.$
(3.80)

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 3.25), $F$ and $\rho$ are two useful metrics in evaluating
the significance of the curvature corrections to our models.

####  3.6.2 Fresnel, Frequency Evolution - Conjecture (IIB)

Given the understanding of the previous timing regimes, we now propose the
most general timing residual model to date, one which includes the Fresnel
corrections, allows for frequency evolution, and still recovers all of the
previous three regimes in the appropriate limits. This is not analytically
derived - a complete and careful derivation along the lines of what we have
done previously is still needed here. However, our proposed conjecture here is
motivated by the reasoning explained from our insights of the heuristic model
in Section 3.6.1. This model makes the same assumptions as in the plane-wave
frequency evolution regime (Section 3.5.3) and the Fresnel monochromatic
regime (Section 3.6.1).

In the monochromatic formalism we generalized the plane-wave solution by
modifying the retarded time in the pulsar term to include the Fresnel
correction term, which encoded information about the curvature of the wave. We
propose doing the same thing with the frequency evolution formalism. We begin
with the plane-wave frequency evolution regime timing residual equation 3.47,
and we modify the pulsar term’s retarded time using equation 3.73 as discussed
in the Fresnel monochromatic case. This will in turn modify both the orbital
frequency of the pulsar term $\overline{\omega}_{0P}$ as well as the phase
$\overline{\theta}_{0P}$ from equation 3.48. 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{(}\overline{\omega}_{0P},\overline{\Theta}_{P}-\frac{\pi}{4}\Big{)}}{\overline{\omega}_{0P}}\right],$
(3.81) $\displaystyle\mathrm{for}\quad\textsc{A}\in[+,\times],$
$\displaystyle\underset{\left(t_{0}\ =\
-\frac{R}{c}\right)}{\mathrm{where}}\quad\begin{cases}\omega_{0E}&\equiv\omega_{0},\\\
\overline{\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},\\\
\overline{\theta}_{0P}&=\theta_{0}+\theta_{c}\left(1-\left[1+\frac{\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}+\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{L}{c}\frac{L}{R}}{\Delta\tau_{c}}\right]^{5/8}\right),\\\
\Theta_{E}&=\theta_{0E}+\omega_{0E}t,\\\
\overline{\Theta}_{P}&=\overline{\theta}_{0P}+\overline{\omega}_{0P}t,\end{cases}$
(3.82) $\displaystyle\hskip
45.52458pt\begin{cases}\Delta\tau_{c}\equiv\frac{5}{256}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\frac{1}{\omega_{0}^{8/3}},\\\
\theta_{c}\equiv\frac{8}{5}\Delta\tau_{c}\omega_{0},\end{cases}$ (2.36, 2.39
r) $\displaystyle\hskip
45.52458pt\begin{cases}h_{+}\equiv-h(\omega)\cos\big{(}2\Theta\big{)},\\\
h_{\times}\equiv-h(\omega)\sin\big{(}2\Theta\big{)},\end{cases}$ (2.22 r)
$\displaystyle\hskip
45.52458pth(\omega)\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega^{2/3},$
(2.23 r) $\displaystyle\hskip 45.52458pt\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{\left(\hat{p}\cdot\hat{\theta}\right)^{2}-\left(\hat{p}\cdot\hat{\phi}\right)^{2}}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\\\
$f^{\times}$&$\equiv\frac{\hat{p}^{i}\hat{p}^{j}e^{\hat{r}\times}_{ij}}{\left(1-\hat{r}\cdot\hat{p}\right)}$&$=\frac{2\left(\hat{p}\cdot\hat{\theta}\right)\left(\hat{p}\cdot\hat{\phi}\right)}{\left(1-\hat{r}\cdot\hat{p}\right)}$
,\end{tabular}\end{cases}$ (3.33 r) $\displaystyle\hskip
45.52458pt\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{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)f^{+}+\sin(2\psi)f^{\times}\right]$
,\\\
$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.34 r) $\displaystyle\hskip
45.52458pt\begin{cases}\begin{tabular}[]{r l l
r}$\hat{r}$&$=\big{[}\sin(\theta)\cos(\phi),$&$\sin(\theta)\sin(\phi),$&$\cos(\theta)\big{]}$,\\\
$\hat{\theta}$&$=\big{[}\cos(\theta)\cos(\phi),$&$\cos(\theta)\sin(\phi),$&$\sin(\theta)\big{]}$,\\\
$\hat{\phi}$&$=\big{[}-\sin(\phi),$&$\cos(\phi),$&$0\big{]}$,\\\
$\hat{p}$&$=\big{[}\sin(\theta_{p})\cos(\phi_{p}),$&$\sin(\theta_{p})\sin(\phi_{p}),$&$\cos(\theta_{p})\big{]}$.\end{tabular}\end{cases}$
(2.28, 3.3 r)

This model behaves as expected in the previously established limits,
recovering all previous regimes.

And just like we did with the heuristic Fresnel monochromatic model in Section
3.6.1, 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),$
(3.83)
$\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{(}\overline{\omega}_{0P},\overline{\Theta}_{P}-\frac{\pi}{4}\Big{)}}{\overline{\omega}_{0P}}\right].$
(3.84)

Here $\mathrm{Res}_{(\mathrm{IB})}(t)$ is the plane-wave frequency evolution
model from Section 3.5.3, and equation 3.84 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 3.84 come from the IB model
equation 3.48.

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 2.23), 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.7 Regime & Model Comparisons

Figure 3.10 shows four examples that highlight the differences in these timing
model regimes, and Figure 3.12 emphasizes just how different the predicted
timing results can be depending on which model you use. It is important to
keep in mind the limitations and assumptions under which each of these timing
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).

This has lead us to consider the following question, “by how much will our
parameter estimations be affected when using the wrong model regime for a
particular source?” Consider, for example, a particular source has both
significant Fresnel and frequency evolution effects that contribute to the
timing residual it induces on our PTA. It would be important to know if using
the IB model (which doesn’t capture those Fresnel effects) for parameter
estimation would result in estimates of the wrong source parameters, and by
how much they are wrong. At the time of completing my PhD work and finishing
this dissertation I had begun a preliminary investigation into this question,
so this may result in the topic of future studies beyond the work presented
here.

The four model regimes can be further subdivided into different categories
given the parameter dependencies of the coalesence time and the Fresnel number
(equations 2.36 and 3.25), as shown in Table 3.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 Figure 3.5). Figure 3.9 shows contour plots of
both the coalescence time $\Delta\tau_{c}$ and Fresnel number $F$ in terms of
their parameter dependencies. In general these can help give us a sense of
where in parameter space we can expect a transition from one of our timing
regimes to another.

| 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)$
Table 3.1: The 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. Figure 3.9: Contour plots of the
coalescence time $\Delta\tau_{c}$ and Fresnel number $F$. We expect the
frequency evolution regime to become significant for smaller coalescence
times, where the chirp mass $\mathcal{M}$ and orbital frequency $\omega_{0}$
are sufficiently high. Similarly the Fresnel regime should become significant
for higher Fresnel numbers. Here as an example we fix the source distance $R$
of a gigaparsec and look at contours in the pulsar distance $L$ and orbital
frequency.

Consequently, the orbital frequency parameter has the most direct influence on
the regimes - increasing the orbital frequency pushes the models towards both
the Fresnel and frequency evolution regimes, and vice versa. When studying
these pulsar timing models, special consideration should be taken in this
regard, because studying systems at higher frequencies not only suggests that
the physical effects due to frequency evolution could become important, but
that also the effects of curvature as well. These physical effects will change
the timing model simultaneously, which 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.

Furthermore, it also gives us motivation to look for and time pulsars at
greater distances. Since the Fresnel number scales as $L^{2}$, increasing the
pulsar distance is the easiest way to push the timing 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 in Chapter 6.

Figure 3.10: 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\\}$.
Figure 3.11:  Continued
(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\\}$.
Figure 3.12: Here we compute the predicted timing residual for Source 2 (Table
LABEL:tab:_source_parameters_reference) at $R=100$ Mpc and the pulsar
J1747-4036 (given the values in Table LABEL:tab:_pta_reference, pulsar #40)
using the four pulsar timing regime models. This source has a coalescence time
of $12.793$ kyr, and the Fresnel number of this source and pulsar is
$F=0.837$. This is meant to show just how different the results can be in the
predicted timing models for a given source and pulsar across each regime. It
is important to consider the strength of the frequency evolution in the signal
and how large the Fresnel number is, in order to decide the accuracy of each
model regime’s predicted results, since certain regimes have physical
limitations. For example, the monochromatic models are insensitive to any
frequency evolution over the pulsar-Earth baseline that the photon travels,
which may otherwise produce a beat frequency in the signal (as seen in the two
frequency evolution regimes here). And the plane-wave regimes are not
sensitive to Fresnel corrections in the phase, which may otherwise both shift
the phase, produce different amplitudes, and further change the underlying
beat frequency.

###  3.8 Regime Search Parameters & Degeneracies

Each of the model regimes which we have discussed allow for the measurement of
different sets of source parameters experimentally depending on which
parameters are degenerate in the model. However certain parameters are
degenerate or highly covariant with others in certain regimes, which has
motivated us to parametrize the model in specific ways. For reference, Table
3.2 indicates the parameters in each of our model regimes.

| 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\\}$
Table 3.2: Our base parametrization choice 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
3.85). In general, frequency evolution allows the direct measurement of chirp
mass $\mathcal{M}$ from the frequency $\omega(t)$ and phase $\Theta(t)$ terms
in the timing residual models, and Fresnel corrections allow for direct
measurement of source distance $R$ from these frequency and phase terms.
Fisher matrices with alternative parametrization choices were computed through
the Jacobian transformation in equation 4.31.

Starting with the most basic model, the plane-wave monochromatic regime, the
source coordinate distance $R$ and chirp mass $\mathcal{M}$ are fully
degenerate, appearing only in the amplitude of the timing residual (equation
3.37). We use equations 3.37, 3.39, and 3.40 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}},$
(3.85)

which is the “Earth term timing residual amplitude,” a similar quantity to
characteristic amplitudes (equations 3.42 and 3.78). 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
3.38), 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 models 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. To clarify, this does mean
that in the two frequency evolution regimes IB and IIB, we wrote the amplitude
of the pulsar term as the combination of parameters
$A_{P,\mathrm{res}}=A_{E,\mathrm{res}}\left(\frac{\omega_{0}}{\omega_{0P}}\right)^{1/3}$
for the plane-wave regime, and
$A_{P,\mathrm{res}}=A_{E,\mathrm{res}}\left(\frac{\omega_{0}}{\overline{\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.

From this we see one novel result of this study is that the Fresnel monochrome
regime allows for the recovery of the coordinate distance even if the
gravitational wave is effectively monochrome. Furthermore, this study shows
that going into the full Fresnel frequency evolution regime allows separate
measurements of both the coordinate distance and the chirp mass from the
frequency and phase terms in the timing residual. These key points are
investigated in detail in Section 6.4.

##  Chapter 4 Fisher Matrix Analysis

The Fisher matrix provides a powerful tool for analyzing the predictive
capabilities of a Bayesian-based model of some physical process. At its core,
the beauty and purpose of the Fisher matrix formalism is that it gives a
“best-case-scenario” of how well the model parameters will be measured and
constrained by the experiment itself, given the inherent uncertainties in
making the observations, before the experiment is even performed. An excellent
motivational summary of the Fisher matrix formalism is:

> “The whole point of the Fisher matrix formalism is to predict how well the
> experiment will be able to constrain the model parameters, before doing the
> experiment and in fact without even simulating the experiment in any detail.
> We can then forecast the results of different experiments and look at
> tradeoffs such as precision versus cost. In other words, we can engage in
> experimental design. …The beauty of the Fisher matrix approach is that there
> is a simple prescription for setting up the Fisher matrix knowing only your
> model and your measurement uncertainties; and that under certain standard
> assumptions, the Fisher matrix is the inverse of the covariance matrix. So
> all you have to do is set up the Fisher matrix then invert it to obtain the
> covariance matrix (that is, the uncertainties on your model parameters)… the
> inverse of the Fisher matrix is the best you can possibly do given the
> information content of your experiment.” — David Wittman (Wittman)
>

Here we will develop the derivation of the Fisher matrix itself from Bayesian
principles, and apply it to our pulsar timing model.

###  4.1 High Statistical Limit Applied to Bayes’ Theorem

From Bayes’ Theorem we know:

$\underbrace{p\left(\vec{\theta}\middle|\vec{X}\right)}_{\mathrm{``posterior"}}=\frac{\overbrace{p\left(\vec{X}\middle|\vec{\theta}\right)}^{\mathrm{``likelihood"}}\overbrace{p\left(\vec{\theta}\right)}^{\mathrm{``prior"}}}{\underbrace{p\left(\vec{X}\right)}_{\mathrm{``evidence"}}},$
(4.1)

where $\vec{\theta}$ are the model parameters and $\vec{X}$ are the observed
random variables (i.e. the data/evidences). Conceptually, the “evidence” is
the total probability of the data $\vec{X}$ regardless of any other
information. The “prior” represents what we know about the parameters
$\vec{\theta}$ before we see the data, and the “likelihood” represents what we
have learned from the observations in our experiment. In this work we will use
the notation
$p\left(\vec{X}\middle|\vec{\theta}\right)\equiv\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)$
to help more quickly differentiate between posterior and the likelihood.

Often when applying Bayesian statistics to a model the evidence is something
that is unknown. However, since it is only a normalizing constant (which
normalizes integration of the posterior to $1$), the actual shape and
“information” contained in the posterior distribution comes entirely from the
likelihood and the prior. This is why it is often written as shorthand that:

$\mathrm{posterior}\propto\mathrm{likelihood}\times\mathrm{prior}.$

Once we have the posterior distribution we can normalize it numerically in the
end if we desire, hence often we don’t have a direct need to compute
$p\left(\vec{X}\right)$ ahead of time.

Typically we find ourselves performing an experiment where we record
observations of some quantity. This quantity we try to predict ahead of time
by mathematically describing it using a model which may be dependent on a
number of parameters, our $\vec{\theta}$. However, we often also do not know
the true values of the parameters governing the thing we are observing in our
experiment. So the problem, then, is to figure out how to infer and estimate
these parameters $\vec{\theta}$ given a set of observation data $\vec{X}$ from
our experiment, and this can be done through a Fisher matrix analysis.

Now, the key behind the Fisher matrix and this problem is rooted in how we
would expect our experiment to behave if we were we able to repeat the
experiment over and over for many realizations (tending towards infinity).
From Bayes’ Theorem we see that the posterior and the likelihood are directly
proportional to each other, and in this so called “high statistical limit” or
“asymptotic limit” of repeating our experiment many times we expect the
posterior will tend to depend less and less on the prior information.

Moreover, often the tendency is that asymptotically the posterior will become
Gaussian distributed (Nikulin, 2011a). These ideas come from the Bernstein-von
Mises theorem, and require certain assumptions such as Cromwell’s rule (Glen,
2018), which says that the prior probability on a parameter should never be
fully $0$ or $1$. The rationale on this assumption is that if the parameter is
fixed absolutely at always being observed or never being observed, then the
posterior can’t be influenced by the actual data and observations taken into
the model described by the likelihood.

The Fisher matrix is by no means a perfect estimator, and this type of
analysis does have its limitations which must be considered. One is the role
of the previously mentioned prior information. It is not universally true that
in the high asymptotic limit the prior information will have no influence on
the resulting posterior. In fact, the prior can be quite informative. One key
idea is that in the region of interest of parameter space we are focusing on
(near the true values), the prior information should remain approximately
constant/unchanging (Vallisneri, 2008). This is because, as will be shown in
the next sections, the Fisher matrix is effectively a second-order Taylor
expansion of the likelihood function which is being used to approximate the
posterior. Because this quadratic approximation describes the local
“curvature” in parameter space about the true parameter values, the Fisher
matrix should remain decent so long as the curvature of the prior is less than
the curvature of the likelihood.

So with these ideas in mind we are effectively looking for a maximum
likelihood estimation (MLE), that is a way of predicting the values of
$\vec{\theta}$ in parameter space which maximize our likelihood. But something
important to always remind ourselves of (which we emphasize explicitly here)
is that our observations $\vec{X}$ are random variables (R.V.s). Now what we
expect is that the observations contain some “true” value produced by nature
(the thing we are really trying to observe), combined with some inherent noise
(also produced by nature… thanks, nature). In other words
$\vec{X}=\vec{X}_{\mathrm{true}}+\vec{n}$, where $\vec{n}$ is the noise in
each measurement. This is what makes our observation an R.V. - here $\vec{X}$
is an R.V. because $\vec{n}$ is an R.V. What this means is that when we
perform an experiment, even if we were using the correct model with the true
parameters for the thing being observed, $\vec{\theta}_{\mathrm{true}}$, the
likelihood won’t necessarily be at its maximum value. If the noise is exactly
zero, then the likelihood is indeed maximized, but every time we repeat the
experiment there will be different random noise added into the observation (as
shown in Figure 4.1(a)). The random noise spreads the likelihood out into a
distribution in terms of the noise. If we happen to pick the wrong set of
parameters, some $\bar{\theta}\neq\theta_{\mathrm{true}}$, then the
distribution will shift and peak for some non-zero value of the noise, as
shown in Figure 4.1(b). This is why we must treat the problem statistically,
and consider what is asymptotically the most likely set of parameters for our
model, and what is our statistical uncertainty of their measurements?

(a) If we happen to choose the correct parameter parameter $\theta$ in our
model, then the likelihood distribution will peak for zero noise. Every time
we repeat the experiment the noise will be different, so we will compute a
different likelihood, as what is being exemplified by the first five
iterations of our experiment here. Eventually, if we repeat this experiment
many many times, we will see the full distribution’s shape.

(b) If we choose a different parameter
$\bar{\theta}\neq\theta_{\mathrm{true}}$, then the distribution will peak for
some non-zero value of the noise.

Figure 4.1: The likelihood function $\mathscr{L}\left(X\middle|\theta\right)$
gives us a “goodness-of-fit” distribution of our data $X$ for a given set of
parameter $\theta$ in our model. The data is a random variable since it is a
true “signal” plus some random noise $n$.

###  4.2 Conceptually Understanding the Fisher Information

In order to focus on motivating all of the logic and ideas which go into
understanding and deriving the Fisher matrix (Albrecht et al., 2009), we
consider here first on a very simple example where we imagine that we perform
an experiment which collects one datum $X_{\mathrm{obs}}$ with Gaussian
uncertainty $\sigma$, and our functional model for this observation only
depends on one parameter $X(\theta)$. The ideas here will extend to a fully
multidimensional data and parameter-space, but for the moment we keep it
simple so as to focus entirely on the logic.

We begin by defining a Chi-squared test as one way of trying to quantify a
sense of the “goodness-of-fit” of our model (given its parameter) to the
observed data:

$\chi^{2}\equiv\frac{\tilde{X}^{2}}{\sigma^{2}}=\frac{\Big{(}X(\theta)-X_{\mathrm{obs}}\Big{)}^{2}}{\sigma^{2}},$
(4.2)

where $X_{\mathrm{obs}}=X_{\mathrm{true}}+n$ with inherent noise $n$. If this
model and its parameters is correct and really does describe the physics of
our experiment, then the likelihood should be a Gaussian, which can be written
in terms of our Chi-squared expression:

$\mathscr{L}\left(X\middle|\theta\right)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\chi^{2}\right].$
(4.3)

Now as we discussed in Section 4.1, our goal is to infer the parameters
$\theta$ from our observation $X_{\mathrm{obs}}$. Due to the inherent noise in
our observations, $X_{\mathrm{obs}}$ is an R.V., which means that simply
minimizing the Chi-squared value, i.e. maximizing our likelihood, for this one
realization of our experiment won’t necessarily give us the “true” parameters.
This would only happen if the noise were exactly zero, otherwise the noise
pushes our observation away from the true $X_{\mathrm{true}}$ (see Figure
4.1(a)).

From Bayes’ theorem equation 4.1 we know the posterior is proportional to the
likelihood,
$p\left(\theta\middle|X\right)\propto\mathscr{L}\left(X\middle|\theta\right)$.
Additionally we may or may not know prior information about our parameter
which we could supply to help inform our experiment. As mentioned already in
Section 4.1, in the approximation we are about to make we are sort of assuming
that the prior is not very informative over the small region of parameter
space that we are going to focus in on. Therefore let’s focus on the high
statistical limit of just our Chi-squared value in the likelihood - what is
the value we expect it to be as we repeat the experiment ad infinitum?
Mathematically, we are going to solve for the expectation value
$\left\langle\chi^{2}\right\rangle$ in order to ask what is the expected
likelihood.

In addition to this, however, we are also interested in focusing our attention
on what is the expectation of our Chi-squared and likelihood near the true
parameter $\theta_{\mathrm{true}}$, the true value of $\theta$ which we are
trying to infer from our data. Hence first we Taylor expand Chi-squared for
small $\Delta\theta=\theta-\theta_{\mathrm{true}}$ about the true parameter
as:

$\chi^{2}=\quad\chi^{2}\bigg{\rvert}_{\theta_{\mathrm{true}}}\quad+\quad\frac{d}{d\theta}\left[\chi^{2}\right]\bigg{\rvert}_{\theta_{\mathrm{true}}}\Delta\theta\quad+\quad\frac{1}{2}\frac{d^{2}}{d\theta^{2}}\left[\chi^{2}\right]\bigg{\rvert}_{\theta_{\mathrm{true}}}\Delta\theta^{2}\quad+\quad\mathcal{O}\left(\Delta\theta^{3}\right).$
(4.4)

Next we want to take the expectation value of this. As a reminder, the
expectation value is over our R.V., which here is $X_{\mathrm{obs}}$ or $n$
(or $\tilde{X}$). Since we are treating these continuously, the expectation
value weights the argument against the probability density we have, which in
this case is our likelihood. Furthermore, we are going to take the expectation
value evaluated at the true parameter $\theta_{\mathrm{true}}$. This is
because again, we want to focus our attention on what happens to Chi-squared
near the true parameter over an (ideally) infinite number of experiment
iterations (see Figure 4.1(b)). So let the following notation mean:

$\displaystyle\langle\ \cdot\ \rangle=\int\limits^{\infty}_{-\infty}\left[\
\cdot\ \right]\ \
\mathscr{L}\left(X\middle|\theta_{\mathrm{true}}\right)dX_{\mathrm{obs}}$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\ \cdot\ \right]\ \
\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\chi\left(X_{\mathrm{obs}},\theta_{\mathrm{true}}\right)^{2}\right]dX_{\mathrm{obs}},$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\ \cdot\ \right]\ \
\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\chi\left(n,\theta_{\mathrm{true}}\right)^{2}\right]dn,$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\ \cdot\ \right]\ \
\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\chi\left(\tilde{X}\left(\theta_{\mathrm{true}}\right)\right)^{2}\right]d\\!\\!\left[\tilde{X}\left(\theta_{\mathrm{true}}\right)\right],$
(4.5)

where the three integral expressions here are all equivalent since they are
all directly related to each other through a simple change of variables given
their definitions. So now we can take the expectation value of equation 4.4
and write:

$\left\langle\chi^{2}\right\rangle=\ \
\left\langle\chi^{2}\bigg{\rvert}_{\theta_{\mathrm{true}}}\right\rangle\ \ +\
\
\left\langle\frac{d}{d\theta}\left[\chi^{2}\right]\bigg{\rvert}_{\theta_{\mathrm{true}}}\right\rangle\Delta\theta\
\ +\ \
\frac{1}{2}\left\langle\frac{d^{2}}{d\theta^{2}}\left[\chi^{2}\right]\bigg{\rvert}_{\theta_{\mathrm{true}}}\right\rangle\Delta\theta^{2}\
\ +\ \ \mathcal{O}\left(\Delta\theta^{3}\right).$ (4.6)

Working this out piece-by-piece we have for the first two terms:

$\displaystyle\left\langle\chi^{2}\bigg{\rvert}_{\theta_{\mathrm{true}}}\right\rangle\hskip
24.18501pt$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\frac{\tilde{X}\left(\theta_{\mathrm{true}}\right)^{2}}{\sigma^{2}}\right]\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\chi\left(\theta_{\mathrm{true}}\right)^{2}\right]d\\!\\!\left[\tilde{X}\left(\theta_{\mathrm{true}}\right)\right],$
$\displaystyle=\frac{1}{\sigma^{2}}\frac{1}{\sqrt{2\pi}\sigma}\int\limits^{\infty}_{-\infty}\tilde{X}\left(\theta_{\mathrm{true}}\right)^{2}\exp\left[-\frac{1}{2}\frac{\tilde{X}\left(\theta_{\mathrm{true}}\right)^{2}}{\sigma^{2}}\right]d\\!\\!\left[\tilde{X}\left(\theta_{\mathrm{true}}\right)\right],$
$\displaystyle=\frac{1}{\sigma^{2}}\frac{1}{\sqrt{2\pi}\sigma}\sqrt{2\pi}\sigma^{3}\equiv
1,$ (4.7)
$\displaystyle\left\langle\frac{d}{d\theta}\left[\chi^{2}\right]\bigg{\rvert}_{\theta_{\mathrm{true}}}\right\rangle$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\frac{2}{\sigma^{2}}\tilde{X}\left(\theta_{\mathrm{true}}\right)\frac{\partial}{\partial\theta}\left[\tilde{X}\right]\bigg{\rvert}_{\theta_{\mathrm{true}}}\right]\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{1}{2}\chi\left(\theta_{\mathrm{true}}\right)^{2}\right]d\\!\\!\left[\tilde{X}\left(\theta_{\mathrm{true}}\right)\right],$
$\displaystyle=\frac{2}{\sigma^{2}}\frac{\partial}{\partial\theta}\Big{[}X\left(\theta_{\mathrm{true}}\right)\Big{]}\frac{1}{\sqrt{2\pi}\sigma}\int\limits^{\infty}_{-\infty}\tilde{X}\left(\theta_{\mathrm{true}}\right)\exp\left[-\frac{1}{2}\frac{\tilde{X}\left(\theta_{\mathrm{true}}\right)^{2}}{\sigma^{2}}\right]d\\!\\!\left[\tilde{X}\left(\theta_{\mathrm{true}}\right)\right],$
$\displaystyle\equiv 0.$ (4.8)

The result of this shows us that the first non-zero terms of our expectation
of Chi-squared near the true parameter value is:

$\displaystyle\left\langle\chi^{2}\right\rangle\quad$ $\displaystyle=\quad
1\quad+\quad
F\Delta\theta^{2}\quad+\quad\mathcal{O}\left(\Delta\theta^{3}\right),$
$\displaystyle\mathrm{where}\quad
F=F\left(\theta_{\mathrm{true}}\right)\equiv\frac{1}{2}\left\langle\frac{d^{2}}{d\theta^{2}}\left[\chi^{2}\right]\bigg{\rvert}_{\theta_{\mathrm{true}}}\right\rangle.$
(4.9)

We can make a conceptual understanding of this in the following way. When
there is no noise, Chi-squared is zero at $\theta_{\mathrm{true}}$. However,
the average or expected value of Chi-squared evaluated at the true parameter
$\left\langle\chi^{2}\big{\rvert}_{\theta_{\mathrm{true}}}\right\rangle$ when
weighted against the likelihood probabilities of $\chi^{2}$ which are also
evaluated at the true parameter (this is important) is $1$. Next looking
locally around $\theta_{\mathrm{true}}$ for small deviations
$\mathcal{O}\left(\Delta\theta\right)$, we find the that linear term
$\Delta\theta$ vanishes, which tells us we are at an extremum here in our
parameter space. However, the quadratic $\Delta\theta^{2}$ term does not
vanish and can be computed, and since it is quadratic this quantity tells us
the local curvature around our extremum. We call this quantity the “Fisher
information,” denoted by $F$, and it is a value which depends on the true
parameter $\theta_{\mathrm{true}}$.

So we have used the likelihood evaluated at the true parameter
$\theta_{\mathrm{true}}$ to calculate the expectation value of our model’s
Chi-squared for small deviations off of the true parameter. This is our high
statistical limit, and now we take that and use it to give us our new high-
statistical limit likelihood, and hence also our high statistical limit
approximation of our posterior:

$\displaystyle
p\left(\theta\middle|X\right)\propto\mathscr{L}\left(X\middle|\theta\right)=\textsf{C}\exp\left[-\frac{1}{2}\left\langle\chi^{2}\right\rangle\right]$
$\displaystyle=\textsf{C}\exp\left[-\frac{1}{2}\Big{(}1+F\Delta\theta^{2}\Big{)}\right],$
$\displaystyle=\bar{\textsf{C}}\exp\left[-\frac{1}{2}F\Delta\theta^{2}\right],$
(4.10)

where notice that C is our new normalization constant for this probability
density (here we just haven’t computed it explicity), and $\bar{\textsf{C}}$
has simply absorbed the additional factor of $\exp\left[-\frac{1}{2}\right]$
which is just a constant. Since this is a Gaussian still, the new
normalization constant would be $\bar{\textsf{C}}=\sqrt{\frac{F}{2\pi}}$. Note
that in this asymptotic approximation, we have dropped
$\mathcal{O}\left(\Delta\theta^{3}\right)$ terms.

Consider what we have done here. The posterior is a distribution on the
parameters (given the data), while the likelihood is a distribution of the
data (given the parameters). However, through Bayes’ theorem and a high
statistical limit, averaging over realizations of the data itself, we now view
the likelihood as a function of the parameters. The result is a Gaussian
dependent on this new Fisher information, which is itself a function of the
true parameter values. So the Fisher analysis takes us from a Gaussian
likelihood in data-space to a Gaussian posterior in parameter-space. Looking
at equation 4.10, the inverse of that Fisher information is the variance on
$\Delta\theta$, that is $F^{-1}=\sigma^{2}$. So inverting the Fisher
information here would tell us how large our uncertainty is on the parameter
$\theta$, centered on the true value of $\theta_{\mathrm{true}}$.

The final crucial key to this comes from the Cramér-Rao theorem (Nikulin,
2011b), which states that the variance produced by some unbiased estimator
will be at least as high as the inverse of this Fisher information. Or in
other words, the Fisher information gives the “best-case-scenario” for how
well we can measure the parameter in our model, which is what makes it such an
important quantity.

To get equation 4.10 we took the high-statistical limit of Chi-squared about
our true parameter, then used that to define our new likelihood in the high
statistical limit. Note though we can cast the definition of the Fisher
information in a more general form rather than the form given in equation 4.9
as follows. Start with the likelihood function equation 4.3. Then take the
negative second derivative of the natural log of that quantity:

$\displaystyle-\frac{\partial^{2}}{\partial\theta^{2}}\ln{\mathscr{L}\left(X\middle|\theta\right)}$
$\displaystyle=-\frac{\partial^{2}}{\partial\theta^{2}}\left[\ln{\frac{1}{\sqrt{2\pi}\sigma}}-\frac{1}{2}\chi^{2}\right],$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}}{\partial\theta^{2}}\chi^{2}.$

Next evaluate this expression at $\theta=\theta_{\mathrm{true}}$, take the
expectation value of both sides to integrate over the random variable, and
drop $\mathcal{O}\left(\Delta\theta^{3}\right)$ terms:

$\displaystyle-\left\langle\frac{\partial^{2}}{\partial\theta^{2}}\ln{\mathscr{L}\left(X\middle|\theta_{\mathrm{true}}\right)}\right\rangle$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}}{\partial\theta^{2}}\left\langle\chi^{2}\right\rangle,$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}}{\partial\theta^{2}}\left[1+F\Delta\theta^{2}+\mathcal{O}\left(\Delta\theta^{3}\right)\right],$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}}{\partial\theta^{2}}\left[1+F\left(\theta-\theta_{\mathrm{true}}\right)^{2}+\mathcal{O}\left(\Delta\theta^{3}\right)\right],$
$\displaystyle\approx F.$

Hence we can see that, assuming again we are considering sufficiently small
deviations away from the true parameter
$\Delta\theta=\theta-\theta_{\mathrm{true}}$, in addition to equation 4.9, we
can say our Fisher information is:

$F=-\left\langle\frac{\partial^{2}}{\partial\theta^{2}}\ln{\mathscr{L}\left(X\middle|\theta_{\mathrm{true}}\right)}\right\rangle.$
(4.11)

Most sources say this is the definition of the Fisher information, which we
can choose it to be, but the logic has this built in assumption regarding
small $\Delta\theta$ which we want to point out from our logic and motivated
derivation here.

###  4.3 Generalizing from the Fisher Information to the Fisher Matrix

Now that we have the conceptual understanding of the ideas underlying the
Fisher information, here we present the fully generalized version. Before we
considered one datum $X_{\mathrm{obs}}$ and a functional model to describe it
which only depended on one parameter $\theta$. From our analysis, we found our
“Fisher information” was a scalar value that described the local curvature in
our 1-D parameter space around the true parameter.

However, now suppose we collect many data, represented by a $d$-dimensional
vector $\vec{X}_{\mathrm{obs}}$, each with Gaussian $\sigma$ uncertainties,
and construct a model for that data which depends on multiple parameters,
represented by a $k$-dimensional vector $\vec{\theta}$. In our work and
derivation here we are assuming that there are no covariances between observed
data measurements - while each data measurement may carry its own uncertainty,
it is not correlated with another data measurement. Finally, in our notation
here, we will use indices $(a,b,c)$ in our data space, and indices $(i,j)$ in
our parameter space. We will use both vector/matrix notation and Einstein
index notation (and implied summation conventions) in our work below.

As we did before, let’s define a multidimensional goodness-of-fit Chi-squared
quantity for this observed data set and our model:

$\displaystyle\chi^{2}$
$\displaystyle\equiv\vec{\tilde{X}}^{T}\bm{\Sigma}^{-1}\vec{\tilde{X}}=\Sigma^{-1}_{ab}\tilde{X}^{a}\tilde{X}^{b}=\frac{1}{\sigma^{2}_{c}}\left(\tilde{X}^{c}\right)^{2}=\frac{1}{\sigma^{2}_{c}}\bigg{(}X^{c}\left(\vec{\theta}\right)-X^{c}_{\mathrm{obs}}\bigg{)}^{2},$
(4.12)
$\displaystyle\mathrm{where}\quad\Sigma^{-1}_{ab}\equiv\frac{1}{\sigma^{2}_{c}}\delta^{c}_{ab},$
$\displaystyle\mathrm{and}\hskip
19.91684pt\delta^{c}_{ab}\equiv\begin{cases}1,\quad a=b=c\\\
0,\quad\text{else}\end{cases}$

Here $\bm{\Sigma}$ is our covariance matrix for the data observations,
however, as we noted already we are assuming here there are no covariances
between data measurements (this is really just a “variance matrix”). Note that
our definition of the delta notation here is just a clever way of expressing
that
$\bm{\Sigma}=\mathrm{diag}\left(\sigma^{2}_{1},\sigma^{2}_{2},\dots,\sigma^{2}_{d}\right)$
is just a diagonal matrix where each element on the diagonal can be a
different value. We still note that
$\vec{X}_{\mathrm{obs}}=\vec{X}_{\mathrm{true}}+\vec{n}$, where the noise
makes this a multidimensional R.V. Again, if our model is correct and does
truly describe the physics of our experiment, then we expect the likelihood
should be a (now $d$-dimensional multivariate) Gaussian, and written in terms
of this Chi-squared quantity we can express it as:

$\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)=\frac{1}{\sqrt{(2\pi)^{d}\mathrm{det}(\bm{\Sigma})}}\exp\left[-\frac{1}{2}\chi^{2}\right],$
(4.13)

where
$\mathrm{det}(\bm{\Sigma})=\displaystyle\prod_{c}^{d}\sigma_{c}^{2}=\left(\displaystyle\prod_{c}^{d}\sigma_{c}\right)^{2}$.
Following the same logic trail as we did before, now we perform a
$k$-dimensional Taylor expansion of Chi-squared about
$\vec{\theta}_{\mathrm{true}}$ which looks like:

$\chi^{2}=\ \ \chi^{2}\bigg{\rvert}_{\vec{\theta}_{\mathrm{true}}}\ \ +\ \
\frac{\partial}{\partial\theta_{i}}\left[\chi^{2}\right]\bigg{\rvert}_{\vec{\theta}_{\mathrm{true}}}\Delta\theta^{i}\
\ +\ \
\frac{1}{2}\frac{\partial^{2}}{\partial\theta_{i}\theta_{j}}\left[\chi^{2}\right]\bigg{\rvert}_{\vec{\theta}_{\mathrm{true}}}\Delta\theta^{i}\Delta\theta^{j}\
\ +\ \ \mathcal{O}\left(\Delta\theta^{3}\right),$ (4.14)

and then take the expectation value of this expression, where the expectation
value now looks like:

$\displaystyle\langle\ \cdot\ \rangle$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\ \cdot\ \right]\ \
\mathscr{L}\left(\vec{X}\middle|\vec{\theta}_{\mathrm{true}}\right)d^{d}X_{\mathrm{obs}},$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\ \cdot\ \right]\ \
\frac{1}{\sqrt{(2\pi)^{d}\mathrm{det}(\bm{\Sigma})}}\exp\left[-\frac{1}{2}\chi\left(\vec{\tilde{X}}\left(\vec{\theta}_{\mathrm{true}}\right)\right)^{2}\right]d^{d}\\!\\!\left[\tilde{X}\left(\vec{\theta}_{\mathrm{true}}\right)\right].$
(4.15)

Remember, we are still evaluating the expectation value at the true parameter
values $\vec{\theta}_{\mathrm{true}}$. The resulting integrals are certainly
uglier, but keeping careful collection of like terms they all work out in the
same way as before. The zeroth order term looks like:

$\displaystyle\left\langle\chi^{2}\bigg{\rvert}_{\vec{\theta}_{\mathrm{true}}}\right\rangle$
$\displaystyle=\int\limits^{\infty}_{-\infty}\left[\frac{\left(\tilde{X}^{c}\left(\vec{\theta}_{\mathrm{true}}\right)\right)^{2}}{\sigma^{2}_{c}}\right]\frac{1}{\sqrt{(2\pi)^{d}\mathrm{det}(\bm{\Sigma})}}\exp\left[-\frac{1}{2}\frac{\left(\tilde{X}^{c}\left(\vec{\theta}_{\mathrm{true}}\right)\right)^{2}}{\sigma^{2}_{c}}\right]d^{d}\\!\\!\left[\tilde{X}\left(\vec{\theta}_{\mathrm{true}}\right)\right],$
$\displaystyle=\int\limits^{\infty}_{-\infty}\frac{\left(\tilde{X}^{c}\left(\vec{\theta}_{\mathrm{true}}\right)\right)^{2}}{\sigma^{2}_{c}}\frac{1}{\sqrt{2\pi}\sigma_{c}}\exp\left[-\frac{1}{2}\frac{\left(\tilde{X}^{c}\left(\vec{\theta}_{\mathrm{true}}\right)\right)^{2}}{\sigma^{2}_{c}}\right]d\\!\\!\left[\tilde{X}^{c}\left(\vec{\theta}_{\mathrm{true}}\right)\right]$
$\displaystyle\qquad\times\ \
\int\limits^{\infty}_{-\infty}\frac{1}{\sqrt{(2\pi)^{d-1}\displaystyle\prod_{a}^{d-1}\sigma_{a}^{2}}}\exp\left[-\frac{1}{2}\frac{\left(\tilde{X}^{a}\left(\vec{\theta}_{\mathrm{true}}\right)\right)^{2}}{\sigma^{2}_{a}}\right]d^{d-1}\\!\\!\left[\tilde{X}\left(\vec{\theta}_{\mathrm{true}}\right)\right],$
$\displaystyle=1\times 1=1.$ (4.16)

Here we have separated out like terms, the $\tilde{X}^{c}$ integral from the
rest of the $\tilde{X}^{a}$ for $a\neq c$ integrals (so note that the implied
sum doesn’t apply to the $c$ index in the third to last line above). As for
the first order term, we can see that when we take the first partial:

$\frac{\partial}{\partial\theta_{i}}\left[\chi^{2}\right]=2\frac{\tilde{X}^{c}}{\sigma^{2}_{c}}\frac{\partial
X^{c}}{\partial\theta_{i}}=\left(\frac{2}{\sigma^{2}_{c}}\frac{\partial
X^{c}}{\partial\theta_{i}}\right)\tilde{X}^{c},$

we will end up with an odd Gaussian in $\tilde{X}^{c}$ like we did in equation
4.8, which will kill off the first order term (as we would expect, if we are
sitting at an extremum like we found in the 1D case). So our result now looks
like:

$\displaystyle\left\langle\chi^{2}\right\rangle\quad$ $\displaystyle=\quad
1\quad+\quad
F_{ij}\Delta\theta^{i}\Delta\theta^{j}\quad+\quad\mathcal{O}\left(\Delta\theta^{3}\right),$
$\displaystyle\mathrm{where}\quad
F_{ij}=F_{ij}\left(\vec{\theta}_{\mathrm{true}}\right)\equiv\frac{1}{2}\left\langle\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\left[\chi^{2}\right]\bigg{\rvert}_{\vec{\theta}_{\mathrm{true}}}\right\rangle.$
(4.17)

Now rather than just a scalar quantity we have a full matrix, referred to as
the “Fisher (Information) Matrix,” whose elements are given by eqn 4.17. Also
like before, following the same steps and logic that lead up to equation 4.11
(remember, $\Delta\theta=\theta-\theta_{\mathrm{true}}$):

$\displaystyle-\left\langle\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}_{\mathrm{true}}\right)}\right\rangle$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\left\langle\chi^{2}\right\rangle,$
$\displaystyle=\frac{1}{2}\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\left[1+F_{ab}\Delta\theta^{a}\Delta\theta^{b}+\mathcal{O}\left(\Delta\theta^{3}\right)\right],$
$\displaystyle=\frac{1}{2}\frac{\partial}{\partial\theta_{i}}\left[F_{ab}\delta^{aj}\Delta\theta^{b}+F_{ab}\Delta\theta^{a}\delta^{bj}+\mathcal{O}\left(\Delta\theta^{2}\right)\right],$
$\displaystyle=\frac{1}{2}\frac{\partial}{\partial\theta_{i}}\left[F_{jb}\Delta\theta^{b}+F_{aj}\Delta\theta^{a}+\mathcal{O}\left(\Delta\theta^{2}\right)\right],$
$\displaystyle=\frac{1}{2}\left[F_{ji}+F_{ij}+\mathcal{O}\left(\Delta\theta^{2}\right)\right],$
$\displaystyle\approx F_{ij},$

where in the second-to-last line we used the fact that the Fisher matrix, as
we can see from equation 4.17 is a symmetric matrix, so $F_{ji}=F_{ij}$. This
means that we can define (approximately as we pointed out in the previous
section) that:

$F_{ij}=-\left\langle\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}_{\mathrm{true}}\right)}\right\rangle.$
(4.18)

Like with what we found in equation 4.10, in the high-statistical limit we
expect now that our posterior distribution should look like a $k$-dimensional
multivariate Gaussian distribution described using the Fisher matrix as:

$\displaystyle p\left(\vec{\theta}\middle|\vec{X}\right)$
$\displaystyle\sim\sqrt{\frac{\mathrm{det}(\bm{F})}{(2\pi)^{k}}}\exp\left[-\frac{1}{2}F_{ij}\Delta\theta^{i}\Delta\theta^{j}\right]=\sqrt{\frac{\mathrm{det}(\bm{F})}{(2\pi)^{k}}}\exp\left[-\frac{1}{2}\Delta\vec{\theta}^{T}\bm{F}\Delta\vec{\theta}\right],$
$\displaystyle=\sqrt{\frac{1}{(2\pi)^{k}\mathrm{det}(\bm{C})}}\exp\left[-\frac{1}{2}C^{-1}_{ij}\Delta\theta^{i}\Delta\theta^{j}\right]=\sqrt{\frac{1}{(2\pi)^{k}\mathrm{det}(\bm{C})}}\exp\left[-\frac{1}{2}\Delta\vec{\theta}^{T}\bm{C}^{-1}\Delta\vec{\theta}\right].$
(4.19)

Here $\bm{C}$ is the parameter covariance matrix, which as we can see is the
inverse of the Fisher matrix, $\bm{F}^{-1}=\bm{C}$. (In the 1D case, we found
the inverse of the Fisher information was the variance, $F^{-1}=\sigma^{2}$,
the natural generalization of that is that now we are talking in terms of
matrices). So inverting the Fisher matrix would give us all of the covariances
between each of the model parameters, and hence how much uncertainty we could
expect when trying to measure these parameters in an experiment. Once again,
this Fisher analysis takes us from a multivariate Gaussian likelihood in data-
space to a multivariate Gaussian posterior in parameter-space.

Finally, we note that all of the results presented in this generalization
faithfully reduce to the 1D case (in both data and parameters) that we gave in
Section 4.2.

###  4.4 Two Common Fisher Matrix Forms

In the literature the Fisher Matrix is often defined111Again, this is still a
sort of approximate definition, with the small $\Delta\theta$ assumption built
into it as pointed out in Section 4.2. in one of two forms:

$\displaystyle F_{ij}\left(\vec{\theta}_{\mathrm{true}}\right)$
$\displaystyle=-\left\langle\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}_{\mathrm{true}}\right)}\right\rangle,$
$\displaystyle=\Bigg{\langle}\left(\frac{\partial}{\partial\theta_{i}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}_{\mathrm{true}}\right)}\right)\left(\frac{\partial}{\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}_{\mathrm{true}}\right)}\right)\Bigg{\rangle}.$
(4.20)

Here we prove that the second form of the Fisher matrix as stated here is
equivalent to the first.

These two forms give two ways of describing what the Fisher Matrix is. In the
first form in equation 4.20, it is the “expected value of the observed
information,” where the observed information is the negative Hessian matrix of
the log-likelihood. In the second form in equation 4.20, it is the “variance
of the score,” where the score is the gradient of the log-likelihood with
respect to the parameter vector (and the variance is the expectation value of
a squared deviation of an R.V.).

Thinking of the Fisher matrix as the expected value of the observed
information is helpful, since when deriving the Fisher matrix we emphasized
how by Taylor expanding $\chi^{2}$ about the true parameters
$\theta_{\mathrm{true}}$ and taking the expectation value of that also around
$\theta_{\mathrm{true}}$, we were localizing ourselves to the behavior of the
function right around its maximum, that is its curvature right around its
maximum (which is effectively what our Fisher matrix is). The Hessian matrix
is the local curvature of a function of multiple variables, hence the natural
conceptual connection.

Proof First note:
$\displaystyle-\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}$
$\displaystyle=-\frac{\partial}{\partial\theta_{i}}\left[\frac{\partial}{\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\right],$
$\displaystyle=-\frac{\partial}{\partial\theta_{i}}\left[\frac{1}{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\frac{\partial}{\partial\theta_{j}}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)\right],$
$\displaystyle=\frac{1}{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)^{2}}\left(\frac{\partial}{\partial\theta_{i}}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)\right)\left(\frac{\partial}{\partial\theta_{j}}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)\right)$
$\displaystyle\hskip 71.13188pt-\ \
\frac{1}{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right),$
$\displaystyle=\left(\frac{\partial}{\partial\theta_{i}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\right)\left(\frac{\partial}{\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\right)$
$\displaystyle\hskip 71.13188pt-\ \
\frac{1}{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right).$
(4.21) Next we take the expectation value of both sides of this equation, and
consider what happens in that process to the second term on the RHS:
$\displaystyle\left\langle\frac{1}{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)\right\rangle$
$\displaystyle=\int\limits^{\infty}_{-\infty}\frac{1}{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\left(\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)\right)\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)d^{n}x,$

$\displaystyle\hskip
28.45274pt=\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\int\limits^{\infty}_{-\infty}\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)d^{n}x,$
$\displaystyle\hskip
28.45274pt=\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}(1)\equiv
0.$ (4.22) This leaves us with the following expression after taking the
expectation value of equation 4.21:
$-\left\langle\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\right\rangle=\left\langle\left(\frac{\partial}{\partial\theta_{i}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\right)\left(\frac{\partial}{\partial\theta_{j}}\ln{\mathscr{L}\left(\vec{X}\middle|\vec{\theta}\right)}\right)\right\rangle.$
Q.E.D.

###  4.5 Adding Prior Knowledge to the Fisher Matrix

As we stated at the start of Section 4.1, the Fisher analysis makes that
assumption that in the high-statistical limit of repeating our experiment many
times,the effect of the prior on the posterior should become less significant.
Because we are effectively “zooming” in on a local region of parameter space
when we perform our Taylor expansion approximation, one of the main ideas is
that the local curvature of the prior should remain effectively constant,
meaning only the curvature of the likelihood is what is important and
contributing to the posterior. These assumptions and this derivation then
leads to our Fisher matrix, whose inverse gives us the best-case-scenario for
how well we can measure the parameters in our model. Hence the Fisher matrix
itself is effectively a best-case-scenario measurement estimate of the
parameters with no prior knowledge built into its derivation.

In practice, however, we may be performing an experiment where we do know
additional information about the parameters (perhaps something related to
their physical nature), which can help further constrain our best-case-
scenario, and make it even more precise (Wittman). In this case we may want to
build back in some of that prior knowledge, and “update” our Fisher matrix
approximation. To do this we can choose to take our high-statistical
likelihood and multiply by an additional prior. If that prior knowledge is
Gaussian distributed, we can express our updated Fisher matrix posterior as:

$\displaystyle p\left(\vec{\theta}\middle|\vec{X}\right)$
$\displaystyle\sim\textsf{C}\exp\left[-\frac{1}{2}\Delta\vec{\theta}^{T}\bm{F}\Delta\vec{\theta}\right]\times\exp\left[-\frac{1}{2}\Delta\vec{\theta}^{T}\textbf{Cov}^{-1}\Delta\vec{\theta}\right],$
$\displaystyle=\textsf{C}\exp\left[-\frac{1}{2}\Delta\vec{\theta}^{T}\Big{(}\bm{F}+\textbf{Cov}^{-1}\Big{)}\Delta\vec{\theta}\right],$
(4.23)

where Cov is some new known prior covariance matrix, and C is the new
normalizing factor. So in an analysis, if we have the Fisher matrix for our
model and experiment, we can include additional prior constraints by adding to
the Fisher matrix the inverted Covariance matrix of our prior. Inverting this
new Fisher matrix,
$\bm{F}^{\bm{{}^{\prime}}}=\Big{(}\bm{F}+\textbf{Cov}^{-1}\Big{)}$, would give
us even more constrained covariances on our experiment’s parameters,
$\textbf{Cov}^{\bm{{}^{\prime}}}=\Big{(}\bm{F}+\textbf{Cov}^{-1}\Big{)}^{-1}$.

In reality, often we may find that from prior knowledge we can constrain some
of our parameters, but not all of them. This means we cannot write a
covariance matrix Cov with the same dimensions as our Fisher matrix $\bm{F}$,
which can be inverted $\textbf{Cov}^{-1}$, because some of the columns/rows
would contain only zeros. Hence we can’t simply write down equation 4.23.
Instead in practice we simply add the covariances directly to the appropriate
elements of the Fisher matrix.

###  4.6 The Timing Residual Fisher Matrix

Now let’s apply this Fisher matrix analysis to our timing residual model, in
order to predict how well pulsar timing experiments will be able to constrain
the model parameters of our different timing models (from Section 3.4). Once
again, since the Fisher matrix is the inverse of the parameter covariance
matrix, $\bf{F}=\bf{C}^{-1}$, inverting the Fisher matrix gives us the
covariances of all of our parameters, which tells us how well we should be
able to measure them experimentally.

Let our observed pulsar timing residual data (indexed by $a$, $b$, and $c$
below) be represented by a $d$-dimensional R.V.,
$\overrightarrow{\mathrm{Res}}=\left\\{\mathrm{Res}_{c}:c=1,2,\ldots
d\right\\}$. As we did previously, let the parameters of the timing residual
model in Table 3.2 (indexed by $i,j$) be described by a $k$-dimensional
$\vec{\theta}$. Therefore the true timing residual we are looking for is
burried in random Gaussian noise,
$\overrightarrow{\mathrm{Res}}=\overrightarrow{\mathrm{Res}}_{\mathrm{true}}+\vec{n}$,
which is what makes the observation an R.V. The timing data for every pulsar
at every observation can then be modeled as a function of $\vec{\theta}$ by
choosing a timing residual model from Section 3.4, organized into the vector
$\overrightarrow{\mathrm{Res}}\left(\vec{\theta}\right)$.222Note that in our
studies, we will simulate our own data sets rather than use real data.
Therefore we must “inject” a source by simulating our own “true” timing
residual. This is done by taking our desired injection parameters
$\theta_{\mathrm{true}}$ and inputting them into our desired model function,
that is
$\overrightarrow{\mathrm{Res}}_{\mathrm{true}}=\overrightarrow{\mathrm{Res}}\left(\vec{\theta}_{\mathrm{true}}\right)$.
The goal is then to see how well we can recover our injected parameters. Let’s
assume our measurements of the timing residual have uncorrelated Gaussian
timing uncertainties $\sigma$, with covariance matrix given by
$\bm{\Sigma}=\mathrm{diag}\left(\sigma^{2}_{1},\sigma^{2}_{2},\dots,\sigma^{2}_{d}\right)$,
as shown in equation 4.12.

Once again using implied Einstein summation notation as we did earlier, our
Chi-squared goodness-of-fit quantity can be written as:

$\chi^{2}\equiv\overrightarrow{\widetilde{\mathrm{Res}}}^{T}\bm{\Sigma}^{-1}\overrightarrow{\widetilde{\mathrm{Res}}}=\Sigma^{-1}_{ab}\widetilde{\mathrm{Res}}^{a}\widetilde{\mathrm{Res}}^{b}=\frac{1}{\sigma^{2}_{c}}\left(\widetilde{\mathrm{Res}}^{c}\right)^{2}=\frac{1}{\sigma^{2}_{c}}\bigg{(}\mathrm{Res}^{c}\left(\vec{\theta}\right)-\mathrm{Res}^{c}\bigg{)}^{2},$
(4.24)

where
$\overrightarrow{\widetilde{\mathrm{Res}}}\equiv\overrightarrow{\mathrm{Res}}\left(\vec{\theta}\right)-\overrightarrow{\mathrm{Res}}$.
One of the key ideas here is that because our noise $\vec{n}$ is Gaussian
distributed, and because we are modeling our data as the true signal plus this
noise (with no correlations between timing data), then our timing residual
data is also Gaussian distributed. Therefore our likelihood should be a
Gaussian written in terms of this Chi-squared quantity as:

$\displaystyle\mathscr{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{\theta}\right)$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{d}\mathrm{det}(\bm{\Sigma})}}\exp{\left[-\frac{1}{2}\chi^{2}\right]},$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{d}\mathrm{det}(\bm{\Sigma})}}\exp{\Bigg{[}-\frac{1}{2}\frac{1}{\sigma^{2}_{c}}\left(\mathrm{Res}^{c}\left(\vec{\theta}\right)-\mathrm{Res}^{c}\right)^{2}\Bigg{]}},$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{d}\mathrm{det}(\bm{\Sigma})}}\exp{\Bigg{[}-\frac{1}{2}\frac{1}{\sigma^{2}_{c}}\left(\widetilde{\mathrm{Res}}^{c}\right)^{2}\Bigg{]}}.$
(4.25)

Now that we have our likelihood function defined, we need to compute the
Fisher matrix. We will use the first form of equation 4.20 for our derivation.
First note the following:

$\displaystyle\ln{\mathscr{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{\theta}_{\mathrm{true}}\right)}$
$\displaystyle=-\frac{1}{2}\ln\left[(2\pi)^{d}\mathrm{det}(\bm{\Sigma})\right]-\frac{1}{2}\frac{1}{\sigma^{2}_{c}}\left(\widetilde{\mathrm{Res}}^{c}\right)^{2},$
$\displaystyle\frac{\partial}{\partial\theta_{j}}\ln{\mathscr{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{\theta}_{\mathrm{true}}\right)}$
$\displaystyle=0-\frac{1}{2}\frac{1}{\sigma^{2}_{c}}2\left(\widetilde{\mathrm{Res}}^{c}\right)\frac{\partial\widetilde{\mathrm{Res}}^{c}}{\partial\theta_{j}}=-\frac{1}{\sigma^{2}_{c}}\left(\widetilde{\mathrm{Res}}^{c}\right)\frac{\partial\mathrm{Res}^{c}}{\partial\theta_{j}},$
$\displaystyle-\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ln{\mathscr{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{\theta}_{\mathrm{true}}\right)}$
$\displaystyle=+\frac{1}{\sigma^{2}_{c}}\left[\frac{\partial\mathrm{Res}^{c}}{\partial\theta_{i}}\frac{\partial\mathrm{Res}^{c}}{\partial\theta_{j}}+\widetilde{\mathrm{Res}}^{c}\frac{\partial^{2}\mathrm{Res}^{c}}{\partial\theta_{i}\partial\theta_{j}}\right],$

where all of the residual functions in these expressions on the RHS are
evaluated at the true set of parameters $\vec{\theta}_{\mathrm{true}}$ (like
what is shown explicitly on the LHS). Now we want to take the high-statistical
limit and average over the R.V., which here is
$\widetilde{\mathrm{Res}}\left(\vec{\theta}_{\mathrm{true}}\right)$,
remembering that the expectation value is also evaluated at the true
parameter. So the resulting integral that we want to solve for is:

$\displaystyle F_{ij}\left(\vec{\theta}_{\mathrm{true}}\right)$
$\displaystyle=-\left\langle\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\ln{\mathscr{L}\left(\overrightarrow{\mathrm{Res}}\mid\vec{\theta}_{\mathrm{true}}\right)}\right\rangle,$