# Gravitational Wave Timing Residual Models for Pulsar Timing Experiments

Casey McGrath

ABSTRACT

by

The University of Wisconsin-Milwaukee, 2021

Under the Supervision of Professor Jolien Creighton, PhD

The ability to detect gravitational waves now gives scientists and astronomers
a new way in which they can study the universe. So far, the LIGO and Virgo
scientific collaborations have been successful in detecting binary black hole
and binary neutron star mergers (Abbott et al., 2016, 2017). These types of
sources produce gravitational waves with frequencies of the order hertz to
millihertz. But there do exist other theoretical sources which would produce
gravitational waves in different parts of the frequency spectrum. Of these are
the theoretical mergers of supermassive black hole binaries (SMBHBs), which
could occur upon the merging of two galaxies with supermassive black holes at
their cores. Sources like these would produce gravitational waves generally
around the nanohertz regime, and the current main effort for detecting and
measuring these waves comes from pulsar timing experiments. Detection of
gravitational waves in these experiments would come as small fluctuations in
the otherwise extremely regular period of pulsars over a long period of time
(months to decades).

There are numerous goals for this dissertation. The first is to re-present
much of the fundamental physics and mathematics concepts behind the
calculations in this paper. While there exist many reference sources in the
literature, we simply try to offer a fully self-contained explanation of the
fundamentals of this research which we hope the reader will find helpful. It
is often a challenge when jumping into a new field of study to quickly learn
and understand the fundamentals (like the derivations of various formulae and
the assumptions behind the models), so if this dissertation can help future
readers to connect the dots between the blanks not filled in by other
literature sources, then this goal will be accomplished. The pedantic approach
to this dissertation is also helpful since much of the initial work for this
dissertation was theoretical development of the mathematical models used in
pulsar timing.

The next goal broadly speaking has been to combine the efforts of two previous
studies by Deng & Finn (2011) and Corbin & Cornish (2010) to further develop
the mathematics behind the currently used pulsar timing models for detecting
gravitational waves with pulsar timing experiments. Previous timing residual
models have first been derived assuming that the pulsar timing array receives
plane-waves coming from distant sources (with the notable exception of Deng &
Finn). Then these models can either treat the SMBHB as a monochromatic
gravitational wave source, or model the frequency evolution of the
gravitational waves over the thousands to tens of thousands of years it takes
light to travel from the pulsar to the Earth. Our research began by first
generalizing these models by removing the plane-wave assumption. In Chapter 3
we classify four regimes of interest (Figure 3.5), governed by the main
assumptions made when deriving each regime. Of these four regimes the plane-
wave models are well established in previous literature. We add a new regime
which we label “Fresnel,” as we will show it becomes important for significant
Fresnel numbers describing the curvature of wavefronts.

With these mathematical models developed, in Chapter 6 we present the first
main study investigated which was to forecast the ability of future pulsar
timing experiments to probe and measure these Fresnel effects. Here we show
the constraints needed on the pulsar timing experiments themselves (largely
explained by the discussion in Chapter 5), and the types of precision
measurements which could theoretically be achieved.

Then we generalize our models to a cosmologically expanding universe in
Chapter 7. We show that in the fully general Fresnel frequency evolution
regime, the Hubble constant enters the model and can now be measured directly.
In this chapter we investigate what we will need of future experiments in
order to obtain a measurement of this parameter. This offers future pulsar
timing experiments the unique possibility of being able to procure a purely
gravitational wave-based measurement of the Hubble constant.

Finally, Chapter 8 shows the initial steps taken to extend this work in the
future, specifically for Doppler tracking experiments. The main goal of the
inclusion of this final section, which was not the primary focus of this
dissertation, is to point out how the mathematics and models derived in
Chapters 2 and 3 can be applied and extended more generally.

© Copyright by , 2021
All Rights Reserved

###### Table of Contents

  1. List of Abbreviations
  2. 1 Introduction
  3. 2 Gravitational Waves from a Binary Source
    1. 2.1 Multipole Expansion
    2. 2.2 Multipole Potentials
    3. 2.3 The Gravitational Wave Solution - General Form
    4. 2.4 The Gravitational Wave Solution - Binary System
    5. 2.5 Generalized Source Orientation and Location
    6. 2.6 Source Parameters and a Word about Notation
    7. 2.7 Two Models: Monochromatism and Frequency Evolution
      1. 2.7.1 Monochromatism
      2. 2.7.2 Frequency Evolution
  4. 3 The Continuous Wave Timing Residual
    1. 3.1 The Observed Pulsar Period
    2. 3.2 The Observable: The Pulsar Timing Residual
    3. 3.3 Time Retardation
    4. 3.4 Four Pulsar Timing Regimes
    5. 3.5 Plane-Wave Formalism (I)
      1. 3.5.1 The Plane-Wave Antenna Response
      2. 3.5.2 Plane-Wave, Monochromatic (IA)
      3. 3.5.3 Plane-Wave, Frequency Evolution (IB)
    6. 3.6 Fresnel Formalism (II)
      1. 3.6.1 Fresnel, Monochromatic (IIA)
      2. 3.6.2 Fresnel, Frequency Evolution - Conjecture (IIB)
    7. 3.7 Regime & Model Comparisons
    8. 3.8 Regime Search Parameters & Degeneracies
  5. 4 Fisher Matrix Analysis
    1. 4.1 High Statistical Limit Applied to Bayes’ Theorem
    2. 4.2 Conceptually Understanding the Fisher Information
    3. 4.3 Generalizing from the Fisher Information to the Fisher Matrix
    4. 4.4 Two Common Fisher Matrix Forms
    5. 4.5 Adding Prior Knowledge to the Fisher Matrix
    6. 4.6 The Timing Residual Fisher Matrix
    7. 4.7 Breaking the Timing Residual Fisher Matrix into Submatrices
    8. 4.8 Fisher Coordinate Transformation
  6. 5 Pulsar Distance Wrapping Problem
  7. 6 Investigation of the Fresnel Models
    1. 6.1 Simple Signal-to-Noise Ratio Analysis
    2. 6.2 Measuring Source Distance
      1. 6.2.1 Figure-of-Merit (F.O.M)
      2. 6.2.2 Distant Sources
    3. 6.3 Source Localization
    4. 6.4 Measuring Chirp Mass from a Monochromatic Source
    5. 6.5 Conclusions
  8. 7 Generalizing the Fresnel Models & Measuring the Hubble Constant
    1. 7.1 Incorporating Cosmological Expansion & the Hubble Constant
    2. 7.2 Generalizing From a Static to Expanding Universe
    3. 7.3 The Fresnel, Frequency Evolution Model - Revisited
    4. 7.4 Methods
    5. 7.5 Pulsar Distance Wrapping Problem (Revisited) & the Error Envelope
    6. 7.6 The 1 Source Problem
    7. 7.7 The 2 Source Problem
    8. 7.8 Conclusions
  9. 8 General Doppler Tracking
    1. 8.1 The Observed Clock Period and Timing Residual
    2. 8.2 Plane-Wave Formalism (I)
      1. 8.2.1 Plane-Wave, Monochromatic (IA)
      2. 8.2.2 Plane-Wave, Frequency Evolution (IB)
    3. 8.3 Fresnel Formalism (II)
      1. 8.3.1 Fresnel, Monochromatic (IIA)
    4. 8.4 Future Directions

###### List of Figures

  1. 2.1 Binary System: Orientation
  2. (a) Binary system orbiting in the $x$/$y$-plane with the origin placed at the system’s center-of-mass.
  3. (b) Generalized source orientation. Two Euler angles $\iota$ (“inclination”) and $\psi$ (“polarization”), and a phase offset $\theta_{0}$ in the binary system’s motion (see equations 2.34 and 2.38), fully generalize the binary’s orientation.
  4. 2.2 Binary System: Position
  5. (a) Fully generalized geometrical positioning of a gravitational wave source with respect to an observer. Angles $\theta$ and $\phi$ denote source sky position, and $\iota$ and $\psi$ indicate source inclination and polarization. Source appears “face-on” if $\iota=0$, and “edge-on” if $\iota=\pi/2$, and “face-off” if $\iota=\pi$.
  6. (b) The source basis vectors are chosen so that when $\theta=\phi=0$, then $\hat{r}=\hat{z}$, $\hat{\theta}=\hat{x}$, and $\hat{\phi}=\hat{y}$. In this case, if the source is not inclined (i.e. $\iota=0$), then the observer will see the binary orbiting in the counter-clockwise direction.
  7. 3.1 Toy Clock Experiment: Period Shift
  8. 3.2 Toy Clock Experiment: Timing Residual 1
  9. (a) Example of $\mathrm{Res}(t)$. In this case, the measured time intervals spend equal amounts of time being shorter than expected as they do longer than expected.
  10. 3.3 Toy Clock Experiment: Timing Residual 2
  11. 3.4 Time Retardation Approximations
  12. 3.5 Gravitational Wave Timing Regimes Classification
  13. 3.6 PTA Antenna Patterns Example
  14. 3.7 Time Scales Chronology
  15. 3.8 Monochromatic Timing Residual Models Comparison
  16. 3.9 Coalescence Time and Fresnel Number Contours
  17. 3.10 Gravitational Wave Timing Model Regime Comparisons 1
  18. 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\\}$.
  19. 3.12 Gravitational Wave Timing Model Regime Comparisons 2
  20. 4.1 Likelihood Function Visualization
  21. (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.
  22. (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.
  23. 5.1 Pulsar Distance Wrapping Problem Example
  24. 5.2 Pulsar Distance Uncertainties
  25. (a)
  26. (b)
  27. (c)
  28. (d)
  29. 6.1 Pulsar Timing Uncertainty
  30. (a)
  31. (b)
  32. 6.2 SNR Analysis: Source Orbital Frequency vs. Chirp Mass
  33. 6.3 PTA Figure-of-Merit 1
  34. 6.4 PTA Figure-of-Merit 2
  35. 6.5 $\mathrm{CV}_{R}$ vs. $R$, $\omega_{0}$, and $\mathcal{M}$
  36. 6.6 Fisher Triangle Plot Example - Source 1
  37. 6.7 Source Localization 1 (PTA Parameters)
  38. 6.8 Source Localization 2 (Source Parameters)
  39. 6.9 Source Localization 3 (Sky & Volume Localization)
  40. 6.10 Regime IIA - Measuring $\mathcal{M}$
  41. 7.1 Gravitational Wave Timing Regime Classification - Revisited
  42. 7.2 The Pulsar Distance “Error Envelope”
  43. 7.3 Classifying the Error Envelope
  44. 7.4 Uncertainties of $D_{L}$ vs. $D_{\mathrm{par}}$
  45. 7.5 1 Known Source - Pulsar Distance Recovery
  46. 7.6 2 Sources - Common Wrapping Cycle
  47. 7.7 2 Sources - $H_{0}$ Recovery (Fisher Survey)
  48. 7.8 2 Sources - $H_{0}$ Recovery (MCMC Example)

###### List of Tables

  1. List of Abbreviations
  2. 3.1 Gravitational Wave Timing Model Sub-Regimes
  3. 3.2 Model Regime Degeneracies

## List of Abbreviations

PTA |  Pulsar Timing Array
---|---
TOA |  Time-of-arrival
SMBHB |  Supermassive Black Hole Binary
NANOGrav |  North American Nanohertz Observatory for Gravitational Waves (a PTA)
LIGO |  Laser Interferometer Gravitational-wave Observatory
KAGRA |  Kamioka Gravitational Wave Detector
CV |  Coefficient of Variation
SNR |  Signal-to-Noise Ratio

## Acknowledgments

I am grateful to have known many wonderful people throughout my entire life,
who have stimulated my curiosity and nurtured me along the path which has now
culminated in the completion of this PhD work. This path for me began in 5th
grade, when I discovered my love for science and astronomy. The last nearly
two decades have been a remarkable journey of personal growth, and this PhD
marks the end of a large phase of my life. All of the people who I have known
and supported me during this time have my love and thanks for everything they
have done for me.

My time as a graduate student has been extremely trying, and not what I
envisioned. But even still there have been certain individuals who have given
me their support and help during this time, and for that I am very grateful. I
am thankful for the relationships I had with all of the graduate students in
my department - many of whom helped and supported me through their friendship.
Especially Alex McEwen, Siddharth Mohite, Deep Chatterjee, Ian Brown, Caitlin
Rose, and Kristina Islo who always had my back and would help me when I needed
it.

From the encouragement of some of the previously mentioned individuals, during
my PhD work I began seeing a therapist for help with my mental health. This
was the single greatest decision I made during this period of my life. The
friendship I found with Katt Cochran has been one of the greatest I have ever
known, and what I needed the most in my life at this time.

Of the post-docs in the department I am especially grateful to have known and
worked with Angela Van Sistine, Megan DeCesar, Shaon Ghosh, Shasvath Kapadia,
Danielle Berg, Annalisa Citro, Laleh Sadeghian, and Sydney Chamberlin. I am
especially grateful to have worked with them in organizing our public outreach
group “CoffeeShop Astrophysics,” to bring science and a love for astronomy to
the greater Milwaukee community.

In the department I am very thankful to Heidi Matera, Robert Wood, Tonia
Klein, and Swapnil Tripathi for their immense support.

I am very grateful to have met and collaborated with Dan D’Orazio near the end
of my PhD work. Working with him was what I had imagined my PhD was going to
be like when I had started my time here, so I am very glad that we met.
Additionally, Luke Kelley is an incredibly important person to me for the help
and support that he gave me. Being able to open up to him and have honest
conversations was very much needed and refreshing. I am also thankful for
Michael Lam, Stephen Taylor, and Joe Swiggum for their help throughout my
research work, answering my questions and pointing me to helpful resources
when I needed it.

Much of my work required extensive computing capabilities for my numerical
simulations, which was accomplished using the University of Wisconsin-
Milwaukee’s High Performance Computing resources on campus. I’d especially
like to thank Darin Peetz and Daniel Siercks who maintain and run these
resources, as they were always extremely willing to help answer any and all of
my questions, and moreover teach me more about high performance computing so
that I could better improve my skills.

Lastly, my PhD work was supported by the National Science Foundation (NSF)
through the NANOGrav collaboration Physics Frontier Center award NSF
PHY-1430284, as well as the awards NSF PHY-1607585 and NSF PHY-1912649. This
work also received support from the Wisconsin Space Grant Consortium Graduate
and Professional Research Fellowship Program under the NASA Training Grant
#NNX15AJ12H, as well as from the University of Wisconsin-Milwaukee, for which
I am grateful for access to their computational resources supported by NSF
PHY-1626190.

##  Chapter 1 Introduction

The mathematics and physics presented in this chapter are generally true for
any binary black hole system. However, as this work is building towards a
pulsar timing specific study, we will specifically focus our discussion
towards “supermassive black hole binaries” (SMBHBs).

It is theorized that during the collision of galaxies with supermassive black
holes at their centers, if the two black holes come close enough then they
could form a binary pair (a SMBHB), and their orbital interaction would cause
massive distortions in spacetime which would propagate away as gravitational
waves. SMBHB sources would likely be slowly orbiting each other with periods
ranging from months to decades, and hence they are thought to produce
gravitational waves in and around the nanohertz regime.

Current efforts to detect these types of gravitational waves comes from pulsar
timing techniques. Some pulsars in our Galaxy can be used as extremely regular
clocks when their pulses are timed by radio telescopes. By timing a collection
of such pulsars, known as a pulsar timing array (PTA), over long periods of
time, scientists can look for small deviations in their otherwise very precise
“clock ticks.” It is understood that the presence of gravitational waves
hitting the timed pulsars would induce a periodic redshifting and blueshifting
of the timing signals as measured on Earth (discussed in detail in Chapter 3).
Hence by looking for such deviations in the timing signals of pulsars,
scientists hope to find the evidence of the presence of gravitational waves. A
discovery such as this would be highly complimentary to ground-based detectors
such as LIGO (LIGO Scientific Collaboration et al., 2015), Virgo (Acernese et
al., 2015), and KAGRA (Akutsu et al., 2018), as it would be in a gravitational
frequency regime that these experiments cannot probe.

SMBHBs would produce what are called “continuous waves,” as the black holes
would likely be sufficiently far apart that as they orbited each other the
loss of energy in the system radiated away as gravitational waves would not
cause a significant collapse of their orbit. As such these sources would
maintain a nearly constant or “continuous” gravitational wave frequency, hence
the name of this type of radiation. It would not be until later times in their
binary evolution that the frequency of the orbit would begin to evolve notably
over the observation time scale (a process called “frequency chirping”). And
unlike ground-based experiments, PTAs would not be sensitive to the much later
inspiral or merger of such objects where the frequency chirping is very strong
as the binary coalesces. However, to try and head-off any misconceptions early
on, we are referring to the frequency of the gravitational waves measured on
observational time scales. As we will see later in Chapter 3, the
gravitational wave frequency can evolve in the measured timing residual of a
pulsar due to the nature of the incredibly large physical scale of a PTA
experiment.

Currently, pulsar timing collaborations such as the North American Nanohertz
Observatory for Gravitational Waves (NANOGrav), the Parkes Pulsar Timing
Array, and the European Pulsar Timing Array have been regularly collecting
timing data on numerous pulsars for over a decade, which makes it hopeful that
a confirmed detection of gravitational waves using this method will be
possible (Zhu et al., 2014; Babak et al., 2015; Aggarwal et al., 2019). While
the first detection may be the cumulative effect of gravitational waves coming
from many sources all across the sky as a stochastic background, in this work
we focus on the possibility of detecting individually loud continuous wave
signals coming from SMBHBs.

Publications Connected to this Dissertation
The work presented in Chapters 3, 4, 5, and 6 was published in the following
paper: McGrath, C., & Creighton, J. (2021). Fresnel Models for Gravitational
Wave Effects on Pulsar Timing. Monthly Notices of the Royal Astronomical
Society.
$\rightarrow$ https://doi.org/10.1093/mnras/stab1417
$\rightarrow$ https://arxiv.org/abs/2011.09561 At the time of finishing and
submitting this dissertation, the work in Chapter 7 was being prepared in
multiple manuscripts to submit for publication. $\longrightarrow$ Update: The
following paper presents the final results of this work: McGrath, C.,
D’Orazio, D., & Creighton, J. (2022). Measuring the Hubble Constant with
Double Gravitational Wave Sources in Pulsar Timing. Monthly Notices of the
Royal Astronomical Society.
$\rightarrow$ https://doi.org/10.1093/mnras/stac2593
$\rightarrow$ https://arxiv.org/abs/2208.06495

##  Chapter 2 Gravitational Waves from a Binary Source

###  2.1 Multipole Expansion

In order to build a mathematical base and help gain deeper physical insights
into the quadrupole formula (equation 2.17) and gravitational waves, we begin
by looking at the multipole expansion in general. The mathematics here are
extremely important and useful in physics, especially in fields like
electromagnetism. However, understanding the multipole expansion can often be
very difficult the first time you learn it, as the fundamental ideas can
easily be lost in the sometimes tedious and difficult math. Here, we aim to
cut directly through that math and show that fundamentally the multipole
expansion is a rather simple idea.

At its heart all the multipole expansion is, is a Taylor expansion of the
function:

$\displaystyle f\left(\vec{x}\right)$
$\displaystyle=\frac{1}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|},$ (2.1)
$\displaystyle=\frac{1}{\sqrt{\left(x-x^{{}^{\prime}}\right)^{2}+\left(y-y^{{}^{\prime}}\right)^{2}+\left(z-z^{{}^{\prime}}\right)^{2}}},$
$\displaystyle\equiv\frac{1}{\sqrt{\Delta x^{2}+\Delta y^{2}+\Delta z^{2}}},$
$\displaystyle=\frac{1}{\sqrt{\delta_{ab}\Delta x^{a}\Delta
x^{b}}}=\frac{1}{\sqrt{\Delta x_{a}\Delta x^{a}}},$

where $\vec{x}$ is the field point of interest, and $\vec{x}^{{}^{\prime}}$ is
the “source” location. In the final line, Einstein index notation (and implied
summations) are used. This is a function of three variables (in Cartesian
coordinates, $x$, $y$, and $z$), so the Taylor expansion will be multivariable
Taylor expansion. We are going to imagine that the source is far from the
origin, and Taylor expand the function about the origin, that is about
$\vec{x}=0$.

Looking at equation 2.1, notice the many different ways in which we can
express this function. We can write it in terms of the vectors $\vec{x}$ and
$\vec{x}^{{}^{\prime}}$, in terms of the variables $x$, $y$, and $z$, or by
using index notation over the variables $x^{i}$. Where this derivation
typically becomes difficult is choosing the mathematical tool you are going to
use to perform the Taylor series expansion. Arguably the most straight forward
approach is to write everything in terms of $x$, $y$, and $z$, and Taylor
expand in all of those variables near $0$. Conceptually this is the easiest,
but mathematically this ends up getting really ugly, really fast.
Alternatively we could use the vectorized version of the Taylor series formula
and apply vector derivatives, but personally I find this formula conceptually
harder to grasp. In my opinion, the easiest and best way to do this derivation
in a manner that keeps in sight of all of the important conceptual and
mathematical ideas is to exploit index notation. We often don’t learn index
summation notation until after first learning the multipole expansion, which I
believe is why this derivation is often not presented using this notation.
However, it helps to retain the most basic expression of the Taylor series
expansion, removes some of the ugliness of using vector notation, and requires
the least amount of written work. The trade off is simply first learning how
to use the notation properly, but after that this derivation is pretty
straight forward.

We need to Taylor expand equation 2.1, which can be written incredibly
concisely (with index notation!) as:

$f\left(\vec{x}\right)\quad\approx\quad
f\Big{\rvert}_{\vec{x}=0}\quad+\quad\partial_{i}f\Big{\rvert}_{\vec{x}=0}x^{i}\quad+\quad\frac{1}{2!}\partial_{i}\partial_{j}f\Big{\rvert}_{\vec{x}=0}x^{i}x^{j}\quad+\quad\ldots,$
(2.2)

where $\partial_{i}=\frac{\partial}{\partial x^{i}}$. The math for the needed
derivatives is given explicitly in the “Working out the terms” box below. When
we evaluate the terms at $\vec{x}=0$, note that we define the distance to the
source as
$R\equiv\sqrt{x^{{}^{\prime}}_{a}x^{{}^{\prime}a}}=\left|\vec{x}^{{}^{\prime}}\right|$.
Therefore, the multipole (Taylor) expansion of equation 2.1 is:

$f\left(\vec{x}\right)\quad\approx\quad\underbrace{\frac{1}{R}}_{\text{``monopole''}}\quad+\quad\underbrace{\frac{x^{{}^{\prime}}_{i}}{R^{3}}x^{i}}_{\text{``dipole''}}\quad+\quad\underbrace{\frac{3}{2R^{5}}\left[x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}-\frac{1}{3}\delta_{ij}R^{2}\right]x^{i}x^{j}}_{\text{``quadrupole''}}\quad+\quad\ldots,$
(2.3)

where we have given the names of each order term as they will be referred to
in what follows.

Working out the terms
Here we work out explicitly the needed derivatives in our Taylor series.
Really it is just careful application of the chain and product rules. Note
that $\partial_{i}\Delta
x^{a}=\frac{\partial\left(x^{a}-x^{{}^{\prime}a}\right)}{\partial
x^{i}}=\delta_{i}^{a}$, and $\partial_{i}\Delta
x_{a}=\frac{\partial\left(x_{a}-x_{a}^{{}^{\prime}}\right)}{\partial
x^{i}}=\delta_{ia}$. $\displaystyle\partial_{i}f$
$\displaystyle=\partial_{i}\frac{1}{\sqrt{\Delta x_{a}\Delta
x^{a}}}=-\frac{1}{2}\frac{\delta_{ia}\Delta x^{a}+\Delta
x_{a}\delta^{ia}}{\left[\Delta x_{a}\Delta
x^{a}\right]^{3/2}}=-\frac{1}{2}\frac{2\Delta x_{i}}{\left[\Delta x_{a}\Delta
x^{a}\right]^{3/2}}=-\frac{\Delta x_{i}}{\left[\Delta x_{a}\Delta
x^{a}\right]^{3/2}},$ $\displaystyle\partial_{i}\partial_{j}f$
$\displaystyle=-\partial_{i}\frac{\Delta x_{j}}{\left[\Delta x_{a}\Delta
x^{a}\right]^{3/2}}=\frac{-\delta_{ij}}{\left[\Delta x_{a}\Delta
x^{a}\right]^{3/2}}-\Delta x_{j}\cdot-\frac{3}{2}\frac{2\Delta
x_{i}}{\left[\Delta x_{a}\Delta x^{a}\right]^{5/2}},$
$\displaystyle=\frac{-\delta_{ij}}{\left[\Delta x_{a}\Delta
x^{a}\right]^{3/2}}+3\frac{\Delta x_{i}\Delta x_{j}}{\left[\Delta x_{a}\Delta
x^{a}\right]^{5/2}},$ $\displaystyle=\frac{3}{\left[\Delta x_{a}\Delta
x^{a}\right]^{5/2}}\left[\Delta x_{i}\Delta x_{j}-\frac{1}{3}\delta_{ij}\Delta
x_{a}\Delta x^{a}\right].$ Now evaluating these derivatives at $\vec{x}=0$,
and using
$R\equiv\sqrt{x^{{}^{\prime}}_{a}x^{{}^{\prime}a}}=\left|\vec{x}^{{}^{\prime}}\right|$,
we find: $\displaystyle\partial_{i}f\Big{\rvert}_{\vec{x}=0}$
$\displaystyle=+\frac{x^{{}^{\prime}}_{i}}{R^{3}},$
$\displaystyle\partial_{i}\partial_{j}f\Big{\rvert}_{\vec{x}=0}$
$\displaystyle=\frac{3}{R^{5}}\left[x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}-\frac{1}{3}\delta_{ij}R^{2}\right].$

###  2.2 Multipole Potentials

The multipole expansion is often used to write out the electric or
gravitational potentials, due to a distribution of charge or mass,
respectively (Creighton & Anderson, 2011; Harrison, 2020). For continuous
distributions, these two expressions are computed in the same way, simply with
different scaling factors out front:

$V\left(\vec{x}\right)=\begin{cases}-G{\int}\frac{dm\left(\vec{x}^{{}^{\prime}}\right)}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|}&=-G{\int}\frac{\rho\left(\vec{x}^{{}^{\prime}}\right)}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|}d^{3}x^{{}^{\prime}},\qquad\text{(Gravitational
Potential)}\\\\[14.0pt]
\frac{1}{4\pi\epsilon_{0}}{\int}\frac{dm\left(\vec{x}^{{}^{\prime}}\right)}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|}&=\frac{1}{4\pi\epsilon_{0}}{\int}\frac{\rho\left(\vec{x}^{{}^{\prime}}\right)}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|}d^{3}x^{{}^{\prime}}.\qquad\text{(Electric
Potential)}\end{cases}$ (2.4)

Focusing our attention on the gravitational potential (an analogous approach
can be taken for the electric potential), we can multipole expand the
denominator using our result in equation 2.3 and write:

$\displaystyle V\left(\vec{x}\right)$
$\displaystyle\quad=\quad-\frac{GM}{R}\quad-\quad\frac{Gx^{i}D_{i}}{R^{3}}\quad-\quad\frac{3Gx^{i}x^{j}\textit{\sout{I}}_{ij}}{2R^{5}}\quad-\quad\ldots,$
(2.5) $\displaystyle\mathrm{where}\quad\begin{cases}\begin{tabular}[]{l l
c}$M$&$\equiv{\int}\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}$
,&(Monopole)\\\\[6.0pt]
$D_{i}$&$\equiv{\int}x^{{}^{\prime}}_{i}\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}$
,&(Dipole)\\\\[6.0pt]
$\textit{\sout{I}}_{ij}$&$\equiv{\int}\left[x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}-\frac{1}{3}\delta_{ij}R^{2}\right]\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}$
,&(Quadrupole)\end{tabular}\end{cases}$ (2.6)

The monopole term in our expression is the standard Newtonian result for a
point mass, but with more complicated mass distributions we gain higher order
corrections.

There is a crucially important difference, however, for a mass distribution in
gravity as compared to a charge distribution in electromagnetism. We can
always choose a center-of-mass coordinate system for mass distributions. The
consequence of this is that the dipole term vanishes identically when working
in this chosen coordinate system. Conceptually this is because in a center-of-
mass coordinate system, by definition the mass is equally distributed in all
directions. Therefore integrating up all of the mass along each direction in
$D_{i}={\int}x^{{}^{\prime}}_{i}\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}={\int}x^{{}^{\prime}}_{i}dm\left(\vec{x}^{{}^{\prime}}\right)$
will just cancel out and net $0$.

We cannot always guarantee the same though for electromagnetism - that is, we
cannot always choose a center-of-charge coordinate system. As a simple
example, just consider two equal and opposite charges lying at $\pm x$. There
is no center-of-charge for this distribution, and thus there will always be
some non-zero dipole moment, no matter the chosen set of coordinates. This
helps point out one of the fundamental differences between gravity and
electromagnetism. In gravity there is no “negative mass particle” like in
electromagnetism with both positive and negative charges.

Another important point to take notice of here is the physical significance of
the quadrupole tensor. The quadrupole tensor in equation 2.6 is the negative
traceless version of the familiarly defined “moment of inertia tensor.” The
normal moment of inertia tensor describes all of the moments of inertia about
an object’s chosen orthogonal basis axes, and is defined as:

$\text{I}_{ij}\equiv{\int}\left[\delta_{ij}R^{2}-x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}\right]\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}.$
(2.7)

The trace of the moment of inertia tensor is:

$\displaystyle\text{Trace}\left(\text{I}_{ij}\right)=\delta^{ij}\text{I}_{ij}$
$\displaystyle={\int}\left[\delta^{ij}\delta_{ij}R^{2}-\delta^{ij}x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}\right]\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}={\int}\left[\delta^{i}_{i}R^{2}-x^{{}^{\prime}}_{i}x^{{}^{\prime}i}\right]\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}},$
$\displaystyle={\int}\Big{[}(3)R^{2}-R^{2}\Big{]}\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}\equiv{\int}2R^{2}\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}.$

If we remove the trace of this tensor from itself by subtracting $1/3$ of this
value off of the diagonal elements, we end up with the negative of our
quadrupole tensor. Mathematically, this statement is written as:

$\displaystyle\text{I}_{ij}-\frac{1}{3}\delta_{ij}\text{Trace}\left(\text{I}_{ij}\right)$
$\displaystyle={\int}\left[\delta_{ij}R^{2}-x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}-\frac{1}{3}\delta_{ij}\Big{(}2R^{2}\Big{)}\right]\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}},$
$\displaystyle={\int}\left[\frac{1}{3}\delta_{ij}R^{2}-x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}\right]\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}},$
$\displaystyle\equiv-\ \textit{\sout{I}}_{ij}.$

The tensor $\textit{\sout{I}}_{ij}$ is typically referred to as the “reduced
quadrupole tensor.” A similarly defined “quadrupole tensor” will appear later
and the two are defined as follows:

$\begin{cases}\begin{tabular}[]{l l
c}$\textit{I}_{ij}$&$\equiv{\int}x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}$
,&(Quadrupole Tensor)\\\\[6.0pt]
$\textit{\sout{I}}_{ij}$&$\equiv{\int}\left[x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}-\frac{1}{3}\delta_{ij}R^{2}\right]\rho\left(\vec{x}^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}$
.&(Reduced Quadrupole Tensor)\end{tabular}\end{cases}$ (2.8)

The regular quadrupole tensor does not have the traceless property like the
reduced quadrupole tensor.

###  2.3 The Gravitational Wave Solution - General Form

For the gravitational waves produced by the sources we are interested in (i.e.
massive binaries) we solve the Einstein field equations under the following
assumptions:

Assumptions: Field Equation Solutions 1. Cosmologically static universe. 2.
Cosmologically flat empty universe (hence the background metric is Minkowski).
3. Weak-field limit (gravitational waves are a metric perturbation on top of
Minkowski). 4. Small source compared to the distance to the observer and the
wavelength of the wave. 5. Slow-moving source (non-relativistic, no post-
Newtonian analysis required). 6. Transverse-traceless (“TT”) gauge. 7. Far
field approximation of the metric perturbation amplitude \- the binary is
sufficiently far from the field point of interest that
$|\vec{x}_{\mathrm{field}}-\vec{x}_{\mathrm{source}}|\approx R$ in the
amplitude of the metric perturbation. This assumption won’t be made when
evaluating the retarded time (which won’t affect the amplitude of the wave,
but rather the phase of the wave).

Assumptions 3-5 together make the “weak-small-slow” assumption for our source
- or in other words, we are in the “zeroth order Newtonian” regime. We will
not be considering sources near coalescence when this assumption might
otherwise be broken. There are numerous helpful resources including Maggiore
(2008); Creighton & Anderson (2011); Moore (2013); Carroll (2013) and Zee
(2013) that this discussion is based upon.

Under assumption 3, we begin by proposing a metric solution that looks like
flat spacetime with a small perturbation $h_{\mu\nu}$ on top of it:

$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},\quad\mathrm{where}\
\left|h_{\mu\nu}\right|\ll 1.$ (2.9)

The next step in the solution process is typically to re-express the metric
perturbation in a different form, known as the trace-reversed metric
perturbation:

$H_{\mu\nu}\equiv
h_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}h\qquad\longleftrightarrow\qquad
h_{\mu\nu}=H_{\mu\nu}-\frac{1}{2}\eta_{\mu\nu}H,$ (2.10)

where $h=\text{Trace}(h_{\mu\nu})=\eta^{\mu\nu}h_{\mu\nu}$. The reason for
this choice is simply mathematical convenience, as it helps further reduce the
complexity of the Einstein equations into a form more easily solved. If we
were to find a solution for $H_{\mu\nu}$, then equation 2.10 would give us the
conversion back to $h_{\mu\nu}$. Under the Lorenz gauge choice (which is part
of assumption 6), the resulting Einstein equations can be expressed in this
weak-field regime as:

$\begin{cases}\Box^{2}H_{\mu\nu}&=\frac{-16\pi G}{c^{4}}T_{\mu\nu},\\\\[6.0pt]
\partial^{\mu}H_{\mu\nu}&=0,\end{cases}$ (2.11)

where $T_{\mu\nu}$ is the “effective stress energy tensor” describing the
source of our gravitational wave solution. Here the “effective” differentiates
the tensor from the typically defined “stress energy tensor” in that it
contains $\mathcal{O}\left(h^{2}\right)$ terms which results above when
substituting equation 2.9 into the Einstein equations.

In general, we know that the sourced (inhomogeneous) wave equation solution
can be constructed by integrating the wave equation Green’s function.
Specifically, with the wave equation written in the following form for some
sourcing function $f\left(\vec{x},t\right)$, the solution can be written as:

$\displaystyle\Box^{2}\psi=-4\pi
f\left(\vec{x},t\right)\qquad\xrightarrow{\text{solution}}\qquad\psi\left(\vec{x},t\right)$
$\displaystyle=\int
G\left(\vec{x},\vec{x}^{{}^{\prime}},t,t^{{}^{\prime}}\right)f\left(\vec{x}^{{}^{\prime}},t^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}dt^{{}^{\prime}},$
$\displaystyle=\int\frac{\delta\left(t-t_{\mathrm{ret}}\right)}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|}f\left(\vec{x}^{{}^{\prime}},t^{{}^{\prime}}\right)d^{3}x^{{}^{\prime}}dt^{{}^{\prime}},$
$\displaystyle=\int\frac{f\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|}d^{3}x^{{}^{\prime}},$
(2.12)

where $t_{\mathrm{ret}}$ is the retarded time (see equation 3.21). Comparing
this to equation 2.11, we can see that the solution in terms of the trace-
reversed metric perturbation is:

$H_{\mu\nu}=\frac{4G}{c^{4}}\int\frac{T_{\mu\nu}\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)}{\left|\vec{x}-\vec{x}^{{}^{\prime}}\right|}d^{3}x^{{}^{\prime}}.$
(2.13)

At this point we have a few options. We could multipole expand this solution
using the results and ideas developed in Sections 2.1 and 2.2. But as we will
see, for our level of precision we actually only need the first term in the
multipole expansion, the “monopole” or “far-field” term (this is our
assumption 7). However, the reason we developed the previous two sections so
extensively is to show that the final answer we will find will actually be
expressed in terms of the quadrupole tensor!

So with this in mind, we can approximate the solution in equation 2.13 as:

$H_{\mu\nu}\approx\frac{4G}{c^{4}R}\int
T_{\mu\nu}\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}}.$
(2.14)

Now we are in a position to being evaluating all of the individual terms of
the metric tensor solution here. Since this is a wave solution to the wave
equation, we are looking for components that are behaving like a “wave.” The
$tt$\- and $it$/$ti$-components of the effective stress energy tensor
represent the energy density and the $i$-momentum density solutions. The
integrals of these components over all space which encompasses the source are:

$\displaystyle\int
T_{tt}\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}}$
$\displaystyle=\int\rho\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}}\equiv
M,$ $\displaystyle\int
T_{it}\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}}$
$\displaystyle=\int
T_{ti}\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}}\equiv
P_{i}\equiv 0,$

The $i$-momentum $P_{i}$ vanishes identically here because if we move the
coordinate system to the center-of-mass frame, then the source itself is not
moving - it has zero momentum. And integrating up the energy content of the
source just gives us its mass $M$, so neither of these terms will behave like
waves. So these two non-waving sets of solutions look like:

$\displaystyle H_{tt}$ $\displaystyle=\frac{4GM}{c^{4}R},$ (2.15)
$\displaystyle H_{it}=H_{ti}$ $\displaystyle=0.$ (2.16)

At this point we just need to focus our attention on the $ij$-components of
the metric perturbation, which will produce the final wave solution.

Thanks to the conservation of energy, the divergence theorem, and some clever
algebra/calculus, the $ij$-components of the metric perturbation solution can
actually be re-expressed in terms of the energy density ($tt$-component)
alone. Details of this are given explicitly in the “A Helpful Identity” box
below. The result is:

$\displaystyle H_{ij}\approx\frac{4G}{c^{4}R}\int
T_{ij}\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}}$
$\displaystyle=\frac{4G}{c^{4}R}\cdot\frac{1}{2}\frac{\partial^{2}}{\partial
t^{2}}\int
x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}T_{tt}\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}},$
$\displaystyle=\frac{2G}{c^{4}R}\frac{\partial^{2}}{\partial t^{2}}\int
x^{{}^{\prime}}_{i}x^{{}^{\prime}}_{j}\rho\left(\vec{x}^{{}^{\prime}},t_{\mathrm{ret}}\right)d^{3}x^{{}^{\prime}}.$

Now we can make our connection back to the multipole potentials of Section
2.2, and note that the final answer can be expressed in terms of the
quadrupole tensor from equation 2.8 as:

$H_{ij}\approx\frac{2G}{c^{4}R}\frac{\partial^{2}}{\partial
t^{2}}I_{ij}\left(t_{\mathrm{ret}}\right).$ (2.17)

This is known as the “quadrupole formula” - the gravitational wave solution to
the weak-field Einstein equations. So the gravitational waves produced by an
object are proportional to the acceleration of its quadrupole tensor, which as
we discussed can sort of be thought of as a different version of the moment of
inertia tensor of that object.

Fundamentally, this is different from electromagnetism. An insightful
discussion of the significance of equation 2.17 is given in the following
quote:

> “In contrast, the leading contribution to electromagnetic radiation comes
> from the changing dipole moment of the charge density. The difference can be
> traced back to the universal nature of gravitation. A changing dipole moment
> corresponds to motion of the center of density - charge density in the case
> of electromagnetism, energy density in the case of gravitation. While there
> is nothing to stop the center of charge of an object from oscillating,
> oscillation of the center of mass of an isolated system violates
> conservation of momentum. (You can shake a body up and down, but you and the
> earth shake ever so slightly in the opposite direction to compensate.) The
> quadrupole moment, which measures the shape of the system, is generally
> smaller than the dipole moment, and for this reason, as well as the weak
> coupling of matter to gravity, gravitational radiation is typically much
> weaker than electromagnetic radiation.” — Sean Carroll (Carroll, 2013)
>

With regards to Carroll’s conservation of momentum statement, we already
explained that the dipole moment in center-of-mass coordinates is necessarily
zero, but to further illustrate this consider a simple example of a system of
point particles with mass $m$. The dipole moment in equation 2.6 would be
$D_{i}=\sum mx^{{}^{\prime}}_{i}$. Imagine if the metric perturbation now were
proportional to the acceleration of the dipole moment instead of the
quadrupole moment. Then $\frac{d^{2}}{dt^{2}}D_{i}=\frac{d^{2}}{dt^{2}}\sum
mx^{{}^{\prime}}_{i}=\frac{d}{dt}\sum mv^{{}^{\prime}}_{i}=\frac{d}{dt}0=0$,
because of the conservation of momentum of the entire system. So in general,
gravitational radiation is such a small perturbation to spacetime largely
because gravity is such a weak force, but also since there is no dipole and
the quadrupole correction tends to be a smaller effect.

A Helpful Identity
First consider the following quantity and carefully apply the product rule to
expand it out term by term in the following way:
$\displaystyle\partial_{i}\partial_{j}\left(T^{ij}x^{m}x^{n}\right)$
$\displaystyle=\partial_{i}\Big{[}x^{m}x^{n}\partial_{j}T^{ij}+x^{n}T^{ij}\delta^{m}_{j}+x^{m}T^{ij}\delta^{n}_{j}\Big{]},$
$\displaystyle=\partial_{i}\Big{[}x^{m}x^{n}\partial_{j}T^{ij}+\left(x^{n}T^{im}+x^{m}T^{in}\right)\Big{]},$
$\displaystyle=x^{m}x^{n}\partial_{i}\partial_{j}T^{ij}+x^{n}\partial_{j}T^{mj}+x^{m}\partial_{j}T^{nj}+\partial_{i}\left(x^{n}T^{im}+x^{m}T^{in}\right),$
$\displaystyle=x^{m}x^{n}\partial_{i}\partial_{j}T^{ij}+\partial_{j}\left(x^{n}T^{mj}+x^{m}T^{nj}\right)-2T^{mn}+\partial_{i}\left(x^{n}T^{im}+x^{m}T^{in}\right),$
$\displaystyle=x^{m}x^{n}\partial_{i}\partial_{j}T^{ij}+2\partial_{i}\left(x^{n}T^{im}+x^{m}T^{in}\right)-2T^{mn}.$
Going from the 2nd to 3rd line we used the following twice:
$\partial_{i}\left(x^{n}T^{mj}\right)=\delta^{n}_{j}T^{mj}+x^{n}\partial_{j}T^{mj}=T^{mn}+x^{n}\partial_{j}T^{mj}\
\implies\
x^{n}\partial_{j}T^{mj}=\partial_{i}\left(x^{n}T^{mj}\right)-T^{mn}.$ Then
going from the 3rd to the 4th line we recognize that in the second term the
$j$-index is summed over, so we can re-label it to anything we want, namely we
can relabel it from $j\rightarrow i$. We can re-arrange the terms in this
expression and write:

$T^{mn}=\frac{1}{2}\Big{[}-\partial_{i}\partial_{j}\left(T^{ij}x^{m}x^{n}\right)+2\partial_{i}\left(x^{n}T^{im}+x^{m}T^{in}\right)+x^{m}x^{n}\partial_{i}\partial_{j}T^{ij}\Big{]}.$
Next let’s integrate this expression for $T^{mn}$ over all space. The first
crucial mathematical exploit we are going to make comes from the divergence
theorem. The generalized form of the divergence theorem states:
$\int\partial_{i}F^{ij}dV=\oint F^{ij}dA_{i}.$ The key idea here is that when
we apply the divergence theorem, we will be completely enclosing our source as
we integrate over all space, and far from our source, the effective stress
energy tensor will be zero, which is where the surface integral is evaluated.
So looking at the first two terms in this integral expression we have terms
that look like:
$\displaystyle\int\partial_{i}\partial_{j}\left(T^{ij}x^{m}x^{n}\right)dV$
$\displaystyle=\oint\oint\left(T^{ij}x^{m}x^{n}\right)dA_{i}dA_{j}\equiv 0$
$\displaystyle\int\partial_{i}\left(x^{n}T^{im}\right)dV$
$\displaystyle=\oint\left(x^{n}T^{im}\right)dA_{i}\equiv 0$ (It is important
here to understand that quantities like $T^{ij}x^{m}x^{n}$ can be thought of
as a single tensor objects - for example, we could just re-define this
quantity as some tensor $T^{ij}x^{m}x^{n}\equiv A^{ijmn}$). The final critical
step comes from the conservation of energy, which tells us
$\partial_{\mu}T^{\mu\nu}=0$. Expanding this statement out slightly means we
can write:
$\partial_{t}T^{tj}+\partial_{i}T^{ij}=0\quad\implies\quad\partial_{t}T^{tj}=-\partial_{i}T^{ij}.$

Two applications of this conservation statement (along with the symmetry
$T^{\mu\nu}=T^{\nu\mu}$), means we can replace the third term in the
expression for $T^{mn}$ with the following:
$x^{m}x^{n}\partial_{i}\partial_{j}T^{ij}=+x^{m}x^{n}\partial_{t}\partial_{t}T^{tt}=\partial_{t}\partial_{t}\left(x^{m}x^{n}T^{tt}\right).$
With this we are replacing spatial derivatives over the $ij$-components of our
effective stress energy tensor with two time derivatives over the energy
density component alone. This gives us our final result, that: $\int
T^{mn}dV=\frac{1}{2}\int\partial_{t}\partial_{t}\left(x^{m}x^{n}T^{tt}\right)dV=\frac{1}{2}\partial_{t}\partial_{t}\int\left(x^{m}x^{n}T^{tt}\right)dV.$

At this point we could specify an objects energy density, and solve equation
2.17 (we could then convert that back to the original metric perturbation
using equation 2.10). However we still need to get the solution into the
transverse-traceless gauge (“TT-gauge”) as stated in assumption 6. We know
that the reduced quadrupole tensor in equation 2.8 is the traceless version of
the quadrupole tensor, so that would get us half-way there. But while the
reduced quadrupole tensor is traceless, it may not necessarily be transverse
to the direction of interest along which the gravitational wave is
propagating.

We can cast any matrix into the TT-gauge by use of the special TT-projection
operator:

$\displaystyle\textbf{M}^{TT}$
$\displaystyle=\Lambda(\textbf{M})=\textbf{PMP}-\frac{1}{2}\mathrm{Trace}(\textbf{PM})\textbf{P},$
$\displaystyle\longrightarrow
M^{TT}_{ij}=\Lambda_{ijkl}M^{kl}=\left(P_{ik}P_{jl}-\frac{1}{2}P_{ij}P_{kl}\right)M^{kl},$
(2.18) $\displaystyle\text{where}\quad
P_{ij}=\delta_{ij}-\hat{n}_{i}\hat{n}_{j}.$

Here $\hat{n}$ is the desired direction of propagation of the gravitational
wave that we want to project onto. The TT-gauge gives us a number of useful
conceptual and mathematical results, namely:

$\begin{cases}\textit{I}^{TT}_{ij}&=\textit{\sout{I}}^{TT}_{ij},\\\\[6.0pt]
H^{TT}_{ij}&=h^{TT}_{ij}.\end{cases}$ (2.19)

With regards to the quadrupole tensors, the regular quadrupole tensor is
equivalent to the reduced quadrupole tensor after it is made traceless, so if
both of these are made traceless then transverse, then necessarily they are
equivalent. In practice we will often work with the reduced quadrupole tensor,
so this realization is helpful. As for the regular and trace-reversed metric
perturbations, the proof is given in the “TT Metric Projection” box.

With this we finalize our wave solution to the original proposed solution
(equation 2.9 to the Einstein equations. In the transvserse-traceless gauge,
the wave-like weak-field solutions are:

$h^{TT}_{ij}=\frac{2G}{c^{4}R}\frac{\partial^{2}}{\partial
t^{2}}\textit{\sout{I}}^{TT}_{ij}\left(t_{\mathrm{ret}}\right).$ (2.20)

TT Metric Projection
Conveniently, when both the regular and the trace-reversed metric
perturbations are cast into the TT-gauge they are equivalent, meaning once we
have the solution to one we have the solution to both. To see this we simply
apply the TT operator to equation 2.10 and see that:

$\displaystyle H^{TT}_{ij}=\Lambda_{ijkl}H^{kl}$
$\displaystyle=\Lambda_{ijkl}h^{kl}-\frac{1}{2}\Lambda_{ijkl}\eta^{kl}h,$
$\displaystyle=h^{TT}_{ij}-\frac{1}{2}{\Lambda_{ijk}}^{k}h,$
$\displaystyle=h^{TT}_{ij}-0.$ The trace of our projection operator along the
$kl$-indices is identically zero, which we can get by simply working out all
of the terms explicitly. First noting that
$P_{k}^{k}=\delta_{k}^{k}-n_{k}n^{k}=3-1=2$, we can write:
$\displaystyle{\Lambda_{ijk}}^{k}$
$\displaystyle=P_{ik}P_{j}^{k}-\frac{1}{2}P_{ij}P_{k}^{k}=P_{ik}P_{j}^{k}-P_{ij},$
$\displaystyle=\left(\delta_{ik}-n_{i}n_{k}\right)\left(\delta_{j}^{k}-n_{j}n^{k}\right)-P_{ij},$
$\displaystyle=\delta_{ij}-n_{j}n_{i}-n_{i}n_{j}+n_{i}n_{j}(1)-P_{ij},$
$\displaystyle=P_{ij}-P_{ij}=0.$

###  2.4 The Gravitational Wave Solution - Binary System

We begin with a simplified problem which we will then generalize. First
consider an arbitrary binary system aligned in the configuration shown in
Figure 2.1(a). The binary orbits in the $x/y$-plane, “face-on” with no
inclination to the $\hat{z}$-axis. Gravitational waves produced by the binary
propagate out in all directions, but we begin by focusing our attention to the
gravitational waves which propagate along the $\hat{z}$-axis. Therefore let’s
start by writing the solution in the TT-gauge for the gravitational waves
along that direction. To remind ourselves of the axis of propagation for our
TT-gauge, we will include the unit vector next to the “TT” symbol in the
notation that follows. (We start with $TT\hat{z}$, then later generalize to an
arbitrary $\hat{r}$ vector, denoted by $TT\hat{r}$).

(a) Binary system orbiting in the $x$/$y$-plane with the origin placed at the
system’s center-of-mass.

(b) Generalized source orientation. Two Euler angles $\iota$ (“inclination”)
and $\psi$ (“polarization”), and a phase offset $\theta_{0}$ in the binary
system’s motion (see equations 2.34 and 2.38), fully generalize the binary’s
orientation.

Figure 2.1:

The following assumptions are made about the binary system:

Assumptions: Binary Mechanics 1. The binary is circular. 2. The source is
located at a fixed angular sky position and distance from the Earth.

Under these assumptions, the metric perturbation solution to equation 2.20 can
be written compactly in the following form:

$\displaystyle h^{TT\hat{z}}_{ij}$
$\displaystyle=e^{\hat{z}+}_{ij}h_{+}+e^{\hat{z}\times}_{ij}h_{\times}=e^{\hat{z}\textsc{A}}_{ij}h_{\textsc{A}}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(2.21)
$\displaystyle\mathrm{where}\quad\begin{cases}h_{+}(t)\equiv-h(t)\cos\big{(}2\Theta(t)\big{)},\\\
h_{\times}(t)\equiv-h(t)\sin\big{(}2\Theta(t)\big{)},\end{cases}$ (2.22)
$\displaystyle\hskip
51.21504pth(t)\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega(t)^{2/3},$
(2.23)
$\displaystyle\qquad\qquad\begin{cases}e^{\hat{z}+}_{ij}\equiv\left(\begin{array}[]{ccc}1&0&0\\\
0&-1&0\\\
0&0&0\end{array}\right)=\hat{x}_{i}\hat{x}_{j}-\hat{y}_{i}\hat{y}_{j},\\\
\\\\[-25.0pt] e^{\hat{z}\times}_{ij}\equiv\left(\begin{array}[]{ccc}0&1&0\\\
1&0&0\\\
0&0&0\end{array}\right)=\hat{y}_{i}\hat{x}_{j}+\hat{x}_{i}\hat{y}_{j},\end{cases}$
(2.24)

where the metric decomposes into two polarizations, labeled “plus” and
“cross,” $\omega(t)$ and $\Theta(t)$ are the orbital frequency and phase
(specified in Section 2.7), $R$ is the Earth-source coordinate distance, $M$
is the total binary system mass, and $\mathcal{M}$ is the “chirp mass,”
defined as
$\mathcal{M}\equiv\frac{\left(m_{1}m_{2}\right)^{3/5}}{\left(m_{1}+m_{2}\right)^{1/5}}$.
Note here we emphasize in equation 2.23 that the amplitude or “strain” is
time-dependent, which will be important to consider later in Chapter 3. The
polarization tensors $e^{\hat{z}+}_{ij}$ and $e^{\hat{z}\times}_{ij}$ can be
conveniently decomposed into outer products of the basis vectors $\hat{x}$ and
$\hat{y}$, as shown in equation 2.24.

###  2.5 Generalized Source Orientation and Location

This solution represents a highly simplified and specific geometry - here the
orbital rotation axis is along the $\hat{z}$-direction, and the TT-gauge has
also been chosen for the gravitational waves propagating along the
$\hat{z}$-axis. However, in general we want a formalism for:

  1. 1.

a source oriented with it’s rotation axis in any arbitrary direction,

  2. 2.

and the TT-gauge to be for a gravitational wave coming from the source and
going along any arbitrary axis (not necessarily the same as the rotation
axis).

This can be achieved in two steps.

##### Step one

\- we want the same TT-gauge as before (along the $\hat{z}$-axis), but now for
a binary source in an arbitrary orientation, as shown in Fig 2.1(b). Apply two
active (“alibi”) rotation transformations to the above result - first about
the $\hat{y}$-axis by angle $\iota$ and then about the $\hat{z}$-axis by angle
$\psi$ to construct
$h^{TT\hat{z}^{{}^{\prime}}}_{ij}=R_{z}(\psi)R_{y}(\iota)h^{TT}_{ij}R_{y}^{\mathrm{Transpose}}(\iota)R_{z}^{\mathrm{Transpose}}(\psi)$.

It is important to note that we actually only need two Euler angles to fully
generalize our problem, not the usual three. There is no unique way of
defining how you rotate the axes, but in our case here let us choose to
perform the “extrinsic” rotations (about the $\hat{x}$,$\hat{y}$, and
$\hat{z}$ axes) in the following order:
$\hat{z}\rightarrow\hat{y}\rightarrow\hat{z}$ by the angles $\theta_{0}$,
$\iota$, and $\psi$, respectively. (This is one of six “proper Euler angle”
definitions). An initial rotation about the $\hat{z}$ axis by angle
$\theta_{0}$ can be equivalently achieved by beginning the binary’s orbit with
some initial orbital phase angle, which is exactly what we will later have the
equations 2.34 and 2.38 which describe the orbital motion. So in this problem
the first Euler angle is already accounted for by interpreting it as a phase
factor in the orbit, which is why we don’t need to explicitly perform this
rotation here, and why we simply perform the other two rotations by $\iota$
and $\psi$, which are the inclination and polarization angles, respectively.

The result $h^{TT\hat{z}^{{}^{\prime}}}_{ij}$ will be in the TT-gauge for a
gravitational wave traveling in the
$\hat{z}^{{}^{\prime}}=R_{z}(\psi)R_{y}(\iota)\hat{z}$ direction. To get back
the TT-gauge along the $\hat{z}$-direction, we can apply the TT-projection
operator equation 2.18 once again. Setting $\hat{n}=\hat{z}$ we write out
$\Lambda_{ijkl}h_{TT\hat{z}^{{}^{\prime}}}^{kl}=h^{TT\hat{z}}_{ij}$. The
resulting metric perturbation solution can be written compactly in a couple of
different ways:

$\displaystyle h^{TT\hat{z}}_{ij}$
$\displaystyle=e^{\hat{z}+}_{ij}H_{+}+e^{\hat{z}\times}_{ij}H_{\times}=e^{\hat{z}\textsc{A}}_{ij}H_{\textsc{A}},$
$\displaystyle=E^{\hat{z}+}_{ij}h_{+}+E^{\hat{z}\times}_{ij}h_{\times}=E^{\hat{z}\textsc{A}}_{ij}h_{\textsc{A}}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(2.25) $\displaystyle\mathrm{where}\quad\begin{cases}H_{+}\equiv
h_{+}\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\cos(2\psi)-h_{\times}\cos(\iota)\sin(2\psi),\\\
H_{\times}\equiv
h_{+}\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\sin(2\psi)+h_{\times}\cos(\iota)\cos(2\psi),\end{cases}$
(2.26) $\displaystyle\qquad\quad\ \
\begin{cases}E^{\hat{z}+}_{ij}\equiv\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)e^{\hat{z}+}_{ij}+\sin(2\psi)e^{\hat{z}\times}_{ij}\right],\\\
E^{\hat{z}\times}_{ij}\equiv\cos(\iota)\left[-\sin(2\psi)e^{\hat{z}+}_{ij}+\cos(2\psi)e^{\hat{z}\times}_{ij}\right].\end{cases}$
(2.27)

In words, we can either look at this as a transformation to the original
$h_{+}$ and $h_{\times}$ terms, or a transformation to the polarization
tensors, which can be useful interpretations.

##### Step two

\- we now want the TT-gauge to be along any arbitrary axis. Considering our
problem, we would like to be able to anchor the observer at the origin of the
coordinate system, and place the gravitational wave source anywhere in the
sky. The axis along which the gravitational wave travels which connects the
observer to the source is then the desired TT-gauge axis.

There is no unique way of doing this, but perhaps an intuitive and straight-
forward approach is as follows. Let $\hat{x}$, $\hat{y}$, and $\hat{z}$ be the
basis vectors for the observer, and let the unit vectors in spherical
coordinates be the basis vectors of the source:

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

where the angles $\theta$ and $\phi$ indicate the sky position of the source
relative to the observer. The source will be a distance $R$ away from the
observer, and the vector $\hat{r}$ points toward the source, so
$\vec{r}=R\hat{r}$.

Now we simply define a priori new polarization tensors from the outer product
of these basis vectors:

$\displaystyle\begin{cases}e^{\hat{r}+}_{ij}\equiv\hat{\theta}_{i}\hat{\theta}_{j}-\hat{\phi}_{i}\hat{\phi}_{j}\longleftrightarrow\textbf{e}^{\hat{r}+}\equiv\hat{\theta}\otimes\hat{\theta}-\hat{\phi}\otimes\hat{\phi},\\\
e^{\hat{r}\times}_{ij}\equiv\hat{\phi}_{i}\hat{\theta}_{j}+\hat{\theta}_{i}\hat{\phi}_{j}\longleftrightarrow\textbf{e}^{\hat{r}\times}\equiv\hat{\phi}\otimes\hat{\theta}+\hat{\theta}\otimes\hat{\phi}.\end{cases}$
(2.29)

Notice that equations 2.29 are actually just a generalization of equations
2.24, where we let $\hat{x}\rightarrow\hat{\theta}$,
$\hat{y}\rightarrow\hat{\phi}$, and $\hat{z}\rightarrow\hat{r}$. The
definitions give these new polarization tensors all of the required
properties. They are traceless, transverse along the $\hat{r}$-axis, and
reduce to our original $\hat{z}$-axis polarization tensors 2.24 when
$\theta=\phi=0$. The final step is a simply replacement. Since these
polarization tensors 2.29 are a generalization of the original polarization
tensors 2.24, we can substitute these new tensors into the expressions in the
metric perturbation from before, 2.25 and 2.27.

The final fully general result for our problem of interest can now be
expressed compactly in a couple of different ways:

$\displaystyle h^{TT\hat{r}}_{ij}$
$\displaystyle=e^{\hat{r}+}_{ij}H_{+}+e^{\hat{r}\times}_{ij}H_{\times}=e^{\hat{r}\textsc{A}}_{ij}H_{\textsc{A}},$
$\displaystyle=E^{\hat{r}+}_{ij}h_{+}+E^{\hat{r}\times}_{ij}h_{\times}=E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}}\qquad\mathrm{for}\quad\textsc{A}\in[+,\times],$
(2.30) $\displaystyle\mathrm{where}\quad\begin{cases}H_{+}\equiv
h_{+}\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\cos(2\psi)-h_{\times}\cos(\iota)\sin(2\psi),\\\
H_{\times}\equiv
h_{+}\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\sin(2\psi)+h_{\times}\cos(\iota)\cos(2\psi),\end{cases}$
(2.26 r) $\displaystyle\qquad\quad\ \
\begin{cases}E^{\hat{r}+}_{ij}\equiv\frac{1}{2}\left(1+\cos^{2}(\iota)\right)\left[\cos(2\psi)e^{\hat{r}+}_{ij}+\sin(2\psi)e^{\hat{r}\times}_{ij}\right],\\\
E^{\hat{r}\times}_{ij}\equiv\cos(\iota)\left[-\sin(2\psi)e^{\hat{r}+}_{ij}+\cos(2\psi)e^{\hat{r}\times}_{ij}\right],\end{cases}$
(2.31) $\displaystyle\qquad\quad\ \
\begin{cases}h_{+}(t)\equiv-h(t)\cos\big{(}2\Theta(t)\big{)},\\\
h_{\times}(t)\equiv-h(t)\sin\big{(}2\Theta(t)\big{)},\end{cases}$ (2.22 r)
$\displaystyle\hskip
51.21504pth(t)\equiv\frac{4(G\mathcal{M})^{5/3}}{c^{4}R}\omega(t)^{2/3},$
(2.23 r) $\displaystyle\qquad\quad\ \
\begin{cases}e^{\hat{r}+}_{ij}\equiv\hat{\theta}_{i}\hat{\theta}_{j}-\hat{\phi}_{i}\hat{\phi}_{j},\\\
e^{\hat{r}\times}_{ij}\equiv\hat{\phi}_{i}\hat{\theta}_{j}+\hat{\theta}_{i}\hat{\phi}_{j},\end{cases}$
(2.29 r) $\displaystyle\qquad\quad\ \ \begin{cases}\begin{tabular}[]{c l c
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{]},$\end{tabular}\end{cases}$
(2.28 r)

As was the case before, we can either look at this as a transformation to the
original $h_{+}$ and $h_{\times}$ terms, or a transformation to the
polarization tensors. We can test that these expressions are indeed a
generalization of the original problem with all of the required properties.
The metric perturbation is finally in the TT-gauge now for a gravitational
wave propagating along the $\hat{r}$-axis. It is transverse to the direction
of propagation of the gravitational wave and traceless, both of which follow
since the polarization tensors are transverse and traceless, and the metric
perturbation directly decomposes into a linear combination of these tensors.
The sky angles $\theta$ and $\phi$ point toward the location of the source,
and the angles $\iota$ and $\psi$ indicate the source’s own orientation. All
of this is shown in Figure 2.2. Furthermore, our motivation for this approach
and these definitions stems from our comfort and familiarity with spherical
coordinate unit basis vectors.

(a) Fully generalized geometrical positioning of a gravitational wave source
with respect to an observer. Angles $\theta$ and $\phi$ denote source sky
position, and $\iota$ and $\psi$ indicate source inclination and polarization.
Source appears “face-on” if $\iota=0$, and “edge-on” if $\iota=\pi/2$, and
“face-off” if $\iota=\pi$.

(b) The source basis vectors are chosen so that when $\theta=\phi=0$, then
$\hat{r}=\hat{z}$, $\hat{\theta}=\hat{x}$, and $\hat{\phi}=\hat{y}$. In this
case, if the source is not inclined (i.e. $\iota=0$), then the observer will
see the binary orbiting in the counter-clockwise direction.

Figure 2.2:

###  2.6 Source Parameters and a Word about Notation

We take a moment here simply to point out the physical significance behind the
notation we have employed so far, as it will greatly simplify our work later
and lead to solutions which are written in a way that are easier to understand
physically.

A circular binary system produces gravitational waves which in the zeroth
order Newtonian regime is described by a total of eight parameters:

$\mathrm{Source}\
\mathrm{Parameters:}\quad\\{R,\theta,\phi,\iota,\psi,\theta_{0},\mathcal{M},\omega_{0}\\}\quad\longrightarrow\quad\begin{cases}\mathrm{Extrinsic:}\quad&\\{R,\theta,\phi,\iota,\psi,\theta_{0}\\}\\\
\mathrm{Intrinsic:}\quad&\\{\mathcal{M},\omega_{0}\\}\end{cases}$

(Note that $\theta_{0}$ and $\omega_{0}$ will be shown explicitly in the next
section 2.7, but they appear in $\Theta(t)$ and $\omega(t)$ above).

Of these eight parameters, the angles $\theta$, $\phi$, $\iota$, and $\psi$
are the binary’s “geometrical orientation” angles - they describe the
orientation of the orbit of the binary, and the binary’s location on the sky
in reference to the Earth’s frame. The parameter $\theta_{0}$ we’ll think of
as the initial phase of the binary, although it was pointed out that this is
also the third Euler angle for the binary’s orientation. The rest of the
parameters - $R$, $\mathcal{M}$, $\theta_{0}$, and $\omega_{0}$ \- influence
only the strain $h$ of the gravitational wave. Lastly, it should be pointed
out too that the strain is also influenced by the $\theta$ and $\phi$, because
the strain is evaluated at the retarded time, and the retarded time will
depend on these two parameters (as discussed in Section 3.3).

In terms of notation, as we mentioned before we can look at the metric
perturbation in equation 2.30 in either of two ways. In this paper we prefer
to use the notation
$h^{TT\hat{r}}_{ij}=E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}}$ because it
groups all of the geometrical orientation and location angles into the
definition of the polarization tensor, and keeps the strain in the simple form
that we started with in equation 2.22. We will point out that many sources
seem to prefer a notation which more closely resembles the alternative form of
writing $h^{TT\hat{r}}_{ij}=e^{\hat{r}\textsc{A}}_{ij}H_{\textsc{A}}$ (see for
example, Corbin & Cornish, 2010; Zhu et al., 2014; Aggarwal et al., 2019).
Perhaps the only disadvantage of this is that while the polarization tensor is
a function only of the sky position angles, the strain $H_{\textsc{A}}$
becomes a mixture of both $h_{+}$ and $h_{\times}$. But since:

$e^{\hat{r}\textsc{A}}_{ij}H_{\textsc{A}}=E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}},$

we can use these two notations interchangeably as is most convenient.

Also as an aside, we point out that when reviewing pulsar timing literature,
one should pay close attention to the use of notation and definitions of
quantities, namely the vectors that go into building the polarization tensors.
Often one will find the use of the vectors $\hat{\Omega}$, $\hat{m}$, and
$\hat{n}$ for the source basis vectors (instead of $\hat{r}$, $\hat{\theta}$,
and $\hat{\phi}$ as we have used. What becomes tricky and sometimes
frustrating, however, is that these vectors won’t always be defined in the
same way. Typically they differ in the direction they are pointed, i.e. by a
$\pm$ sign (see for example, Arzoumanian et al., 2014; Chamberlin et al.,
2015; Taylor et al., 2020). Noting these differences is crucial, as they will
necessarily introduce $\pm$ signs in the definitions of the polarization
tensors and antenna patterns (see Section 3.5.1).

Mathematically this is perfectly fine, in the sense that all of the math is
still consistent when working within a given choice of notation/definitions.
However, the conceptualization quantities such as angular parameters between
notations may differ. This is why we propose and explain the motivation of the
notation and conventions established here in our work, as they simply adopt
the often used and familiar spherical coordinate basis vectors as their
building point.

###  2.7 Two Models: Monochromatism and Frequency Evolution

There are two important models used for computing the evolution of the metric
perturbation and the timing residual. In short, the frequency evolution model
is a more generalized version of the timing residual, and reduces to the
monochromatic model in the appropriate limit (described below).

Here we call attention to $t_{0}$, the “fiducial time” for our two models
defined below. As the model’s reference time, we define $t_{0}$ as the time at
which $\theta_{0}$ and $\omega_{0}$ are measured:

$\displaystyle\Theta(t_{0})\equiv\theta_{0},$
$\displaystyle\omega(t_{0})\equiv\omega_{0},$ (2.32)

We can choose any time to be the fiducial time of our model, but perhaps the
two most comfortable choices are either $t_{0}=0$, or $t_{0}=-\frac{R}{c}$.
Conceptually choosing $t_{0}=0$ means these parameters are the Earth frame
“present-day” values of the orbital phase and frequency of the binary (when
our experiment first begins), and choosing $t_{0}=-\frac{R}{c}$ are the
“retarded-day” values, that is the values of these parameters that produced
the gravitational waves which are currently arriving at the Earth (again, at
the time our experiment first begins). Later in this work we will choose the
fiducial time to be $t_{0}=-\frac{R}{c}$.

####  2.7.1 Monochromatism

In the most simple model of the gravitational waves emitted by a binary
system, one assumes no energy is lost from the emission of gravitational
waves. Under this assumption, the binary’s energy will remain conserved, the
orbit will not coalesce (it has “infinite” coalescence time), and hence the
orbital period/frequency will remain constant. This assumption means that the
gravitational waves emitted by the binary system will have a constant
frequency for all time. Hence the name “monochromatic” gravitational waves.

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

$\displaystyle\omega(t)=$ $\displaystyle\omega_{0},$ (2.33)
$\displaystyle\Theta(t)=$ $\displaystyle\theta_{0}+\omega_{0}(t-t_{0}).$
(2.34)

####  2.7.2 Frequency Evolution

More realistically, however, energy in the binary system is lost due to the
emission of gravitational wave radiation. This means that over time the binary
will collapse and coalesce. As the binary collapses, the orbital frequency
will increase, meaning that the gravitational wave frequency will also
increase. This evolution towards higher frequencies is called “frequency
chirping.”

Under assumption 5, gravitational wave radiation from a binary system causes
the frequency of the orbit to change as:

$\dot{\omega}=\frac{d\omega}{dt}=\frac{96}{5}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{5/3}\omega^{11/3}.$
(2.35)

As the orbit coalesces, the orbital frequency increases without bound.
Formally we define the time the binary coalesces to be when it reaches
infinite orbital frequency, that is $\omega(t)\rightarrow\infty$ as
$t\rightarrow\tau_{c}$. We can compute the time of coalescence of the binary
by integrating equation 2.35:

$\int\limits^{\tau_{c}}_{t_{0}}dt=\int\limits^{\infty}_{\omega_{0}}\frac{5}{96}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\omega^{-11/3}d\omega,$

which gives:

$\displaystyle\tau_{c}$ $\displaystyle=t_{0}+\Delta\tau_{c},$
$\displaystyle\mathrm{where}\quad\Delta\tau_{c}\equiv\frac{5}{256}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\frac{1}{\omega_{0}^{8/3}}.$
(2.36)

Here in our notation $\tau_{c}$ is the “time of coalescence” and
$\Delta\tau_{c}\equiv\tau_{c}-t_{0}$ is the “time of coalescence measured from
our chosen fiducial time.” We make the pedantic distinction of this here as
the difference between the two will become more significant later on. This
difference is represented in Figure 3.7, and the dependence of
$\Delta\tau_{c}$ on the chirp mass and orbital frequency parameters can be
seen in the contour plot in Figure 3.9.

For arbitrary times, we can repeat the integration:

$\int\limits^{t}_{t_{0}}dt=\int\limits^{\omega}_{\omega_{0}}\frac{5}{96}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\omega^{{}^{\prime}\
-11/3}d\omega^{{}^{\prime}},$

which leads us to an expression for the orbital frequency as a function of
time for our binary system:

$\omega(t)=\omega_{0}\left[1-\frac{t-t_{0}}{\Delta\tau_{c}}\right]^{-3/8}.$
(2.37)

As for the phase of the binary, we know $\frac{d\Theta}{dt}=\omega(t)$. We can
integrate this expression over time, or similary we can first write
$\omega(t)=\frac{d\Theta}{d\omega}\frac{d\omega}{dt}$, so using equation 2.35
again we can write
$\frac{d\Theta}{d\omega}=\frac{5}{96}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\omega^{-8/3}$
and integrate this for arbitary $\omega$:

$\int\limits^{\Theta}_{\theta_{0}}d\Theta^{{}^{\prime}}=\int\limits^{\omega}_{\omega_{0}}\frac{5}{96}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\omega^{{}^{\prime}\
-8/3}d\omega^{{}^{\prime}},$

which leads us to an expression for the orbital phase as a function of time
for our binary system (expressed here in a number of useful forms):

$\displaystyle\Theta(t)$
$\displaystyle=\theta_{0}+\theta_{c}\left[1-\left(\frac{\omega(t)}{\omega_{0}}\right)^{-5/3}\right],$
$\displaystyle=\theta_{0}+\theta_{c}\left[1-\left(1-\frac{t-t_{0}}{\Delta\tau_{c}}\right)^{5/8}\right],$
$\displaystyle=\theta_{0}+\frac{1}{32}\left(\frac{c^{3}}{G\mathcal{M}}\right)^{5/3}\left[\omega_{0}^{-5/3}-\omega(t)^{-5/3}\right],$
(2.38)

where we define the quantity:

$\theta_{c}\equiv\frac{8}{5}\Delta\tau_{c}\omega_{0}=\frac{1}{32}\left(\frac{c^{3}}{G\mathcal{M}\omega_{0}}\right)^{5/3}.$
(2.39)

This is a rather important quantity itself, and conceptually it is the number
of radians that the binary will sweep out before it coalesces (i.e.
$\frac{\theta_{c}}{2\pi}$ is the number of orbital revolutions before the
system coalesces). That is, in the limit the binary coalesces where
$\omega(t)\rightarrow\infty$, we see that
$\Theta\rightarrow\theta_{0}+\theta_{c}$. For this reason we will refer to
$\theta_{c}$ as the “coalescence angle.” Remember, both $\Delta\tau_{c}$ and
$\theta_{c}$ are measured with respect to the chosen fiducial time.

Note that in the “large coalescence time” limit that
$\Delta\tau_{c}\rightarrow\infty$, the frequency evolution model (equations
2.37 and 2.38) reduces to the monochromatic model (equations 2.33 and 2.34) as
we would expect. In practice, the monochromatic limit for pulsar timing
experiments merely requires that the gravitational wave frequency stays nearly
constant as pulsar’s radiation travels between the pulsar and the Earth, so
that $\Delta\tau_{c}\gg\left(1-\hat{r}\cdot\hat{p}\right)\frac{L}{c}$ (see
ahead to equation 3.48). However in this work since we are specifically
investigating and classifying different models based on their mathematical
differences, we will use the term “monochromatic” specifically when working
under the assumption of Section 2.7.1.

Finally, using this information we can also express equation 2.35 as:

$\displaystyle\dot{\omega}$
$\displaystyle=\frac{d\omega}{dt}=\frac{96}{5}\left(\frac{G\mathcal{M}}{c^{3}}\right)^{5/3}\omega^{11/3},$
$\displaystyle=\frac{3}{8}\frac{\omega^{11/3}}{\Delta\tau_{c}\omega_{0}^{8/3}},$
$\displaystyle=\frac{3}{8}\frac{\omega}{\tau_{c}-t}.$ (2.40)

##  Chapter 3 The Continuous Wave Timing Residual

A binary system of two massive objects orbiting each other will, under general
relativity, cause the metric perturbation to vary in time sinusoidally, and
this perturbation will propagate away from the source as a wave. The result is
that spacetime will stretch and compress as the wave propagates through the
universe, and will change the path (geodesic) that external objects move along
as the wave interferes with them. For gravitational wave experiments, it is
especially important to try and understand how the gravitational wave will
affect null geodesics, i.e. the path that light moves along, because we can
design experiments that can actually test and observe this effect on light.

With this in mind we now turn to the question of using pulsars to design a
gravitational wave experiment. When observed with radio telescopes here on
Earth, highly regularly spinning pulsars act as incredibly precise clocks,
with each pulse acting as a clock tick. Now consider the effect of a
gravitational wave on the photons in a pulse of light beamed from a pulsar, as
it travels from the pulsar towards our observatories here on Earth. In the
absence of any gravitational wave disturbances, the photons would move along
an essentially flat spacetime background as they traveled from the pulsar to
the Earth. Since the spacetime the photons traverse is static, the interval
between the light pulses will remain constant, or in other words, the pulsar’s
pulse period $T$ will remain constant as observed on Earth. However, a
gravitational wave will cause that flat spacetime to perturb sinusoidally, and
in this way the photon spacetime path from the pulsar to the Earth will
change. The question we want to ask is how is the observed period of the
pulsar affected by a gravitational wave, and can we find an observable that
allows us to turn this theoretical observation into an experimental one?

For the work below we will make the following assumptions:

Assumptions: Observed Pulsar Period 1. The universe is flat and static.
Therefore the expression for the retarded time is the familiar:
$t_{\mathrm{ret}}=t-\frac{|\vec{x}-\vec{x}^{{}^{\prime}}|}{c}$ . 2. Earth is
at the center of our coordinate system. 3. The pulsar is located at a fixed
angular sky position and distance from the Earth. 4. The local background
spacetime between the pulsar and the Earth is flat Minkowski,
${\eta_{\mu\nu}=\mathrm{diag}(-c^{2},\ 1,\ 1,\ 1)}$. 5. Transverse-traceless
(“TT”) gauge. 6. The antenna patterns are assumed to remain constant over the
Earth-pulsar baseline. 7. The pulsar’s rotational period timescale must be
appropriately small in comparison to the timescale of the gravitational wave.
Specifically we have two case requirements: $\Rightarrow$ Monochromatic
Source: $\omega_{0}T\ll 1$ , $\Rightarrow$ Frequency Evolving Source:
$\frac{T}{\Delta\tau_{c}}\ll 1$ .

###  3.1 The Observed Pulsar Period

So the first question we ask is, what is the path of a photon traveling from
the pulsar to the Earth?111The derivation presented here follows the logic
presented in Maggiore (2018). To determine this we begin with the spacetime
interval and our metric, a flat background spacetime with a gravitational wave
perturbation coming from some binary source, as given by equation 2.30:

$ds^{2}=-(cdt)^{2}+\left[\eta_{ij}+h^{TT\hat{r}}_{ij}\right]dx^{i}dx^{j}.$
(3.1)

It’s important to remember here that the metric perturbation
$h^{TT\hat{r}}_{ij}$ is a function of the retarded time, which itself is a
function of $t$ and $\vec{x}$, that is, our notation for the functional
dependence on $t_{\mathrm{ret}}$, $t$, and $\vec{x}$ can be indicated as:

$h^{TT\hat{r}}_{ij}\equiv
h^{TT\hat{r}}_{ij}\left(t_{\mathrm{ret}}\right)\equiv
h^{TT\hat{r}}_{ij}\big{(}t_{\mathrm{ret}}(t,\vec{x})\big{)}\equiv
h^{TT\hat{r}}_{ij}\left(t,\vec{x}\right).$ (3.2)

Remembering this order, that the metric perturbation is a function of the
retarded time, which itself depends on the field point of interest $\vec{x}$,
will become important shortly, which is why we make special emphasis of it
here before continuing.

Let the pulsar be at a fixed angular position given by $\theta_{p}$ and
$\phi_{p}$; the vector pointing from the Earth towards the pulsar will be
denoted as $\hat{p}$:

$\hat{p}\equiv\big{[}\sin(\theta_{p})\cos(\phi_{p}),\
\sin(\theta_{p})\sin(\phi_{p}),\ \cos(\theta_{p})\big{]}$ (3.3)

The photon will travel along the radial spatial path connecting the pulsar and
the Earth (hence $d\theta=d\phi=0$), so using spherical coordinates we can
write $dx^{i}=\hat{p}^{i}dr$. Moreover, since we are considering the motion of
a photon, $ds=0$. So equation 3.1 becomes:

$\displaystyle(cdt)^{2}$
$\displaystyle=\left[\eta_{ij}+h^{TT\hat{r}}_{ij}\right]\hat{p}^{i}\hat{p}^{j}dr^{2},$
$\displaystyle=\left[1+\hat{p}^{i}\hat{p}^{j}h^{TT\hat{r}}_{ij}\right]dr^{2},$
$\displaystyle=\left[1+\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}}\right]dr^{2}\quad\quad\longrightarrow\quad\pm\
dr=\frac{cdt}{\left[1+\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}}\right]^{1/2}},$
(3.4)

where the second line just comes from writing out the Einstein notation:
$\hat{p}^{i}\hat{p}^{j}\eta_{ij}=(\hat{p}^{x})^{2}+(\hat{p}^{y})^{2}+(\hat{p}^{z})^{2}=|\hat{p}|^{2}=1$.
(Or alternatively, we can write the spatial part of the Minkowski metric as a
decomposed sum of the outer products of the unit Cartesian vectors
$\eta_{ij}\equiv\left(\hat{x}_{i}\hat{x}_{j}\right)+\left(\hat{y}_{i}\hat{y}_{j}\right)+\left(\hat{z}_{i}\hat{z}_{j}\right)$
and then recognize that this is just the dot product of $\hat{p}$ with itself:
$\hat{p}^{i}\hat{p}^{j}\eta_{ij}=|\hat{p}|^{2}=1$). In the third line we
switch to the notation established and motivated in equation 2.30, and note
that here the quantity $\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}$ is
function only of the source’s orientation and sky location angles, and
importantly, doesn’t have any time dependence. We will find later that this
combines with another term to become what we will call the “antenna response”
function (discussed in detail in Section 3.5.1).

The $\pm$ from the square root denotes a radially outbound ($+$) versus
radially inbound ($-$) photon. We now integrate both sides of this expression
to get the photon’s path. The photon leaves the pulsar at the “emitted” time
$t_{\mathrm{em}}$ at a distance $L$ from the Earth, and travels radially
inbound arriving at the Earth at the “observed” time $t_{\mathrm{obs}}$:

$\displaystyle
L=-\int\limits^{0}_{L}dr=\int\limits^{t_{\mathrm{obs}}}_{t_{\mathrm{em}}}\frac{cdt}{\left[1+\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}}\right]^{1/2}}$
$\displaystyle\approx\int\limits^{t_{\mathrm{obs}}}_{t_{\mathrm{em}}}\left[1-\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}h_{\textsc{A}}\right]cdt,$
$\displaystyle=c\left(t_{\mathrm{obs}}-t_{\mathrm{em}}\right)-\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{obs}}}_{t_{\mathrm{em}}}h_{\textsc{A}}\left(t,\vec{x}\right)cdt.$

We can see in the final equality here that to “zeroth order” in the metric
perturbation the path that the photon takes is $L\approx
c\left(t_{\mathrm{obs}}-t_{\mathrm{em}}\right)$, or $t_{\mathrm{obs}}\approx
t_{\mathrm{em}}+L/c$, which is exactly as we would expect for flat unperturbed
spacetime. Since we are only interested in the solution to first order in the
metric perturbation, we can use this to replace the upper limit of integration
and then evaluate the metric perturbation along the zeroth order path (the
overall integrand is still first order). We can also use this to write the
spatial path of our photon to zeroth order as:

$\vec{x}_{0}(t)\equiv\left[ct_{\mathrm{em}}+L-ct\right]\hat{p}\approx\left[ct_{\mathrm{obs}}-ct\right]\hat{p}.$
(3.5)

This describes the desired path - i.e. $\vec{x}_{0}=L\hat{p}$ at the time the
photon leaves the pulsar and $\vec{x}_{0}=0$ at the time it arrives at the
Earth. Now we can write:

$t_{\mathrm{obs}}\approx
t_{\mathrm{em}}+\frac{L}{c}+\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{em}}+L/c}_{t_{\mathrm{em}}}h_{\textsc{A}}\big{(}t,\vec{x}_{0}(t)\big{)}dt.$
(3.6)

So we have considered the path that a photon will travel between the pulsar
and the Earth, and have found an expression equation 3.6 which gives the
difference in emitted and observed times of the photon. Next imagine the
pulsar emits a flash, then rotates once and emits a second flash. We’ll denote
the true period of the pulsar as $T$, but now want to find what is the
observed period here on Earth. We repeat the exact same steps above for a
photon one period later. Now the emitted time is
$t_{\mathrm{em}}^{{}^{\prime}}=t_{\mathrm{em}}+T$ and the observed time is
$t_{\mathrm{obs}}^{{}^{\prime}}$. The the spatial path then becomes:

$\displaystyle\vec{x}^{{}^{\prime}}_{0}(t)$
$\displaystyle\equiv\vec{x}_{0}(t,t_{\mathrm{em}}\rightarrow
t_{\mathrm{em}}^{{}^{\prime}})=\left[ct_{\mathrm{em}}+cT+L-ct\right]\hat{p}=\left[ct_{\mathrm{em}}+L-c(t-T)\right]\hat{p},$
$\displaystyle=\vec{x}_{0}(t-T).$ (3.7)

The expression one period later changes to:

$\displaystyle t_{\mathrm{obs}}^{{}^{\prime}}$ $\displaystyle\approx
t_{\mathrm{em}}+T+\frac{L}{c}+\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{em}}+T+L/c}_{t_{\mathrm{em}+T}}h_{\textsc{A}}\big{(}t,\vec{x}_{0}(t-T)\big{)}dt,$
$\displaystyle=t_{\mathrm{em}}+T+\frac{L}{c}+\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{em}}+L/c}_{t_{\mathrm{em}}}h_{\textsc{A}}\big{(}t+T,\vec{x}_{0}(t)\big{)}dt,$
(3.8)

where in the second line we performed the substitution/coordinate-shift
$\tilde{t}=t-T$. This shifts the limits of integration, and since the label of
the coordinate of integration in an integral is arbitrary, we can let
$\tilde{t}\rightarrow t$. The integrals in expressions 3.6 and 3.8 now only
differ by a factor of $T$ in the time-coordinate, but not in the spatial path.
So the difference between equation 3.6 and 3.8 is the observed pulsar pulse
period here on Earth:

$\displaystyle T_{\mathrm{obs}}\equiv
t^{{}^{\prime}}_{\mathrm{obs}}-t_{\mathrm{obs}}$ $\displaystyle\approx
T+\Delta T,$ $\displaystyle\mathrm{where}\quad\Delta T$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{em}}+L/c}_{t_{\mathrm{em}}}\bigg{[}h_{\textsc{A}}\big{(}t+T,\vec{x}_{0}(t)\big{)}-h_{\textsc{A}}\big{(}t,\vec{x}_{0}(t)\big{)}\bigg{]}dt,$
(3.9)

with equation 3.9 being the general expression for the resulting perturbation
of a gravitational wave on the observed pulsar period.

In Section 3.2 we will explain that a more useful quantity to consider is the
fractional change in the pulsar’s period, so dividing equation 3.9 by $T$ we
write:

$\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{em}}+L/c}_{t_{\mathrm{em}}}\left[\frac{h_{\textsc{A}}\big{(}t+T,\vec{x}_{0}(t)\big{)}-h_{\textsc{A}}\big{(}t,\vec{x}_{0}(t)\big{)}}{T}\right]dt,$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{em}}+L/c}_{t_{\mathrm{em}}}\left[\frac{h_{\textsc{A}}\big{(}t+T,\vec{x}\big{)}-h_{\textsc{A}}\big{(}t,\vec{x}\big{)}}{T}\right]\Bigg{\rvert}_{\vec{x}=\vec{x}_{0}(t)}dt,$
$\displaystyle\approx\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{\mathrm{em}}+L/c}_{t_{\mathrm{em}}}\frac{\partial
h_{\textsc{A}}\left(t,\vec{x}\right)}{\partial
t}\Bigg{\rvert}_{\vec{x}=\vec{x}_{0}(t)}dt.$ (3.10)

The approximation of the integrand that we made in the final line is crucial,
and makes use of the definition of a derivative. However, to arrive at this
approximation we must consider the two fundamentally different cases of
interest. Regardless if the source is truly monochromatic or if its frequency
is evolving, the result here is the same, but the approximation statements
will be different as shown in the box below.

Monochromatic Source
For a monochromatic source we see that from the functional form of the phase,
frequency, and strain (Section 2.7.1 and equation 2.23), $t$ will always
appear in the combination $\omega_{0}t$, a dimensionless quantity, in our
equations. Let’s use this to non-dimensionalize our expression so that we can
properly take the limit definition of our derivative. It is tempting to take
the limit definition of our derivative of the integrand in equation 3.10 in
terms of $T$ alone, but this is a dimensional quantity, and therefore this
statement would not be sufficiently rigorous. So let:
$\tilde{t}\equiv\omega_{0}t\quad\implies\tilde{T}=\omega_{0}T.$ This means we
can write:
$\displaystyle\frac{h_{\textsc{A}}\big{(}t+T,\vec{x}\big{)}-h_{\textsc{A}}\big{(}t,\vec{x}\big{)}}{T}$
$\displaystyle=\frac{h_{\textsc{A}}\big{(}\tilde{t}+\tilde{T},\vec{x}\big{)}-h_{\textsc{A}}\big{(}\tilde{t},\vec{x}\big{)}}{T}\times\frac{\omega_{0}}{\omega_{0}},$

$\displaystyle=\left[\frac{h_{\textsc{A}}\big{(}\tilde{t}+\tilde{T},\vec{x}\big{)}-h_{\textsc{A}}\big{(}\tilde{t},\vec{x}\big{)}}{\tilde{T}}\right]\omega_{0},$
$\displaystyle\underset{\tilde{T}\ll 1}{\approx}\left[\frac{\partial
h_{\textsc{A}}\big{(}\tilde{t}\big{)}}{\partial\tilde{t}}\right]\omega_{0},$
$\displaystyle=\left[\frac{\partial t}{\partial\tilde{t}}\frac{\partial
h_{\textsc{A}}\big{(}t\big{)}}{\partial t}\right]\omega_{0},$
$\displaystyle=\left[\frac{1}{\omega_{0}}\frac{\partial
h_{\textsc{A}}\big{(}t\big{)}}{\partial t}\right]\omega_{0},$
$\displaystyle=\frac{\partial h_{\textsc{A}}\big{(}t\big{)}}{\partial t}.$ The
first line is true because again, $t$ and $T$ appear everywhere in our
expression for $h_{\textsc{A}}$ in the combinations $\omega_{0}t$ and
$\omega_{0}T$. In the third line we have applied our derivative definition
approximation, namely that in the limit $\tilde{T}\rightarrow 0$ the
approximation becomes an exact equality. And in the fourth line we apply the
chain rule to write the derivative with respect to $t$. So the crucial
assumption for a monochromatic source is that we require
$\tilde{T}=\omega_{0}T\ll 1$. For our sources of interest, SMBHBs, the orbital
period will be on the order of months to decades, so the orbital angular
frequency $\omega_{0}$ will be on the order of $10^{-7}$ to $10^{-9}$ Hz. And
the pulsar rotation period will be on the order of milliseconds. So
$\omega_{0}T\sim\mathcal{O}(10^{-10})$ or so, which indeed will be a very
small quantity, so this approximation is safe. Frequency Evolving Source
For a source whose orbit evolves with time, we see that from the functional
form of the phase, frequency, and strain (Section 2.7.2 and equation 2.23),
$t$ will always

appear in the combination $\frac{t}{\Delta\tau_{c}}$ in our equations, which
is again a dimensionless quantity. For the same reasons as in the
monochromatic cases, let:
$\tilde{t}\equiv\frac{t}{\Delta\tau_{c}}\quad\implies\tilde{T}=\frac{T}{\Delta\tau_{c}}.$
Once again, we carry through the same series of steps:
$\displaystyle\frac{h_{\textsc{A}}\big{(}t+T,\vec{x}\big{)}-h_{\textsc{A}}\big{(}t,\vec{x}\big{)}}{T}$
$\displaystyle=\frac{h_{\textsc{A}}\big{(}\tilde{t}+\tilde{T},\vec{x}\big{)}-h_{\textsc{A}}\big{(}\tilde{t},\vec{x}\big{)}}{T}\times\frac{1/\Delta\tau_{c}}{1/\Delta\tau_{c}},$
$\displaystyle=\left[\frac{h_{\textsc{A}}\big{(}\tilde{t}+\tilde{T},\vec{x}\big{)}-h_{\textsc{A}}\big{(}\tilde{t},\vec{x}\big{)}}{\tilde{T}}\right]\frac{1}{\Delta\tau_{c}},$
$\displaystyle\underset{\tilde{T}\ll 1}{\approx}\left[\frac{\partial
h_{\textsc{A}}\big{(}\tilde{t}\big{)}}{\partial\tilde{t}}\right]\frac{1}{\Delta\tau_{c}},$
$\displaystyle=\left[\frac{\partial t}{\partial\tilde{t}}\frac{\partial
h_{\textsc{A}}\big{(}t\big{)}}{\partial t}\right]\frac{1}{\Delta\tau_{c}},$
$\displaystyle=\left[\Delta\tau_{c}\frac{\partial
h_{\textsc{A}}\big{(}t\big{)}}{\partial t}\right]\frac{1}{\Delta\tau_{c}},$
$\displaystyle=\frac{\partial h_{\textsc{A}}\big{(}t\big{)}}{\partial t}.$ So
the crucial assumption for a source whose frequency does evolve over time is
that we require $\tilde{T}=\frac{T}{\Delta\tau_{c}}\ll 1$. Again for our
sources of interest, SMBHBs, the time of coalescence (measured from the
fiducial time) will typically be on the order of kiloyears all of the way to
gigayears. Combined with a millisecond pulsar period this requirement again
will normally be safely met. This might only become a problem if we are
considering an agressively evolving source, wherein we actually capture the
moment of coalescence. If this is the case then this assumption would break
down.

Next, let’s consider the functional form of this integrand expression in
equation 3.10. From equations 2.22, 2.23, and 3.2 the functional dependence of
the metric perturbation on time looks like:
$h_{\textsc{A}}=h_{\textsc{A}}\Bigg{(}h\Big{(}\omega\big{(}t_{\mathrm{ret}}(t,\vec{x})\big{)}\Big{)},\Theta\big{(}t_{\mathrm{ret}}(t,\vec{x})\big{)}\Bigg{)}$.
Furthermore, under assumption 1, $\frac{\partial t_{\mathrm{ret}}}{\partial
t}=1$. Therefore taking the time derivative of the metric perturbation, the
chain rule yields:

$\displaystyle\frac{\partial
h_{\textsc{A}}\left(t_{\mathrm{ret}}\left(t,\vec{x}\right)\right)}{\partial
t}$
$\displaystyle=\frac{dh_{\textsc{A}}\left(t_{\mathrm{ret}}\right)}{dt_{\mathrm{ret}}}\frac{\partial
t_{\mathrm{ret}}}{\partial
t}=\frac{dh_{\textsc{A}}\left(t_{\mathrm{ret}}\right)}{dt_{\mathrm{ret}}}=\frac{dh_{\textsc{A}}\left(h\left(\omega\left(t_{\mathrm{ret}}\right)\right),\Theta\left(t_{\mathrm{ret}}\right)\right)}{dt_{\mathrm{ret}}},$
$\displaystyle=\frac{\partial h_{\textsc{A}}(h,\Theta)}{\partial
h}\frac{\partial
h(\omega)}{\partial\omega}\frac{d\omega}{dt_{\mathrm{ret}}}+\frac{\partial
h_{\textsc{A}}(h,\Theta)}{\partial\Theta}\frac{d\Theta}{dt_{\mathrm{ret}}},$
$\displaystyle=\left(\frac{h_{\textsc{A}}}{h}\right)\left(\frac{2}{3}\frac{h}{\omega}\right)\dot{\omega}(t_{\mathrm{ret}})+\bigg{(}2h_{\textsc{A}}\left(h,\Theta+\frac{\pi}{4}\right)\bigg{)}\omega(t_{\mathrm{ret}}),$
$\displaystyle=\frac{2}{3}h_{\textsc{A}}(t_{\mathrm{ret}})\frac{\dot{\omega}(t_{\mathrm{ret}})}{\omega(t_{\mathrm{ret}})}+2h_{\textsc{A}}\left(\omega(t_{\mathrm{ret}}),\Theta(t_{\mathrm{ret}})+\frac{\pi}{4}\right)\omega(t_{\mathrm{ret}}),$
(3.11)

where in the third line here we have used the notation described in the
“Notational Aside” box. The “dot” derivative $\dot{\omega}$ here is a
derivative with respect to $t_{\mathrm{ret}}$.

We see here by equation 3.11 that if we carefully consider the order of the
derivatives and then evaluate the expression along the photon path
$\vec{x}=\vec{x}_{0}(t)$, we can say (again see the “Notational Aside” box
below):

$\displaystyle\frac{\partial
h_{\textsc{A}}\left(t_{\mathrm{ret}}\left(t,\vec{x}\right)\right)}{\partial
t}\Bigg{\rvert}_{\vec{x}=\vec{x}_{0}(t)}$ $\displaystyle=\ \
\frac{dh_{\textsc{A}}\left(t_{\mathrm{ret}}\right)}{dt_{\mathrm{ret}}}\Bigg{\rvert}_{t_{\mathrm{ret}}=t_{\mathrm{ret}}^{0}},$
$\displaystyle=\ \
\frac{2}{3}h_{\textsc{A}}(t^{0}_{\mathrm{ret}})\frac{\dot{\omega}(t^{0}_{\mathrm{ret}})}{\omega^{0}}+2h_{\textsc{A}}\left(\omega^{0},\Theta^{0}+\frac{\pi}{4}\right)\omega^{0},$
$\displaystyle=\ \
\frac{dh_{\textsc{A}}\left(t^{0}_{\mathrm{ret}}\right)}{dt^{0}_{\mathrm{ret}}}.$
(3.12)

In words, taking the time derivative of $h_{\textsc{A}}=h_{\textsc{A}}(t)$ and
then evaluating the result at $\vec{x}=\vec{x}_{0}(t)$ is equivalent to just
writing $h_{\textsc{A}}=h_{\textsc{A}}\left(t^{0}_{\mathrm{ret}}\right)$ and
taking a derivative of that with respect to $t^{0}_{\mathrm{ret}}$. This will
end up being a very important result later in Section 3.5.

With each of these established notations we can write the integrand of the
fractional change in the pulsar period in a number of useful forms, all of
which we will exploit later in our studies:

$\displaystyle\frac{\Delta T}{T}(t_{\mathrm{obs}})$
$\displaystyle\approx\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{u}}_{t_{l}}\frac{\partial
h_{\textsc{A}}\left(t_{\mathrm{ret}}\left(t,\vec{x}\right)\right)}{\partial
t}\Bigg{\rvert}_{\vec{x}=\vec{x}_{0}(t)}dt,$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{u}}_{t_{l}}\frac{dh_{\textsc{A}}\left(t^{0}_{\mathrm{ret}}\right)}{dt^{0}_{\mathrm{ret}}}dt,$
$\displaystyle=\frac{1}{2}\hat{p}^{i}\hat{p}^{j}E^{\hat{r}\textsc{A}}_{ij}\int\limits^{t_{u}}_{t_{l}}\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,$
(3.13) $\displaystyle\mathrm{where}\quad\begin{cases}\begin{tabular}[]{l l
l}$t_{u}$&= $t_{\mathrm{em}}+\frac{L}{c}$&$\approx t_{\mathrm{obs}}$ ,\\\
$t_{l}$&= $t_{\mathrm{em}}$&$\approx t_{\mathrm{obs}}-\frac{L}{c}$
.\end{tabular}\end{cases}$

Notational Aside
The following property of the sine and cosine functions will be quite useful
to us in terms of keeping our notation “clean” and convenient throughout this
paper:
$\displaystyle\begin{cases}\sin\left(2\Theta\right)&=-\cos\Big{(}2\left(\Theta+\frac{\pi}{4}\right)\Big{)},\\\
&=\quad\cos\Big{(}2\left(\Theta-\frac{\pi}{4}\right)\Big{)},\\\\[10.0pt]
\cos\left(2\Theta\right)&=\quad\sin\Big{(}2\left(\Theta+\frac{\pi}{4}\right)\Big{)},\\\
&=-\sin\Big{(}2\left(\Theta-\frac{\pi}{4}\right)\Big{)}.\end{cases}$ (3.14)

These relations allow us to write: $\frac{\partial
h_{\textsc{A}}(h,\Theta)}{\partial\Theta}=\left\\{\begin{tabular}[]{r l
r}$2h\sin\left(2\Theta\right)$&$=-2h\cos\left(2\Theta+\frac{\pi}{2}\right)$&$(+)$\\\\[5.0pt]
$-2h\cos\left(2\Theta\right)$&$=-2h\sin\left(2\Theta+\frac{\pi}{2}\right)$&$(\times)$\end{tabular}\right\\}=2h_{\textsc{A}}\left(\Theta+\frac{\pi}{4}\right).$
(3.15) This notation is motivated in the sense that it will allow us to
continually express our formulae in terms of the original $h_{\textsc{A}}$
expressions in equation 2.22, with the inclusion of an additional phase
factor. Additionally, we define the following streamlined notation:
$\begin{cases}t^{0}_{\mathrm{ret}}&\equiv
t_{\mathrm{ret}}\big{(}\vec{x}=\vec{x}_{0}(t)\big{)},\\\
\omega^{0}&\equiv\omega\left(t^{0}_{\mathrm{ret}}\right),\\\
\Theta^{0}&\equiv\Theta\left(t^{0}_{\mathrm{ret}}\right).\end{cases}$ (3.16)
In words this is the retarded time, the orbital frequency, and the orbital
phase as a function of time along the photon’s path.

Here we emphasize the limits of integration. The motion of the photon as it
travels from the pulsar to the Earth can be referenced in terms of the emitted
time or the observation time. So far we have explained this derivation with
reference to the emitted time since it is perhaps more natural. But for more
conceptual convenience going forward we will express the limits with respect
to the observation time. Again, since overall we are working on a solution
which is good only to first order in the metric perturbation, we can use the
zeroth order path of the photon to interchange the limits in this way (recall
back to equation 3.6, higher order corrections would introduce additional
factors of the metric perturbation). Additionally we make it explicitly clear
that that the metric perturbation $h_{\textsc{A}}$ is still a function of both
the orbital frequency (which appears in the amplitude, equation 2.23) and the
phase. We do this because other derivations in the literature seem to ignore
the amplitude’s dependence on the frequency or don’t properly explain the
assumptions they are making about the amplitude’s frequency dependence.

Once again, this expression is the fractional change in the period with time
due to the presence of a gravitational wave. As a reminder, conceptually we
are integrating along the photon’s path between the pulsar and the Earth. This
is achieved by integrating over the dummy time variable $t$ from the time the
photon leaves the pulsar, $t_{\mathrm{obs}}-L/c$, to the time it arrives at
the Earth and is observed, $t_{\mathrm{obs}}$. Hence the result is a function
of the observation time - i.e. it is the observed fractional change in the
period of the pulsar at some given observation time.

###  3.2 The Observable: The Pulsar Timing Residual

It may seem like the expression for $\frac{\Delta T}{T}$ in equation 3.13 is a
quantity which we could try to observe experimentally. However, as we will see
in the coming sections when we actually explicitly solve this equation, the
resulting quantity will still be extraordinarily small. For example, jumping
ahead to Section 3.5 and considering the solution of $\frac{\Delta T}{T}$ in
the plane-wave formalism in equation 3.30, aside from the geometrical terms
out front we see that the solution is effectively on the order of the metric
perturbation itself, $\mathcal{O}\left(h\right)$. If the period of our pulsar
is something like $1\mathrm{ms}$ and the magnitude of the metric perturbation
produced by our SMBHB is something like $10^{-17}$, then we would be looking
at trying to measure a fluctuation in the period of our otherwise extremely
regular pulsar clock ticks on the order of $\Delta T\sim 10^{-20}\mathrm{s}$!
So this itself is not useful to us experimentally, but it does get us one step
closer to an experimental quantity which we can measure. 222See the box at the
end of this section for a discussion of the difference between measuring the
pulsar periods vs. pulse time-of-arrivals.

Consider for a moment the integral of this quantity:

$\int\frac{\Delta T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}.$

What does this represent physically, and is it useful to us? For concreteness
let’s imagine the following toy problem. There is a clock which we are
observing, which we expect has a true period of $T=2\mathrm{s}$, but some
unknown physical effect is causing that clock period to shift by a constant
amount such that $\Delta T\left(t_{\mathrm{obs}}\right)=3\mathrm{s}$. We are
going to study this clock with a reference clock of our own for some amount of
“observation time.” Even if the observed clock doesn’t “tick” at the same rate
as ours, as clocks they should still be synchronized in that a given time
interval on the observed clock should always correspond to the same interval
of time on our reference clock.

Figure 3.1: A toy experiment. Imagine an observed clock with period
$T=2\mathrm{s}$, with some effect which causes its period to shift by $\Delta
T=3\mathrm{s}$. An interval of time lasting $6\mathrm{s}$ our time we expect
to be recorded by three periods on the observed clock if they are
synchronized. However, we observe three periods (“clock ticks”) lasting
$15\mathrm{s}$, i.e. $9\mathrm{s}$ longer than our expectation.

Imagine we begin our experiment as shown in Figure 3.1. Focus on our
observation of the time interval lasting three observed clock ticks. If we
expect the period of this observed clock to be $T=2\mathrm{s}$ our time, then
we expect three periods of the observed clock to measure out $t=6\mathrm{s}$.
However, in this example we observe this interval to last for
$t=15\mathrm{s}$. So the difference in the expected time interval and the
observed time interval of the third clock tick is $9\mathrm{s}$. This is also
what our expression tells us. If we integrate the fractional change in the
period of the observed clock over the interval of interest from our time of
$t=0\mathrm{s}$ to $t=6\mathrm{s}$, then this gives us the difference in the
expected and observed times for the three ticks:

$\int\limits^{6\mathrm{s}}_{0}\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\int\limits^{6\mathrm{s}}_{0}\frac{3\mathrm{s}}{2\mathrm{s}}dt_{\mathrm{obs}}=\frac{3}{2}(6\mathrm{s})=\quad
9\mathrm{s}\quad=15\mathrm{s}-6\mathrm{s}.$

In fact, this quantity tells us that for any interval of time that is observed
and measured:

$\int\limits\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\frac{3}{2}t,$

so if I expect to measure an interval of time $t$ on the observed clock, I
will actually observe it to be different from my expectation by
$\frac{3}{2}t$. For example,

$\displaystyle\int\limits^{2\mathrm{s}}_{0}\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\quad
3\mathrm{s}\quad=5\mathrm{s}-2\mathrm{s},$
$\displaystyle\int\limits^{4\mathrm{s}}_{0}\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\quad
6\mathrm{s}\quad=10\mathrm{s}-4\mathrm{s},$
$\displaystyle\int\limits^{4\mathrm{s}}_{1\mathrm{s}}\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\quad
4.5\mathrm{s}\quad=7.5\mathrm{s}-3\mathrm{s}.$

(In the last example here, if our clocks are in sync then the $3\mathrm{s}$
time interval between $1\mathrm{s}$ and $4\mathrm{s}$ on our own reference
clock would correspond to $1.5$ periods on our observed clock. But here we
observe $1.5$ periods to be $7.5\mathrm{s}$ long, or $4.5\mathrm{s}$ longer
than expected.)

This is a very simple toy model, but it illustrates another key concept. Say
at this point we come to some realization that the intrinsic period of the
observed clock $T$ is actually $T=5\mathrm{s}$, and that there is no unknown
physical effect causing that clock period to shift like we originally
suspected. Then if we change our model to reflect this, the observed clock’s
measurement of time will now always match our expectation (e.g. we will now
expect three periods on the observed clock to correspond to a time interval of
$15\mathrm{s}$, which is what we indeed observe). Moreover, since $\Delta
T\left(t_{\mathrm{obs}}\right)=0$ now our integral quantity will always be
zero, i.e. there will be no difference between any measured time interval on
our reference clock and on our observed clock.

We also gain more insight on this quantity by writing it in the following way,
using the definition that $\Delta T\equiv T_{\mathrm{obs}}-T$:

$\int\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\int\frac{T_{\mathrm{obs}}\left(t_{\mathrm{obs}}\right)-T}{T}dt_{\mathrm{obs}}=\int\left[\frac{T_{\mathrm{obs}}\left(t_{\mathrm{obs}}\right)}{T}-1\right]dt_{\mathrm{obs}}=\int\frac{T_{\mathrm{obs}}\left(t_{\mathrm{obs}}\right)}{T}dt_{\mathrm{obs}}-t.$

The very final term $t$ on the right hand side is just the time interval we
would expect to the observed clock to measure if it was synchronized with our
reference clock. We also know that our term on the left hand side is the
difference between the observed and expected measured time intervals on the
observed clock, so this means the quantity
$\int\frac{T_{\mathrm{obs}}}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}$
is the actual observed time interval. We can see this directly from our
earlier examples too. From the first example for the interval $t=0\mathrm{s}$
to $t=6\mathrm{s}$, the expected time is simply $t=6\mathrm{s}$, and the
observed time is
$\int\limits^{6\mathrm{s}}_{0}\frac{T_{\mathrm{obs}}}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\int\limits^{6\mathrm{s}}_{0}\frac{5\mathrm{s}}{2\mathrm{s}}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}=\frac{5}{2}(6\mathrm{s})=15\mathrm{s}$,
which is what we had earlier.

Therefore, since this quantity is a residual of the observed and expected time
intervals measured by our observed clock, we call it the “timing residual,”
which can be written as:

$\displaystyle\begin{array}[]{rll}\mathrm{Res}(t)&={\int}\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}&=\
\mathrm{Obs}(t)-\mathrm{Exp}(t),\\\
\mathrm{Res}(t,t_{0})&={\int}\limits^{t}_{t_{0}}\frac{\Delta
T}{T}\left(t_{\mathrm{obs}}\right)dt_{\mathrm{obs}}&=\
\mathrm{Obs}(t,t_{0})-\mathrm{Exp}(t,t_{0}),\end{array}$ (3.19)
$\displaystyle\mathrm{where}\quad\begin{cases}\mathrm{Exp}(t,t_{0})=t-t_{0},\\\
\mathrm{Obs}(t,t_{0})={\int}\limits^{t}_{t_{0}}\frac{T_{\mathrm{obs}}(t_{\mathrm{obs}})}{T}dt_{\mathrm{obs}}.\end{cases}$
(3.20)

Here we denote $\mathrm{Obs}(t)$ as the “observed time interval” and
$\mathrm{Exp}(t)$ as the “expected time interval” of the observed clock. If we
want to measure the residual time intervals beginning from a specific time,
then we specify $t_{0}$ in the definite integral as our initial start time.
Therefore the residual at the time $t_{0}$ will always be zero because this is
the start of the time interval and the residual measures the deviation away
from this specific time. Otherwise the indefinite integral gives us the
general expression for how the residual changes over time, without “zeroing”
it out at a specific time $t_{0}$. If we want to know the difference in the
observed time and the expected time for a specific interval, then we need to
explicitly give those limits. Another example of a $\Delta T$ which varies
sinusoidally in time, as we would expect for a time shift caused by a
gravitational wave, is given in Figure 3.2.

(a) Example of $\mathrm{Res}(t)$. In this case, the measured time intervals
spend equal amounts of time being shorter than expected as they do longer than
expected.

(b) Example of $\mathrm{Res}(t,t_{0})$, zeroed at $t=t_{0}=4$. In this case,
nearly all measured time intervals will appear to be shorter than we would
expect them to be. Figure 3.2: Here is another toy example of our experiment,
where the true period of the clock we are watching is $T=2$, and $\Delta
T(t)=A\sin\left(2\pi t/T_{\mathrm{gw}}+\theta_{0}\right)$, with the parameters
set to $A=0.1$, $T_{\mathrm{gw}}=15$, and $\theta_{0}=3\pi/8$. Both figures
are the same data set, but show how the measured time intervals can change
when we “zero” the expected and observed time intervals at a specific time
$t_{0}$. Above each timing residual is a visualization of the arrivals of the
observed and the expected clock ticks, for the times $t$ indicated on the
horizontal axis.

Physically, we can understand what is happening as follows. Although the
fluctuations of the pulsar’s period $\Delta T$ are too small to measure
directly, over time, after every pulsar period, these period shifts begin to
stack. If the values of $\Delta T>0$, after each pulsar period the additional
time added pushes the expected time-of arrival (commonly referred to as a
“TOA”) and the observed time-of-arrival of the subsequent ticks further and
further apart. If the values of $\Delta T<0$, the pulsar periods are shorter
than normal so the observed TOA of subsequent ticks comes earlier than
expected. So while an individual $\Delta T$ may not be measurable, this
measurement allows us to exploit the cumulative effect of the pulsar’s period
shift over a long period of time. If enough of these $\Delta T$’s stack on top
of each other over a long enough sequence of clock ticks, then the resulting
shift in time on the observed clock will become measurable.

This brings us to another crucial insight into the feasibility of making this
measurement possible. Ideally, we want the period of our observed clock, which
in this case is our pulsar, to be much smaller than the period of our induced
$\Delta T$ shift, which is the period of our gravitational wave, that is:

$T_{\mathrm{gw}}>T_{\mathrm{clock}}\quad\longleftrightarrow\quad
1>\omega_{\mathrm{gw}}T_{\mathrm{clock}}.$

This is because longer gravitational wave periods mean more time for the
$\Delta T$ shifts to constructively add on top of each other after each period
of the clock. This will give more time for the observed time intervals to
increase with respect to the expected time intervals, and thus will make the
timing residuals larger. If the observed clock period and the gravitational
wave period are close in magnitude to each other, then fewer clock cycles
$\Delta T$ shifts will be able to compound, and the timing residual will be
smaller. In fact we see this in our earlier assumption for a monochromatic
source, assumption 7. So this is the key physical motivation in turning this
theoretical effect (an otherwise immeasurable deviation to a pulsar’s clock
period) into one which can be measured experimentally. An example of this is
shown in Figure 3.3.

Figure 3.3: Continuation of the example shown in Figure 3.2. Here we plot two
timing residuals from two clocks with different intrinsic periods. What this
aims to show is that ideally we need to use a clock (i.e. pulsar) which has a
period $T$ much smaller than the period $T_{\mathrm{gw}}$ of the gravitational
wave inducing the $\Delta T$ shift, so that more clock cycles will allow the
$\Delta T$ shifts to constructively add and make the timing residuals larger
in magnitude, to the point where they can be measured (even if the $\Delta T$
shift itself is too small to be measurable). The smaller the $T$ (compared to
$T_{\mathrm{gw}}$), the greater the magnitude of the resulting residual.

Now we can design our experiment:

  * •

Observe a “clock tick” from a pulsar and record the time that the tick
occurred, using some master reference clock here on Earth.

  * •

Make subsequent observations like this of the pulsar’s clock ticks over the
course of years to decades. This will be our $\mathrm{Obs}(t)$ data.

  * •

For a highly regular pulsar, we will have a good idea of what the pulsar’s
clock period $T$ is from one observation. However, if a gravitational wave is
in fact adding some $\Delta T$ to it, we don’t know how many cycles may have
passed where period deviations have stacked on top of the true period.
Therefore we average a long timescale period for the pulsar from our collected
data. With this period we can then compute our $\mathrm{Exp}(t)$ data.

  * •

Given the observed time interval data and the expected time interval data we
can now compute the residual from their difference.

  * •

If there is nothing affecting the period of our pulsar then experimentally we
would expect the residual to simply look like white noise centered about
$\mathrm{Res}=0$. If, however, there is something affecting the period and our
observations are sensitive enough to detect it, then we should see some
deviation away from zero in the residual over time. We can then compare the
data to our model for a timing residual induced by a gravitational wave. For
example, this analysis could be done under a Bayesian framework in order to
try and estimate the parameters of the source from the data itself.

  * •

This entire process is iterative. As we collect more data over longer
observation timescales, we both increase our $\mathrm{Obs}(t)$ data and refine
our $\mathrm{Exp}(t)$ data, which in turn refines our $\mathrm{Res}(t)$ data
and hopefully leads to more accurate measurements of the gravitational wave
parameters.

Measuring Pulsar Periods vs. TOAs (Lam, 2020)
One potential for misconception when discussing the timing of a pulsar is how
well astronomers can currently measure pulsar periods vs. the actual TOAs of
pulsar pulses. The precision to which we can measure the periods of many
millisecond pulsars is truly astonishing - down to period uncertainties around
the order of attoseconds! However, our current ability to measure TOAs (and
hence, timing residuals) of pulsars is currently around the order of hundreds
of nanoseconds (Arzoumanian et al., 2018). It may seem strange that despite
our incredible ability to measure a pulsar’s period, our uncertainties on the
arrival times of the actual pulses from a pulsar are around nine orders of
magnitude larger! If we know a pulsar’s period $T$ down to a $\sigma_{T}$
uncertainty of order attoseconds, then given our previous discussion, why is
it perhaps not possible to measure the period shift $\Delta T$ induced by a
gravitational wave directly? Why is it we must try to measure the integrated
effect of these period shifts, i.e. the timing residuals, instead? The reason
for this discrepancy lies in how we make the two different measurements
themselves. Because millisecond pulsars are spinning so rapidly, when we point
our telescopes at them and record their light pulses, we observe many pulsar
period cycles in a relatively short amount of time. In order to measure the
pulsar’s period, all we basically need to do is divide the total observation
time by the number of pulsar pulses that occurred during that time. And this
type of observation is easily repeatable - every time we look back at that
pulsar we can repeat this

measurement, counting more pulses over more time. Notice that by doing this we
are not really trying to time the exact moment each of those pulses arrived.
All we care about is how many pulses arrived in the time we were observing the
pulsar - we simply need to be able to see that there is a pulse, and then
count that pulse. Calculating the pulsar’s period by dividing the observation
window by the number of pulses effectively averages everything out, giving us
the average period. Maybe there are gravitational waves causing those pulses
individually to arrive slightly sooner (some $T-\Delta T$) or slightly later
(some $T+\Delta T$) than we would expect, but this measurement is averaging
that out and giving us the underlying “true” period to within some (insanely
small) $\sigma_{T}$ uncertainty. But that said, if during any of those
observations we want to know the exact time every one of those pulsar light
pulses arrived, then we need to have extremely good temporal resolution on all
of the individual pulsar light pulses themselves. Given the averaging
procedure we may now know the true underlying period for that pulsar, but
being able to say just how much sooner or how much later all of the pulsar
pulses are appearing as compared to the expected theoretical TOA is more
difficult. This is why our current capabilities make this measurement many
orders of magnitude larger than measuring simply the pulsar’s underlying
period.

###  3.3 Time Retardation

For a static flat universe, the retarded time is:

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

where $\vec{x}^{{}^{\prime}}$ is the wave source’s position and $\vec{x}$ is
the field point of interest. The expansion we have made here is good when the
source is much further away than our field point. It also gives us a way of
categorizing different wave solution approximation regimes, which will
ultimately be the crux of our work to come, so gaining insight into this is
critical. Keeping successive terms after the $t$ in this expression, the first
term is our “far field” regime, the second term is the “plane-wave” regime,
and the third term is the “Fresnel” regime. One way to visually interpret the
different approximation regimes is by plotting contours of constant
$t_{\mathrm{ret}}$, as shown in Figure 3.4. The shape of the contours trace
the shape of the wavefront as it propagates away from the source. The Fresnel
term gives the first order curvature of the physical wavefront. Conceptually,
it introduces an additional time delay to the arrival time of the wavefront
predicted by the plane-wave approximation (see again Figure 3.4). Our goal is
to understand when this curvature term/time delay will become significant in
our timing residual model.

Figure 3.4: An example contour diagram of the retarded time, in the source-
Earth-pulsar plane at fixed time $t=0$. Labels on the contours give the value
of $t_{\mathrm{ret}}$ multiplied by the speed of light. Thus these values
indicated for our fixed $t=0$ the distance each contour line is from the
source. The approximations indicate up to what term is included in the
expansion in equation 3.21. Since the lines show constant $t_{\mathrm{ret}}$,
this traces the curvature of the wavefront. For illustrative purposes that
exaggerate and show the differences in the approximation regimes, we set the
pulsar at a distance $L=5$ kpc and the source at $R=30$ kpc from the Earth. A
source with frequency $\omega_{0}=1$ nHz along with this pulsar would have a
Fresnel number $F=27.3$ in this example. As we can see here, the curvature of
the wavefront means the wave arrives at the pulsar’s location later than what
is predicted by the plane-wave approximation.

As mentioned earlier, from equation 2.28 we define
$\vec{x}^{{}^{\prime}}\equiv R\hat{r}$ as the source location. In our problem
we want to follow the path of a photon along the radial line connecting the
Earth to the pulsar, i.e. along the direction $\hat{p}$. Therefore, the field
point of interest we will always be working with will be some vector that
points along this line, so we can write $\vec{x}=x\hat{p}$ (where $x\geq 0$).
(Namely the path we are interested in is given by equation 3.5 which between
our limits of integration is always positive). Hence we can write:

$\displaystyle t_{\mathrm{ret}}(t,x\hat{p})$
$\displaystyle=t-\frac{1}{c}\sqrt{R^{2}-2\left(\hat{r}\cdot\hat{p}\right)Rx+x^{2}},$
$\displaystyle=t-\frac{R}{c}+\left(\hat{r}\cdot\hat{p}\right)\frac{x}{c}-\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{x}{c}\frac{x}{R}+\ldots.$
(3.22)

In pulsar timing, typically most sources are approximated as far enough from
the Earth-pulsar baseline that gravitational waves arrive as plane-waves.
However, we are interested in seeing under what circumstances the effects of
the wavefront’s curvature from the Fresnel regime may begin to be important in
the computation of the timing residual.

Namely, the plane-wave and Fresnel regimes become increasingly different away
from the coordinate system origin, so let’s consider the maximum size of our
detector, which is the full Earth-pulsar distance $\vec{x}=L\hat{p}$:

$t_{\mathrm{ret}}(t,L\hat{p})=t\ \ -\
\underbrace{\frac{R}{c}}_{\mathrm{``Far}\ \mathrm{Field"}}\ +\ \
\underbrace{\left(\hat{r}\cdot\hat{p}\right)\frac{L}{c}}_{\mathrm{``Plane-
Wave"}}\ \ -\ \
\underbrace{\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\left(\frac{L}{c}\right)\left(\frac{L}{R}\right)}_{\mathrm{``Fresnel"}}\
\ +\ \ \ldots$ (3.23)

Now in general a wave solution behaves sinusoidally, oscillating with some
angular frequency $\omega$. For our order of magnitude estimate let’s consider
a monochromatic wave - then the wave solution evaluated at the retarded time
will behave like:

$\sin(\omega_{\mathrm{gw}}t_{\mathrm{ret}})\sim\sin\Bigg{(}\omega_{\mathrm{gw}}t-\omega_{\mathrm{gw}}\frac{R}{c}+\left(\hat{r}\cdot\hat{p}\right)\omega_{\mathrm{gw}}\frac{L}{c}-\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\omega_{\mathrm{gw}}\left(\frac{L}{c}\right)\left(\frac{L}{R}\right)\Bigg{)}.$

Since sine functions are cyclic on the interval from $0\rightarrow 2\pi$, the
Fresnel term will become significant if it is an appreciable fraction of
$2\pi$, so let’s write:

$\displaystyle 2\pi$
$\displaystyle\sim\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\omega_{\mathrm{gw}}\left(\frac{L}{c}\right)\left(\frac{L}{R}\right),$
$\displaystyle\equiv\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\pi F,$

or more roughly if we ignore all of the geometric terms and factors of order
unity in this expression we can write this condition as roughly:

$\displaystyle F$ $\displaystyle\sim 1,$ (3.24)
$\displaystyle\mathrm{where}\quad
F\equiv\frac{L^{2}}{\lambda_{\mathrm{gw}}R}=\frac{\omega_{0}}{\pi}\frac{L}{c}\frac{L}{R},$
(3.25) $\displaystyle\mathrm{and}\quad\lambda_{\mathrm{gw}}\equiv\frac{\pi
c}{\omega_{0}}.$ (3.26)

Equation 3.26 comes from the fact that the gravitational wave travels at the
speed of light, so $\lambda_{\mathrm{gw}}=cT=2\pi c/\omega_{\mathrm{gw}}$, and
the gravitational wave frequency is twice the orbital frequency of the source,
so $\omega_{\mathrm{gw}}=2\omega_{0}$. (Again, we emphasize here this is a
rough approximation to our original statement, which doesn’t worry about
factors of order unity). The quantity $F$ is the familiar “Fresnel number”
from diffraction theory, so when $F\sim\mathcal{O}(1)$, this suggests that our
plane-wave approximation may be less accurate and using the Fresnel
approximation would be required. In fact, it is because of this connection to
the Fresnel number that we have been calling this term in the retarded time
expansion the “Fresnel” term. Hence in the sections to come we will develop
the formalisms for both regimes and study how the timing residuals in these
regimes differ. To see visually the dependence of $F$ on the orbital
frequency, source distance, and pulsar distance parameters, see the contour
plot in Figure 3.9.

Of course it is still important to remember that this condition equation 3.24
should just be used to get a quick sense of which regime may be more dominant
- the geometric terms which we ignored can become important, but we will
discuss those in detail in the coming sections.

Natural Plane-Wave Limit
Note there is an important limit here wherein higher order terms past the
plane-wave term in equation 3.23 naturally vanish. Writing out additional
terms in the Taylor expansion in equation 3.22, we would find that in addition
to the “Fresnel” term, all higher terms contain a factor of
$\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)$. This means when
$\hat{r}\cdot\hat{p}=\pm 1$, all terms except for the plane-wave term vanish -
hence we refer to this as the “natural plane-wave limit.” Physically, this
corresponds to the pulsar and the source being either “aligned” ($+1$) or
“anti-aligned” ($-1$) with the Earth. We can see why all effects of curvature
of the wavefront would not be noticeable in this limit simply from considering
Figure 3.4. All points along the vertical line that runs through the Earth
marker (pointing along the arrow) “see” the incoming wavefront as a plane -
there is no curvature along the wavefront in this precise geometrical
alignment. This can also be seen in Figure 3.6, where no induced change in the
pulsar’s period would occur in either of these alignments.

Finally in our problem, in order to solve for the timing residual we will need
to evaluate the retarded time along the path of the photon as it travels
between the pulsar and the Earth, given by equation 3.5. Inserting this into
equation 3.22 gives (using our notation in equation 3.16):

$\displaystyle t_{\mathrm{ret}}^{0}$
$\displaystyle=t-\frac{R}{c}\sqrt{1-2\hat{r}\cdot\hat{p}\left(\frac{ct_{\mathrm{obs}}-ct}{R}\right)+\left(\frac{ct_{\mathrm{obs}}-ct}{R}\right)^{2}},$
$\displaystyle=t\ -\frac{R}{c}\
+\left(\hat{r}\cdot\hat{p}\right)\frac{ct_{\mathrm{obs}}-ct}{c}\
-\frac{1}{2}\left(1-\left(\hat{r}\cdot\hat{p}\right)^{2}\right)\frac{\left(ct_{\mathrm{obs}}-ct\right)^{2}}{cR}\
+\ \ldots.$ (3.27)

###  3.4 Four Pulsar Timing Regimes

To date, pulsar timing literature has been primarily considered timing models
derived assuming that the PTA receives plane-waves coming from distant
sources. They have also been derived assuming that the frequency of the
gravitational wave is monochromatic, or that in some cases the frequency
evolves slightly in the thousands to tens of thousands of years it takes light
to travel from the pulsar to the Earth. We can classify these as two separate
model regimes, “plane-wave monochromatic” and “plane-wave frequency
evolution.”

One of the primary goals of my work has been to further generalize these
models by removing the plane-wave assumption. In my work I add two new and
complimentary timing regimes to the literature, represented in Figure 3.5.
These regimes increase in generality from left-to-right and top-to-bottom. The
new regime I have labeled the “Fresnel” regime as we will show it becomes
important for significant Fresnel numbers. To this I have derived an analytic
formula for the timing residuals in the Fresnel monochromatic regime, and I
propose a physically motivated conjecture as to what the most general timing
model should be, the Fresnel frequency evolution regime. The Fresnel frequency
evolution regime recovers all of the previously predicted results by each of
the other three regimes in the appropriate limits. In general, the frequency
evolution regime reduces to the monochromatic regime in the large coalescence
time limit $\Delta\tau_{c}\rightarrow\infty$ (see Section 2.7). And the
Fresnel regime reduces back to the plane-wave regime in either the natural
plane-wave limit $\hat{r}\cdot\hat{p}=\pm 1$ or the small Fresnel number limit
$F\rightarrow 0$.

Figure 3.5: My classification of the gravitational wave regimes coming from
sources of interest. They increase in generality from left-to-right and top-
to-bottom. The models reduce from frequency evolution to monochromatism in the
large coalescence time limit ($\Delta\tau_{c}\rightarrow\infty$), and from
Fresnel to plane-wave in either the small Fresnel number limit ($F\rightarrow
0$) or the natural plane-wave limit ($\hat{r}\cdot\hat{p}\rightarrow\pm 1$).
The plane-wave regime has been well studied and used in pulsar timing