Realistic sensitivity curves for pulsar timing arrays

TL;DR

This paper develops realistic sensitivity curves for pulsar timing arrays that account for noise and model fitting, enabling quicker assessment of gravitational wave detectability.

Contribution

It introduces a method to construct approximate sensitivity curves incorporating noise and fitting effects, validated against exact calculations.

Findings

Sensitivity curves agree within tens of percent with exact methods.

The approach significantly reduces computation time.

Useful for rapid evaluation of gravitational wave signals.

Abstract

We construct realistic sensitivity curves for pulsar timing array searches for gravitational waves, incorporating both red and white noise contributions to individual pulsar noise spectra, as well as the effect of fitting to a pulsar timing model. We demonstrate the method on both simulated pulsars and a realistic array consisting of a subset of NANOGrav pulsars used in recent analyses. A comparison between the results presented here and exact sensitivity curves shows agreement to tens of percent. The resulting sensitivity curves can be used to assess the detectability of predicted gravitational-wave signals in the nanohertz frequency band in a fraction of the time that it would take to compute the exact sensitivity curves.

Figures

Equations

$$t=t^{\det }(\xi )+n+h.$$

$$\delta t\equiv t-t^{\det }({\xi }_0)=M\delta \xi +n+h,$$

$$M\equiv\left(\frac{\partial t^{\rm det}}{\partial\xi}\right)\bigg{|}_{\xi=\xi_{0}}$$

$$p(\delta t|\delta \xi ,C_n,C_h,\theta )\propto$$

$$\exp \left[-\frac{1}{2}(\delta t-M\delta \xi -h(\theta ){)}^T{C}^{-1}(\delta t-M\delta \xi -h(\theta ))\right]$$

$$C\equiv C_n+C_h$$

$$\delta {\xi }_{\mathrm{ML}}=(M^T{C}^{-1}M{)}^{-1}{M}^T{C}^{-1}\delta t.$$

$$\delta t^{\mathrm{post}}\equiv \delta t-M\delta {\xi }_{\mathrm{ML}}=R\delta t,$$

$$R\equiv 1-M(M^T{C}^{-1}M{)}^{-1}{M}^T{C}^{-1}.$$

$$p(\delta t|C_n,C_h,\theta )\propto$$

$$\exp \left[-\frac{1}{2}(\delta t-h(\theta ){)}^{T}G(G^{T}CG{)}^{-1}{G}^T(\delta t-h(\theta ))\right],$$

$$M=US{V}^T,U=(F,G).$$

$$r\equiv G^T\delta t,$$

$${\sigma }_r^2={\int }_0^{\infty }\mathrm{d}f\mathcal{T}(f)P(f),$$

$$\mathcal{T}(f)\equiv \frac{1}{N}\sum _{k,l}(G{G}^T{)}_{kl}{e}^{i2\pi f(t_k-t_l)}.$$

$${\sigma }_{\mathrm{post}}^2={\int }_0^{\infty }\mathrm{d}f{\mathcal{T}}_R(f)P(f),$$

$${\mathcal{T}}_R(f)\equiv \frac{1}{N}\sum _{k,l}{R}_{kl}{e}^{i2\pi f(t_k-t_l)}.$$

$${\chi }^2\equiv (\delta t-h(\theta ){)}^{T}G(G^{T}CG{)}^{-1}{G}^T(\delta t-h(\theta )).$$

$$h_k(\theta )\equiv h(t_k;\theta )={\int }_{-f_{\mathrm{Nyq}}}^{f_{\mathrm{Nyq}}}\mathrm{d}f\tilde{h}(f;\theta )e^{i2\pi f{t}_k},$$

$$\begin{gathered} {\chi }^2=2T{\int }_{-f_{\mathrm{Nyq}}}^{f_{\mathrm{Nyq}}}\mathrm{d}f{\int }_{-f_{\mathrm{Nyq}}}^{f_{\mathrm{Nyq}}}\mathrm{d}{f}^{\prime }(\widetilde{\delta t}(f)-\tilde{h}(f;\theta )) \\ \times {\mathcal{N}}^{-1}(f,f^{\prime })(\widetilde{\delta t}{}^{\ast }(f^{\prime })-{\tilde{h}}^{\ast }(f^{\prime };\theta )), \end{gathered}$$

$$\begin{gathered} {\mathcal{N}}^{-1}(f,f^{\prime }) \\ \equiv \frac{1}{2T}\sum _{k,l}{e}^{i2\pi f{t}_k}[G(G^{T}CG{)}^{-1}{G}^T{]}_{kl}{e}^{-i2\pi f^{\prime }{t}_l}. \end{gathered}$$

$${\mathcal{N}}^{-1}(f)\equiv \frac{1}{2T}\sum _{k,l}[G(G^{T}CG{)}^{-1}{G}^T{]}_{kl}{e}^{i2\pi f(t_k-t_l)},$$

$${\mathcal{N}}^{-1}(f)=\mathcal{T}(f)/P(f).$$

$$C_{ij}={\int }_0^{f_{\mathrm{Nyq}}}\mathrm{d}f\cos [2\pi f(t_i-t_j)]P(f),$$

$$P(f)=2{\sigma }^2\Delta t+A{f}^{-\gamma },\gamma >0,$$

$${\mathcal{N}}^{-1}(f)\approx \mathcal{T}(f)/P(f).$$

$${\mathcal{D}}_N((f-f^{\prime })\Delta t)\equiv \frac{1}{N}\frac{\sin [N\pi (f-f^{\prime })\Delta t]}{\sin [\pi (f-f^{\prime })\Delta t]}.$$

$${\mathcal{N}}^{-1}(f,f^{\prime })$$

$$\simeq P^{-1}(f){\delta }_{f{f}^{\prime }}.$$

$$\hat{k}$$

$$\hat{l}$$

$$\hat{m}$$

$$\cos \iota \equiv \hat{k}\cdot \hat{L},\cos \psi \equiv \hat{u}\cdot \hat{l},$$

$$\hat{u}\equiv \frac{\hat{L}\times \hat{k}}{|\hat{L}\times \hat{k}|},\hat{v}\equiv \hat{k}\times \hat{u},$$

Full text

I Motivation

Pulsar timing arrays (PTAs) are poised to make the first detection of nanohertz gravitational waves (GWs) in the next 2-5 years Siemens et al. (2013); Rosado et al. (2015); Taylor et al. (2016); Kelley et al. (2017). These galactic-scale GW detectors use the correlated times of arrival (TOAs) from millisecond pulsars to search for GWs Sazhin (1978); Detweiler (1979); Foster and Backer (1990). The recent inception of observational relativity by the advanced LIGO and VIRGO ground-based detectors Abbott et al. (2016, 2018) and the multi-messenger observations of binary neutron stars Abbott et al. (2017) have drastically changed our understanding of stellar-mass compact objects. PTAs are poised to complement these observations by observing GWs from binary systems comprised of super-massive black holes (SMBHs) in the centers of distant galaxies.

A common tool used to assess the observability of GW sources across the spectrum are detection sensitivity curves (see, e.g., Moore et al. (2015a); Robson et al. (2019) and Figure 1). These curves are basic “figures of merit,” constructed by the developers of GW observatories to assess the sensitivity of current detectors and to predict the sensitivity of future, next-generation detectors. The wider astrophysics community uses detection sensitivity curves as an initial estimate of the ability of a given detector to observe GWs from a particular source.

While detailed sensitivity curves for extant detectors are usually published for each observation run, those for PTAs are often simplified Hobbs (2011); Thrane and Romano (2013), only including identical white-noise components and often assuming that all pulsar observation epochs are evenly spaced and have the same baseline of observations. When drawn, these curves are often cut-off at the timespan of the observations and do not include important insensitivities at frequencies of $1/{\rm yr}$ and $2/{\rm yr}$, due to fitting for a pulsar’s astrometric parameters (Figure 1).

It has long been known that the fit to a pulsar’s timing model acts as a filter function Cordes (1980); Blandford et al. (1984), absorbing frequencies in the pulsar timing data in a predictable manner. These effects have been studied in the context of searches for GWs Cutler et al. (2014); Taylor et al. (2016); Blandford et al. (1984). Reference Madison et al. (2013) go one step further, showing how one can mitigate for losses in sensitivity using very-long-baseline interferometry to localize pulsars sky locations without explicitly fitting for their positions using the timing data.

Modern PTA data analysis strategies and algorithms are designed with this complication of the timing model fit in mind Demorest (2007); van Haasteren et al. (2009); Coles et al. (2011); van Haasteren and Levin (2013); Chamberlin et al. (2015); Lee et al. (2012); Lee (2016). This formalism was used e.g., in Moore et al. (2015b) to study PTA sensitivity curves for deterministic and stochastic sources of GWs, calculating sensitivity curves both analytically and numerically, using frequentist and Bayesian methods. The approach in Moore et al. (2015b) is similar in spirit to ours in that they start from the same likelihood function as we do (Section II.2), and they use properties of the expected signal-to-noise ratios for deterministic and stochastic GW signals to start to incorporate the effect of timing model fits. Our analysis differs from theirs in that we explicitly identify a component of the likelihood function that encodes both the noise power spectral density in a given pulsar’s data set and the effects of the timing model fit. This information is combined with known sources of realistic noise in pulsar timing data, including time-correlated (red) noise, to construct sensitivity curves for individual pulsars. (Reference Caballero et al. (2016) also discusses the effect of red noise on the sensitivity of pulsar timing searches for GWs, using a Fisher matrix calculation to estimate the errors.) For an array of pulsars, we use the expected signal-to-noise ratio of detection statistics for both deterministic and stochastic GW signals to construct effective sensitivity curves for the whole array.

I.1 Plan of paper

In Section II, we describe the basic formalism underlying pulsar timing analyses—i.e., timing residuals, timing models, and the effect of fitting to a timing model. This leads us to timing-model-marginalized residuals and their associated transmission functions, which play a key role in the subsequent construction of detection sensitivity curves. In Section III, we describe in detail the response of pulsar timing measurements to both deterministic and stochastic GWs. Then, in Section IV, we introduce detection statistics for both types of signals. The expressions for their corresponding expected signal-to-noise ratios allow us to read off an effective strain-noise power spectral density for the PTA, which has the interpretation of a detection sensitivity curve. As an application of our analysis, we construct sensitivity curves for the NANOGrav 11-yr pulsars using realistic noise properties and timing model fits, and compare our predicted sensitivities to published upper limits. We conclude in Section V. We also include Appendix A, in which we cast the results of an early seminal paper Blandford et al. (1984) into the more modern notation used in recent pulsar timing analyses.

The calculations provided in this work are packaged in a Python package available on the Python Package Inventory (PyPI) and GitHub.

II Pulsar timing analyses

Here we review the formalism underlying pulsar timing analyses used in GW searches. Readers interested in more details should see Blandford et al. (1984); Demorest (2007); Ellis et al. (2012); van Haasteren and Levin (2013); Chamberlin et al. (2015).

II.1 Times of arrival and timing residuals

Let us start with a single pulsar. The measured pulse times of arrival (TOAs) consist of three parts:111To simplify the notation, we have not included indices to label the particular pulsar ($I=1,2,\cdots,N_{\rm p}$), the individual TOAs ($i=1,2,\cdots,N$), or the timing model parameters ($a=1,2,\cdots,N_{\rm par}$). If one wants to include those indices explicitly, one should write $t_{Ii}=t_{Ii}^{\rm det}(\xi_{a})+n_{Ii}+h_{Ii}$.

[TABLE]

The first term gives the expected TOAs due to deterministic processes, which depend on intrinsic properties of the pulsar (e.g., its spin period, period derivative, …), extrinsic properties of the pulsar (e.g., its sky location, proper motion, distance from the solar system barycenter, …), and processes affecting the pulse propagation (e.g., disperion delays due to the interstellar medium, relativistic corrections, …). The timing model parameters are denoted by $\xi$. The second term is (stochastic) noise intrinsic to the pulsar or to the measurement process itself. The third term is a perturbation to the pulse arrival times induced by GWs, which in general will have contributions from both deterministic and stochastic sources, $h=h^{\rm det}+h^{\rm stoch}$.

Timing residuals are then defined by subtracting the expected TOAs (predicted by the timing model for an initial estimate of the model parameters $\xi_{0}$) from the measured TOAs:

[TABLE]

where

[TABLE]

is the design matrix. The above expression for $\delta t$ is obtained by Taylor expanding the timing model $t^{\rm det}(\xi)$ around the initial parameter estimates $\xi_{0}$, assuming that the initial estimates are close enough to the true values that only 1st-order terms in the parameter deviations $\delta\xi$ are needed in the expansion. The design matrix $M$ is a rectangular matrix of dimension $N\times N_{\rm par}$, with components $M_{ia}$. Each column of the design matrix encodes the linearized fit to one parameter in the timing model.

II.2 Fitting to a timing model

From the form of (2), one sees that errors $\delta\xi$ in our orignal estimate $\xi_{0}$ of the timing model parameters lead to deterministic features in the timing residuals. For example, an error in the pulse period leads to timing residuals that grow linearly with time, $\delta t\sim t$, while an error in the period derivative leads to residuals that grow quadratically with time, $\delta t\sim t^{2}$. Thus, we can improve our estimates of the timing model parameters by fitting for $\delta\xi$ in our linear timing model for the residuals.

This can be done in two ways, both of which take the likelihood function

[TABLE]

as the starting point. In the above expression,

[TABLE]

is the noise covariance matrix, which has contributions from both detector noise $C_{n}$ (i.e., noise intrinsic to the pulsar and from the measurement process) and a potential GW background $C_{h}$. The term $h(\theta)$ are the timing residuals induced by a deterministic GW source (e.g., the expected waveform from an individual SMBH binary parametrized by $\theta$).

(i) The first approach to fitting to the timing model is to maximize the likelihood function with respect to the parameter deviations $\delta\xi$. Since $\delta\xi$ appears linearly in the expression for the timing residuals (quadratically in the argument of the exponential), the maximization is easy to do. One obtains the standard result

[TABLE]

From these maximum-likelihood estimates, we can then form post-fit residuals

[TABLE]

Note that $R$ is an $N\times N$ matrix that implements the fit to the linear timing model; it depends in general on both the timing model (via $M$) and the detector noise (via $C$). One can show that $R$ is a projection operator ($R^{2}=R$), and hence not invertible.

(ii) The second approach to fitting to the timing model is to marginalize the likelihood function over the parameter deviations $\delta\xi$, assuming flat priors for $\delta\xi$. The result of this marginalization is the timing-model-marginalized (TMM) likelihood function van Haasteren et al. (2009); van Haasteren and Levin (2013)

[TABLE]

where $G$ is an $N\times(N-N_{\rm par})$ matrix constructed from a singular-value decomposition of the design matrix

[TABLE]

Note that $G$ depends only on the timing model (via $M$) and not on the noise. In terms of components, $G\equiv G_{i\alpha}$, where $\alpha=1,2,\cdots,N-N_{\rm par}$. Using $G$, one can construct associated TMM residuals

[TABLE]

which are orthogonal to the timing model. Since $U$ is a unitary matix, it follows that $[G^{T}G]_{\alpha\beta}=\delta_{\alpha\beta}$. For white noise (i.e., $C$ proportional to the identity matrix), we have the identitiy $R=GG^{T}$.

Although both approaches for fitting to the timing model have been used in the past, in this paper we will use the second approach, given that it is the one used most often for current pulsar timing array searches for GWs.

II.3 Transmission functions

The process of fitting to a timing model removes power from the post-fit or TMM residuals. This can be easily demonstrated by calculating the variance of the TMM residuals $r\equiv G^{T}\,\delta t$. One finds

[TABLE]

where $P(f)$ is the (one-sided) power spectral density of the original (pre-fit) timing residuals $\delta t$, and

[TABLE]

The function ${\cal T}(f)$ has the interpretation of a transmission function, selectively removing power associated with the timing model fit. A plot of ${\cal T}(f)$ for a simple timing model consisting of quadratic spin-down (i.e., fitting to the phase offset, spin period, and period derivative of the pulsar), the pulsar’s sky position, and the distance to the pulsar is shown in Figure 2.

Note that fitting to the sky position absorbs power at and around a frequency of 1/year, corresponding to the Earth’s yearly orbital motion around the Sun. Fitting to the pulsar distance absorbs power at a frequency of 2/year, which corresponds to a parallax measurement. The quadratic spin-down parameter fit acts as a high-pass filter, absorbing frequencies substantially below $1/T$, where $T$ is the time span of the data. The effect of the observing time on the shape of the transmission function is shown in Figure 2.

Pulsars in binaries famously have additional components to the timing model that take into account the various Doppler shifts due to binary motion and relativistic effects, if the line-of-sight passes by the companion (Shapiro delay) or if the binary is in a tight enough orbit to observe the loss of power due to GWs Taylor and Weisberg (1982). These components of the timing model have a minimal effect on sensitivity curves for GWs as the frequencies in question are much higher than those of the sources for which PTAs are searching. We do not include these components when simulating pulsar design matrices, but we will see the (mostly subtle) changes they make when looking at the design matrices of real pulsar data.

Finally, we note that one can also calculate an analogous transmission function associated with the post-fit timing residuals $\delta t^{\rm post}\equiv R\,\delta t$. One finds

[TABLE]

where

[TABLE]

This $R$-matrix transmission function was originally described in Blandford et al. (1984), although from a slightly different perspective. In Appendix A, we cast the approach of Blandford et al. (1984) into the more modern $R$-matrix notation.

II.4 Inverse-noise-weighted transmission function

It turns out that there is another way of obtaining a quantity that behaves like a transmission function by working directly with the TMM likelihood (9). The argument of the exponential can be written as $-\frac{1}{2}\chi^{2}$, where

[TABLE]

If we write this in the Fourier domain by substituting

[TABLE]

where $t_{k}\equiv k\Delta t$ and $f_{\rm Nyq}\equiv 1/(2\Delta t)$, we find

[TABLE]

where

[TABLE]

The quantity ${\cal N}^{-1}(f,f^{\prime})$ is a function of two frequencies, $(f,f^{\prime})$, but it turns out to be diagonally-dominated, with the majority of its support on the diagonal $f=f^{\prime}$, as shown in Figure 3. (The broadening of the diagonal band at low frequencies is an artefact of using log-scale axes for the frequencies.) The diagonal component

[TABLE]

and three off-diagonal cross-sections of ${\cal N}^{-1}(f,f^{\prime})$ are shown in Figure 4. (The fact that the off-diagonal cross-sections are curved in panel (a) of Figure 4 is again due to using log-scale axes for the frequencies.)

A few remarks are in order:

(i) For this particular example, the diagonal component ${\cal N}^{-1}(f)$ is identical in shape with the transmission function ${\cal T}(f)$ shown in Figure 2. The amplitude of ${\cal N}^{-1}(f)$ differs from ${\cal T}(f)$ by a constant factor $1/P(f)=1/(2\sigma^{2}\,\Delta t)$, corresponding to a white noise covariance matrix.222For our white noise simulations, we take $P(f)=2\sigma^{2}\,\Delta t$, with $\sigma=100\leavevmode\nobreak\ {\rm ns}$ and $\Delta t={\rm yr}/20$. These numerical values are often chosen for pulsar timing simulations. Thus, for white noise

[TABLE]

This is illustrated in Figure 5(a). If we also include red noise in the noise covariance matrix $C$ by taking

[TABLE]

then the relationship between ${\cal N}^{-1}(f)$ and ${\cal T}(f)/P(f)$ is only approximate,

[TABLE]

This is illustrated in Figure 5(b).

(ii) Away from the dip at 1/yr, where there is suppression of power due to the timing model fit to the pulsar sky position, the off-diagonal cross-sections are proportional to Dirichlet sinc functions

[TABLE]

When multiplied by $T$, a Dirichlet sinc function can be thought of as finite-time approximation to the Dirac delta function—i.e., $\delta(f-f^{\prime})\simeq T{\cal D}_{N}((f-f)^{\prime}\Delta t)$. Dirichlet sinc functions arise when taking the Fourier transform of a discretely-sampled rectangular window of duration $T=N\Delta t$, see e.g., Romano and Cornish (2017). This diagonally-dominated behavior is what you would expect for ${\cal N}^{-1}(f,f^{\prime})$ if one had only Gaussian-stationary noise. This is the case if one doesn’t have to fit a timing model (Figure 3). Then one can simply replace $G$ by the identity matrix, for which

[TABLE]

The approximate equality in the above equation is a consequence of the Karhunen-Loeve theorem, which states that the discrete Fourier transform operation defined by the unitary matrix $U_{jk}\equiv\frac{1}{N}e^{-i2\pi jk/N}$ approximately diagonalizes a stationary covariance matrix in the limit that the observation time $T=N\,\Delta t$ is much larger than the correlation time of the noise.

(iii) Since fitting to a timing model introduces non-stationarities into the TMM residuals van Haasteren and Levin (2013), one cannot directly appeal to the Karhunen-Loeve theorem for the general expression (19). One needs to explicitly check the validity of the diagonal approximation for ${\cal N}^{-1}(f,f^{\prime})$ as we have done in Figures 3 and 4. We have also numerically computed the sum of ${\cal N}^{-1}(f,f^{\prime})$ over the full two-dimensional array of frequencies $(f,f^{\prime})$ and compared that to the sum of ${\cal N}^{-1}(f,f^{\prime})$ just along the diagonal $f=f^{\prime}$. Even for the more challenging case of a red+white noise covariance matrix (Figure 5(b)) and a fit to the our quadratic spin-down model, the two summations agree to within $\approx 6\%$.

III Timing residual response to gravitational waves

To proceed further in our calculation of pulsar timing sensitivity curves, we need to describe in more detail the timing residual response of a pulsar to an incident GW. We will consider both deterministic and stochastic sources of GWs. Interested readers should see Blandford et al. (1984); Demorest (2007); van Haasteren et al. (2009); Coles et al. (2011); van Haasteren and Levin (2013) for more details.

III.1 Response to a single deterministic source

We will start by writing down the metric perturbations $h_{ab}(t,\vec{x})$ for a single deterministic source emitting plane GWs in the direction $\hat{k}$ (Figure 6). To do this we introduce two coordinate frames: one associated with the solar system barycenter (SSB) and the other associated with the propagation of the GW. We will assume that the source has a symmetry axis (e.g., the direction of the orbital angular momentum vector $\vec{L}$ for a binary system), and that the symmetry axis makes an angle $\iota$ with respect to the line of sight $\hat{k}$ from the GW source to the solar system barycenter, and an angle $\psi$ with respect to the vector $\hat{l}$ when projected onto the plane perpendicular to $\hat{k}$ (Figure 7).

The vectors $\hat{k}$, $\hat{l}$, $\hat{m}$ are defined in the solar system barycenter frame by

[TABLE]

where $(\theta,\phi)$ are the standard polar and azimuthal angles on the 2-sphere in equatorial coordinates, and the origin of coordinates is at the solar system barycenter. The right ascension $\alpha$ and declination $\delta$ of a source are given in terms of $\theta$ and $\phi$ by $\alpha=\phi$ and $\delta=\pi/2-\theta$.

The angles $\iota$ and $\psi$ are the inclination and polarization angles of the source, respectively. They can be written in terms of the unit vectors $\hat{k}$, $\hat{l}$, $\hat{L}\equiv\vec{L}/|\vec{L}|$, and $\hat{u}$ via:

[TABLE]

where

[TABLE]

are two orthogonal unit vectors in the plane perpendicular to $\hat{k}$ (Figure 7). Note that $\iota=0$ or $\pi$ corresponds to the orbital plane being seen face-on or face-off; $\iota=\pi/2$ or $3\pi/2$ corresponds to seing the orbital plane edge on. The unit vectors $\hat{u}$, $\hat{v}$ are related to $\hat{l}$, $\hat{m}$ by a rotation around $\hat{k}$ through the angle $\psi$ as shown in Figure 7.

From $\hat{u}$ and $\hat{v}$, we can construct a preferred set of polarization tensors:

[TABLE]

Using these polarization tensors, we can expand the metric perturbations:

[TABLE]

or, equivalently,

[TABLE]

where $\tilde{h}_{+,\times}(f;\iota)$ are the Fourier transforms of $h_{+,\times}(t;\iota)$. The timing residual response of a pulsar to such a deterministic GW is then Detweiler (1979):

[TABLE]

where

[TABLE]

with

[TABLE]

Here $\hat{p}^{a}$ is a unit vector pointing from the solar system barycenter to the pulsar, and $D$ is the distance to the pulsar. The function $R^{+,\times}(f,\hat{k},\psi)$ is the timing residual response function of a pulsar to a monochromatic plane GW propagating in direction $\hat{k}$, with frequency $f$, polarization $+,\times$, and polarization angle $\psi$. The two terms in the response function are called the ‘Earth term’ and ‘pulsar term’, respectively, since they involve sampling the GW phase at Earth and at the location of the pulsar, a distance $D$ away from Earth. The factor of $1/(i2\pi f)$ comes from the fact that we are working with timing residuals, as opposed to Doppler shifts in the pulse frequency.

For the analyses that we will do in this paper, we will typically ignore the pulsar-term contribution to the timing residual response to GWs, as this term will not contribute to the cross-power when correlating the signal associated with distinct pulsars. (The separation between pulsars ($\sim{\rm kpc}$) is much greater the wavelengths of the GWs that we are sensitive to, which are of order $\lesssim 10\leavevmode\nobreak\ {\rm lyr}$.) There is a contribution, however, to the auto-correlated power for a single pulsar, which comes from the exponential part of $|R^{+,\times}(f,\hat{k},\psi)|^{2}$:

[TABLE]

where we have ignored the cosine term since it is a rapidly-oscillating function of the GW propagation direction $\hat{k}$, and hence does not contribute significantly when summed over the sky. The value ‘2’ corresponds to the sum of the Earth-Earth and pulsar-pulsar auto-correlation terms.

III.1.1 Circular binaries

To proceed further, we need to specify the form of $h_{+,\times}(t;\iota)$ or its Fourier transform $\tilde{h}_{+,\times}(f;\iota)$. For example, for a circular binary

[TABLE]

where $\Phi(t)$ is the orbital phase and $h_{0}(t)$ is a dimensionless amplitude given by

[TABLE]

Here $D_{L}$ is the luminosity distance to the source, $\mathcal{M}_{\rm c}\equiv(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ is the chirp mass of the binary system, and $\omega(t)$ is the instantaneous orbital angular frequency, $\Phi(t)=\int^{t}{\rm d}t^{\prime}\>\omega(t^{\prime})$. For an evolving binary system

[TABLE]

which is a consequence of energy balance between the radiated power in GWs and the orbital energy lost by the binary system. The instantaneous GW frequency $f(t)$ is related to the orbital frequency $\omega(t)$ via $\omega(t)=\pi f(t)$.

The above differential equation for $\omega(t)$ (or, equivalently, for $f(t)$) can be integrated to yield

[TABLE]

where $\tau\equiv t_{\rm col}-t$ is the time to coalescence. Inverting (40), we obtain

[TABLE]

which is the time to coalescence for a binary system currently having GW frequency $f$. Note that for a SMBH binary with $10^{9}$ solar-mass BHs (which is the primary source for PTAs) and GW frequency $f=8\leavevmode\nobreak\ {\rm nHz}$ (which is one of the most sensitive frequencies for the current decade-long PTA searches), the time to coalescence is $\tau\sim 10^{5}\leavevmode\nobreak\ {\rm yr}$, which is four orders of magnitude larger than a decade-long observation $T=10\leavevmode\nobreak\ {\rm yr}$. Over the course of the observation the change in the GW frequency for the above SMBH binary is

[TABLE]

which is four orders of magnitude smaller than the frequency bin width $1/T$, set by the total observation time $T$. Thus, for the purposes of this paper, we will take our deterministic source to be a monochromatic binary with $f(\tau)=f_{0}={\rm const}$.

With this simplification, equations (37) and (38) become

[TABLE]

where $\phi_{0}$ is the initial phase and $h_{0}$ is the (constant) strain amplitude

[TABLE]

The Fourier transforms of $h_{+,\times}(t;\iota,\phi_{0})$ are then

[TABLE]

But since the signals are observed for only a finite duration, the Dirac delta functions $\delta(f\mp f_{0})$ should be replaced by their finite-time equivalents $\delta_{T}(f\mp f_{0})$ defined by

[TABLE]

where $T$ is the observation time for the pulsar. If one wants to also include the discreteness $\Delta t$ of the time-series data, then the Dirac delta functions should be replaced by Dirichlet sinc functions, $T\mathcal{D}_{N}[(f\mp f_{0})\Delta t]$ (see (25)). It turns out that the final (approximate) expressions that we obtain, cf. (48) and (52), are independent of which finite-time approximation we use.

III.1.2 Averaging over inclination, polarization, and sky position

Using the above expressions for $\tilde{h}_{+,\times}(f;\iota,\phi_{0})$ and (35) for $R^{+,\times}(f,\hat{k},\psi)$, we can calculate the squared response $|\tilde{h}(f)|^{2}$ averaged over the inclination of the source (defined by the inclination and polarization angles $\iota$ and $\psi$), initial phase $\phi_{0}$, and sky direction $\hat{n}\equiv-\hat{k}$. This is relevant for the case where these quantities are not known a priori. Defining

[TABLE]

it is fairly easy to show that

[TABLE]

where

[TABLE]

The factor of $4/5$ in (48) comes from the average over inclination angles $(\iota,\psi)$; $\mathcal{R}(f,\hat{k})$ encodes the timing residual response of a pulsar to a plane GW propagating in direction $\hat{k}$ averaged over the $(+,\times)$ polarizations and the polarization angle $\psi$; and $S_{h}(f)$ is the strain power-spectral density of a monochromatic GW having frequency $f_{0}$. The approximate equality in (48) is there because we made the approximation $\delta_{T}^{2}(f\mp f_{0})\simeq T\delta(f\mp f_{0})$ for the product of two finite-time Dirac delta functions. This allows us to write $S_{h}(f)$ in terms of ordinary Dirac delta functions, which are formally singular at $f=\pm f_{0}$. But this is not a problem, as $S_{h}(f)$ will only need to be evaluated under an integral sign for the expected signal-to-noise ratio calculations that we will perform in Section IV.1. This approximation gives answers that are good to within $\lesssim 10\%$ for noise power spectral densities that don’t vary significantly over a frequency bandwidth $\Delta f\sim 1/T$ in the neighboorhood of $\pm f_{0}$.

If we also average over sky location, defining

[TABLE]

we find

[TABLE]

where

[TABLE]

Note that the expression for ${\cal R}(f)$ is independent of the direction $\hat{p}$ to the pulsar. The above expressions will be used later on when defining the detection sensitivity curves in Section IV.

III.2 Response to a stochastic GW background

For a stochastic GW background, the metric perturbations can be written as a superposition of plane GWs having different frequencies $f$, polarizations $\{+,\times\}$, and propagation directions $\hat{k}$:

[TABLE]

where $e^{+,\times}_{ab}(\hat{k})\equiv\epsilon^{+,\times}_{ab}(\hat{k},0)$. This is basically (32) but allowing for contributions from different propagation direction $\hat{k}$. Since we will assume that the sources producing the GW background have no preferred polarization direction or symmetry axis, we have set $\psi=0$ and $\iota=0$ in the expansion for $h_{ab}(t,\vec{x})$. The timing residual response of a pulsar to the background is then

[TABLE]

where

[TABLE]

with $R^{+,\times}(f,\hat{k},0)$ given by (35). As discussed there, we will generally ignore the contribution of the pulsar term to the response function, except when calculating the auto-correlated power, which will have contributions from both the Earth-Earth and pulsar-pulsar auto-correlation terms.

The Fourier components $\tilde{h}_{+,\times}(f;\hat{k})$ that enter the plane-wave expansion of the metric perturbations are random fields. Their quadratic expectation values completely define the statistical properties of the background, under the assumption that it is Gaussian-distributed. For simplicity, we will assume that the GW background is stationary, unpolarized, and isotropic,333See e.g., Romano and Cornish (2017) for a review of analyses that drop these assumptions. for which $\langle\tilde{h}_{P}(f;\hat{k})\rangle=0$ and

[TABLE]

where $P=\{+,\times\}$. Here $S_{h}(f)$ is the (one-sided) strain power spectral density of the background (units of ${\rm strain}{}^{2}/{\rm Hz}$), which is related to the dimensionless energy-density spectrum $\Omega_{\rm gw}(f)$ via

[TABLE]

It is also common to describe the background in terms of it dimensionless characteristic strain defined by

[TABLE]

where the second equality assumes a power-law form for the background. Note that for a background produced by the cosmological population of SMBH binaries, $\alpha=-2/3$.

III.2.1 GW contribution to the noise covariance matrix

Using the above expressions for the timing residual response of a pulsar to a GW background, we can calculate the GW contribution to the noise covariance matrix when cross-correlating timing residuals associated with two Earth-pulsar baselines $I$ and $J$. Denoting the GW contributions to the two sets of timing residuals as $h_{I}(t)$ and $h_{J}(t)$, respectively, one can show that the covariance matrix is block-diagonal with components

[TABLE]

where

[TABLE]

and

[TABLE]

The full noise covariance matrix, which includes contributions instrinsic to the pulsar and to the measurement process, is also block-diagonal with components

[TABLE]

Here $C_{n,I}$ is given by (62), but with the pulsar noise power spectral density $P_{n_{I}}(f)$ replacing $P_{h}(f)$. This last equation assumes that the noise contributions associated with different pulsars are not correlated with one another.

The quantity $\chi_{IJ}\equiv\chi(\zeta_{IJ})$ defined in (61) is the Hellings and Downs factor Hellings and Downs (1983) for a pair of pulsars separated by angle $\zeta_{IJ}=\cos^{-1}(\hat{p}_{I}\cdot\hat{p}_{J})$ (see Figure 8). It arises when cross-correlating the GW-induced timing residuals for an unpolarized, isotropic GW background. Note that $\chi_{IJ}$ has been normalized such that $\chi_{II}=1$ (for a single pulsar).

IV Sensitivity curves

Ultimately, a detection sensitivity curve should tell us how likely it is to detect a particular type of GW signal. So it should depend not only on the properties of the noise in the detector, but also on the type of signal that one is searching for and the method that one uses to search for it. So here we extend the formalism of the previous two sections to define sensitivity curves for searches for a deterministic GW signal from a circular binary and an unpolarized, isotropic stochastic GW background. We begin by writing down expressions for the optimal detection statistics for these two different sources and their corresponding expected signal-to-noise ratios (SNRs). We will see that from these expected SNRs, we can read off an effective strain-noise power spectral density, which has the interpretation of a detection sensitivity curve.

IV.1 Matched filtering for a deterministic GW signal

For a deterministic GW signal, we can use the method of matched filtering to construct an optimal detection statistic. This method has been used extensively in the PTA literature, Babak et al. (2016); Babak and Sesana (2012); Ellis et al. (2012); Yardley et al. (2010) and is also the basis for the approximate deterministic sensitivity curves in Moore et al. (2015b). Letting $Q_{I}$ denote the filter function for pulsar $I$ (where $I=1,2,\cdots,N_{\rm p}$), we define

[TABLE]

where $r_{I}\equiv G_{I}^{T}\delta t_{I}$ are the TMM residuals for pulsar $I$. The filter function is determined by maximizing the expected signal-to-noise ratio, $\rho\equiv\mu/\sigma$, of $\hat{\mathcal{S}}$. The expectation value of $\hat{\mathcal{S}}$ is given by

[TABLE]

and its variance is given by

[TABLE]

where $\Sigma_{n,I}\equiv G_{I}^{T}C_{n,I}G_{I}$ is the noise covariance matrix for $r_{I}$. This result for the variance assumes that the only GW contribution to the timing residuals is from a deterministic GW source, and not from a stochastic GW background. The presence of a stochastic background would contribute to both the diagonal and off-diagonal block matrices (see (64)). In what follows, we will assume that the off-diagonal terms are small compared to the diagonal (auto-correlated) terms. But we will replace $\Sigma_{n,I}$ by $\Sigma_{I}\equiv G_{I}^{T}C_{I}G_{I}$, where $C_{I}\equiv C_{n,I}+C_{h}$, thereby allowing a stochastic background to contribute to the auto-correlated noise (sometimes called GW self-noise).

Using the above results for the mean and variance of $\hat{\mathcal{S}}$, the square of the expected signal-to-noise ratio is

[TABLE]

with the optimal filter given by

[TABLE]

Note that $Q_{I}$ is a noise-weighted version of the TMM signal waveform, as expected for a matched-filter statistic. Using this expression, the expected signal-to-noise ratio becomes

[TABLE]

This last expression can be evaluated in the frequency domain by using (19) for ${\cal N}_{I}^{-1}(f,f^{\prime})$ and restricting to the diagonal component ${\cal N}^{-1}(f)$ as discussed in Section II.4:

[TABLE]

Recall that $\theta$ denote the set of GW parameters. For the case of a circular binary discussed in Section III.1.1, $\theta=\{\hat{k},\iota,\psi,\phi_{0}\}$.

IV.1.1 Detection sensitivity curve for sky and

inclination-averaged sources

To proceed further, we first consider the case of GWs from a single binary system averaged over the initial phase, inclination of the source, as well as its sky location. Using (52) for $|\tilde{h}_{I}(f)|^{2}$, we have

[TABLE]

where

[TABLE]

Here, $S_{I}(f)$ is the strain-noise power spectral density for pulsar $I$, and $S_{\rm eff}(f)$ is an effective strain noise power spectral density for an array of pulsars. Given how $S_{\rm eff}(f)$ appears in the expression for the expected signal-to-noise ratio, we will use it, or its dimensionless characteristic strain,

[TABLE]

as a sensitivity curve for detecting a deterministic GW source averaged over its initial phase, inclination, and sky location. A plot of $S_{\rm eff}(f)$ for the array of pulsars in the NANOGrav 11-year data Arzoumanian et al. (2018) is shown in Figure 9.

Note that for a monochromatic source, $S_{h}(f)$ has a very simple form given by (50), which implies

[TABLE]

IV.1.2 SNR and characteristic amplitude sky maps for

inclination-averaged sources

If we average over initial phase and source inclination, but not over sky location, cf. (48) for $|\tilde{h}_{I}(f,\hat{k})|^{2}$, we obtain

[TABLE]

where

[TABLE]

with $\mathcal{R}_{I}(f,\hat{k})$ given by (49). These expressions are analogous to (73), but with added dependence on the propagation direction $\hat{k}$ of the GW. It turns out that we can factor out the $\hat{k}$ dependence on the right-hand side of the above expression for $S_{\rm eff}(f,\hat{k})$ if we ignore the frequency-dependent part of the pulsar-term contribution to $|R_{I}^{P}(f,\hat{k},0)|^{2}$, as discussed in the context of (36). Making this approximation,

[TABLE]

where $F_{I}^{+,\times}(\hat{k})$ are defined by

[TABLE]

As before, it is easy to do the integral over frequency for a monochromatic source, for which $S_{h}(f)$ is given by (50). The result is

[TABLE]

where the direction $\hat{n}$ of the source on the sky is opposite the direction of GW propagation, $\hat{n}=-\hat{k}$. A plot of $\rho(\hat{n})$ for a pair of $10^{9}$ solar-mass BHs at a luminosity distance of 100 Mpc, emitting monochromatic GWs at the frequency $f_{0}=8\leavevmode\nobreak\ {\rm nHz}$ is shown in Figure 10.

Finally, it is a simple matter to recast the form of the sky map so that we solve (82) for the strain amplitude $h_{0}$ of a monochromatic binary, cf. (44), that would produce a particular value of the signal-to-noise ratio $\rho$:

[TABLE]

A sky map of $h_{0}(\hat{n})$ is shown in panel (a) of Figure 11 for $\rho=2$ using the NANOGrav 11-year data. For comparison, panel (b) shows the actual 95% confidence-level upper limit map taken from the NANOGrav 11-year single-source paper Aggarwal et al. (2018).

IV.2 Single-pulsar characteristic strain noise curves

For an individual pulsar, we will use the characteristic strain

[TABLE]

to characterize its polarization and sky-averaged sensitivity; see (74). Plots of single-pulsar characteristic strain-noise sensitivity curves for the simple quadratic spin-down model described in Section II.1 and for both white and red+white noise are shown in Figure 12.

More realistic single-pulsar strain-noise sensitivity curves can be constructed using a subset of the NANOGrav 11-year pulsars (Figure 13) Arzoumanian et al. (2018).

These pulsars have noise contributions specified by the parameters EQUAD, ECORR, and EFAC Arzoumanian et al. (2018); Lam et al. (2017); Arzoumanian et al. (2014), which are denoted by $Q$, $J_{ij}$, and $F$ in the following expression for the noise covariance matrix:

[TABLE]

Here $\sigma_{i}^{2}$ are individual TOA errors, which are associated with the finite-signal-to-noise ratio determination of the pulse arrival times (obtained by correlating the observed pulses with a pulse template). EQUAD are white noise contributions to the covariance matrix that add in quadrature with the TOA errors. EFAC is an overall scale factor that can be used to adjust the overall uncertainty if necessary. ECORR are noise contributions that are correlated within an observing epoch, but not from epoch to epoch. Hence the ECORR contributions to the covariance matrix are block diagonal. Red noise, modeled as a power law, was added for those pulsars that show significant detections in the NANOGrav 11-year data set Arzoumanian et al. (2018). In Figure 13, B1937+21, J1713+0747 and J1909-3744 have injections of red noise. This can be distinguished by the “flatter” appearance of the sensitivity curves around the minimum, as compared to the other pulsars. For a detailed list of noise parameters, and to see which pulsars have significant detections of red noise, consult Table 2 in Arzoumanian et al. (2018).

The NANOGrav 11-year pulsars also have more complicated timing model fits than the simple quadratic spin-down model described in Section II.1. In Figure 13, one can see that pulsar J1024-0719 is fit to a cubic spin-down model, leading to a steeper frequency-dependence ($\sim\!f^{-5/2}$) at low frequencies. One also sees that J1713+0747 and J1853+1303 are in binary systems: there are additional spikes at the binary orbital frequency and twice the binary orbital frequency for J1853+1303. Finally, these pulsars have timing models that also include fits to a piecewise, time-dependent dispersion measure fluctuation (DMX), which is associated with perturbations of the dispersion of the radio pulses as they propagate through the interstellar medium from the pulsar to a radio receiver on Earth. (The lower-frequency components of a pulse are delayed more than the higher-frequency components.) Fitting to DMX in the timing model leads to broadband absorption of power relative to a timing model that doesn’t fit for DMX. Figure 14 shows plots of the transmission function for NANOGrav pulsar J1944+0907, with and without DMX included in the timing model.

IV.3 Optimal cross-correlation statistic for a

stochastic GW background

The derivation of the optimal cross-correlation statistic for a stochastic GW background is similar to that presented above for a single deterministic GW, expect that we work with data from pairs of pulsars. Starting with a single distinct pair, labeled by $I$ and $J$, we define

[TABLE]

where $r_{I}$ and $r_{J}$ are the TMM residuals for pulsars $I$ and $J$ (assuming that we have already fit for all deterministic GW sources), and $Q$ is an $m_{I}\times m_{J}$ matrix, where $m_{I}\equiv N_{I}-N_{{\rm par},I}$, etc. As before, we determine the filter function $Q$ by maximizing the signal-to-noise ratio of $\hat{\mathcal{S}}_{IJ}$. Similar derivations appear in the literature Allen and Romano (1999); Anholm et al. (2009); Ellis et al. (2012); Chamberlin et al. (2015); Vigeland et al. (2017). The final result for the optimal filter is

[TABLE]

where

[TABLE]

The expected squared signal-to-noise ratio for this optimal choice of $Q$ is then

[TABLE]

The above calculation assumes that we are in the weak-signal limit where the cross-correlation terms are assumed to be negligible compared to auto-correlation terms (i.e., we assume that the GW signal power is much less than that for the intrinsic pulsar and measurement noise).

We can then combine the signal-to-noise ratios for each distinct pair in quadrature since, in the weak-signal limit, there is negligible correlation between these estimators:

[TABLE]

As we saw for deterministic GWs, it is useful to write the above expression for the expected squared signal-to-noise ratio in the frequency domain. Proceeding as we did there, we find

[TABLE]

where $P_{h}(f)=\mathcal{R}(f)S_{h}(f)$, and where $\mathcal{N}_{I}^{-1}(f)$ is defined by (20). This suggests defining the following effective strain-noise power spectral density for the whole PTA:

[TABLE]

which includes contributions from the Hellings and Downs factors $\chi_{IJ}$ and the individual pulsar strain-noise power spectral densities $S_{I}(f)\equiv 1/(\mathcal{N}^{-1}_{I}(f)\mathcal{R}(f))$. Note that $S_{\rm eff}(f)$ has dimensions of strain2/Hz, and that

[TABLE]

in terms of $S_{\rm eff}(f)$.

A plot comparing dimensionless charateristic strain curves $h_{c}(f)\equiv\sqrt{fS_{\rm eff}(f)}$ for stochastic GW backgrounds for the NANOGrav 11-year pulsars is given in Figure 15. The three curves show the effect of including a contribution from the GWB to the auto-power spectra of all the pulsars (blue versus dashed-orange curves) and the false improvement in sensitivity that arises if one fails to include the red-noise component of the individual pulsar noise covariance matrices (green versus dashed-orange curves). Typical PTA sensitivity curves that one sees in the literature incorrectly ignore this red noise component.

IV.3.1 Comparing stochastic and deterministic sensitivity curves

Although one uses different statistics to search for deterministic and stochastic GW signals, it is interesting to compare the sensitivity curves for these two different cases. Figure 16 shows plots of the deterministic and stochastic sensitivity curves for the NANOGrav 11-year pulsars (taken from Figure 9 and Figure 15, dashed-orange curve). Note that the sensitivity curve for a single deterministic source is lower than that for a stochastic background, since the Hellings and Downs factors $\chi_{IJ}$ in (92) reduce the effective number of pulsar pairs that contribute to the stochastic analysis.

To demonstrate this explicitly, compare equations (73) and (92) for $S_{\rm eff}(f)$ assuming that all the pulsars have the same noise characteristics and timing model fits (i.e., $S_{I}(f)\equiv S(f)$ for all $I$), and that all the pulsars are observed for the full observation time (i.e., $T_{I}\equiv T_{IJ}\equiv T_{\rm obs})$. Then

[TABLE]

where $N_{\rm p}$ is the number of pulsars. Since the maximum value of $\chi_{IJ}$ for any pair of pulsars is $1/2$, we have

[TABLE]

which implies

[TABLE]

Thus,

[TABLE]

Although we have compared the full sensitivity curves $S_{\rm eff}(f)$ for deterministic and stochatic GW sources, we note that the corresponding signal-to-noise ratio for a monochromatic deterministic source uses only the value of the sensitivity curve at a single frequency $f=f_{0}$ (see (76)); while that for a stochastic source involves an integral of $S_{\rm eff}(f)$ over all $f$ (see (93) and the discussion in Section IV.3.3).

IV.3.2 Pairwise stochastic sensitivity curves

As a by-product of the stochastic sensitivity curve analysis, we obtain pairwise stochastic sensitivity cuves

[TABLE]

by simply restricting ourselves to a single term in the sum (92). Plots of such curves are useful as a diagnostic for comparing the contribution of different pulsar pairs to the stochastic optimal statistic signal-to-noise ratio. Figure 17 shows pairwise sensitivity curves for a subset of the NANOGrav 11-year pulsars, comparing pairwise correlations of some of the most and least sensitive NANOGrav pulsars.

IV.3.3 Power-law integrated sensitivity curves

For stochastic backgrounds that have a power-law spectrum, cf. (59), it is possible to construct a sensitivity curve that takes into account the improvement in sensitivity that comes from integrating over frequency Thrane and Romano (2013). Given a range of power-law indices, one determines the amplitude of each power-law background that yields a prescribed value of the optimal statistic signal-to-noise ratio $\rho$ (e.g., $\rho=1$). The envelope of these power-law backgrounds defines the power-law-integrated sensitivity curve for the PTA. Figure 18 shows the $\rho=1$ power-law integrated sensitivity curve for the NANOGrav 11-year data set using the dashed-orange characteristic strain-noise curve from Figure 15.

V Discussion

We have presented a method for constructing realistic detection sensitivity curves for pulsar timing arrays, valid for both deterministic and stochastic GW signals. We can include different noise characteristics and the effect of fitting to a timing model via an inverse-noise-weighted transmission function ${\cal N}_{I}^{-1}(f)\approx{\cal T}_{I}(f)/{P(f)}$. Single-pulsar sensitivity curves are then calculated from the strain-noise power spectral density $S_{I}(f)\equiv 1/({\cal N}^{-1}_{I}(f){\cal R}(f))$, where ${\cal R}(f)$ is the polarization and sky-averaged timing residual response of a pulsar to a passing GW. Detection sensitivity curves for multiple pulsars (i.e., a PTA) are similary constructed from an effective strain-noise power spectral density $S_{\rm eff}(f)$, which is a combination of single-pulsar strain-noise power spectral densities $S_{I}(f)$, cf. (73), (78), (92), appropriate for the GW source that one is interested in detecting.

The sensitivity curves that we have calculated can be used to assess the detectability of different GW signals by exisiting or planned PTAs. The computational cost of producing these sensitivity curves is minimal; they can be calculated much faster than doing Monte Carlo simulations using injected signals. By properly incorporating realistic noise properties and the effect of timing model fits into the sensitivity curves, we can produce detectability estimates that agree quite well with the more-computationally-involved predictions.

Acknowledgements

JSH and JDR acknowledge subawards from the University of Wisconsin-Milwaukee for the NSF NANOGrav Physics Frontier Center (NSF PFC-1430284). JDR also acknowledges support from start-up funds from Texas Tech University. TLS acknowledges support from NASA 80NSSC18K0728 and the Hungerford Fund at Swarthmore College. Finally, we thank Robert Caldwell, Rutger van Haasteren and Xavi Siemens for useful discussions and Justin Ellis for sharing some preliminary code.

Appendix A Casting the Blandford et al. analysis Blandford et al. (1984)

in more modern notation

When using pulsar timing data to search for GWs, one needs to take into account the effects of fitting to a deterministic timing model when doing any type of additional signal analysis. Following Blandford et al. (1984), we define the residuals $R\left(t\right)$ as the difference between the observed arrival times of the pulses and the expected arrival times as determined by our best guesses to the parameters. These residuals are fit to an expression linear in the corrections to the unknown parameters, $\alpha_{a}$.444In our notation, $R(t)$ is $\delta t_{i}$ and $\alpha_{a}$ is $\delta\xi_{a}$. (Noise terms are added later in their analysis.) We start in the notation of Blandford et al. (1984), and then translate to expressions in terms of modern PTA GW analyses:

[TABLE]

We will define $R_{i}\equiv R\left(t_{i}\right)$, which is a vector of length $N$, and $\psi_{ia}\equiv\psi_{a}\left(t_{i}\right)$, which is a 2-dimensional matrix with dimensions $N\times N_{{\rm par}}$. (Note we have reversed the order of the indices on $\psi_{ia}$ from that in Blandford et al. (1984), to be consistent with later work.) In more modern PTA data analysis papers, like van Haasteren and Levin (2013) or Ellis et al. (2013); Chamberlin et al. (2015), this matrix is referred to as the design matrix of the timing model (our $M_{ia}$.) The above expression for the residuals can be transformed into an orthonormal basis

[TABLE]

where

[TABLE]

Using these definitions we calculate a relation that will be useful in the next section. To simplify the notation a bit we will use the Einstein convention of summing over repeated indices without including summation symbols, using matrix transposes where necessary. Thus, for example, the orthonormality conditions can be written as

[TABLE]

Since a change of basis change is invertible, we can act with the inverse transformation matrices:

[TABLE]

where $L^{-T}$ denotes the inverse of the transpose matrix $L^{T}$, which is the same as the transpose of the inverse matrix $L^{-1}$. Finally, using the well-known identity for the inverse of a product of two matrices:

[TABLE]

A.1 Least-squares regression

One finds the best fit to a timing model by minimizing a $\chi^{2}$ function, which we will define below. In Blandford et al. (1984) an ordinary least squares (OLS) minimization is used. In subsequent PTA papers a weighted-least-squares (WLS) regression is used, where each residual is weighted by the inverse of the TOA error, $W_{i}\equiv 1/\sigma_{i}$. In the most modern work a generalized least squares (GLS) regression is used where the noise covariance matrix, $N_{ij}$, is used, encoding covariances between all residuals:

[TABLE]

Here we solve the GLS minimization problem, restricting to simpler scenarios if needed—i.e., $N^{-1}_{ij}=\sigma_{i}^{-2}\,\delta_{ij}$ for the case of WLS, and $N^{-1}_{ij}=\delta_{ij}$ for OLS (as noise is not taken into account during the OLS fit). We minimize the expression for $\chi^{2}$ above by finding the root(s) of the derivative with respect to the parameters:

[TABLE]

Solving for $\alpha_{b}^{\prime}$ gives

[TABLE]

In Blandford et al. (1984), they consider OLS fitting. There the noise is taken into account after the fit, but its existence is implicit throughout. For instance the difference between the LHS and RHS side of their Equation (2.9) would be zero if there was no noise. Setting $N_{ij}=\delta_{ij}$ gives

[TABLE]

This is the result that Blandford et al. (1984) reports for the best fit. For WLS fitting, we have

[TABLE]

where $W^{2}_{ij}\equiv\sigma_{i}^{-2}\,\delta_{ij}$.

A.2 Transmission function for ordinary least-squares regression

The transmission function is defined by Blandford et al. (1984) as the transfer function relating the power in the pre-fit residuals $R_{i}$ to that in the post-fit residuals

[TABLE]

where $\alpha_{a}^{\prime}$ are the best-fit values to the parameter deviations, determined by the $\chi^{2}$ minimization procedure discussed above. For the case of OLS fitting, which Blandford et al. (1984) consider, $\alpha_{a}^{\prime}$ is given by (109), implying

[TABLE]

The variance in the post-fit residual is then

[TABLE]

where we used orthogonality of the $\psi^{\prime}_{ja}$ to get the last line. Since the covariance matrix $\langle R_{i}R_{j}\rangle$ is related to its power spectral density $P(f)$ via

[TABLE]

it follows that

[TABLE]

where

[TABLE]

with $\tilde{\psi}_{a}^{\prime}$ the Fourier transforms of the basis functions:

[TABLE]

Making this substitution and transforming $\psi_{ia}^{\prime}$ back to the original basis, we find

[TABLE]

which is an expression for transmission function in terms of the original design matrix $\psi_{ia}$.

A.3 Transmission function for generalized least-squares regression

For the case of GLS fitting, the best-fit values for the timing parameter deviations are given by (108), for which the post-fit residuals are given by

[TABLE]

We can write this in terms of the original basis as

[TABLE]

which has exactly the same form as (119) with $\psi^{\prime}_{ia}$ replaced by $\psi_{ia}$. The variance of the post-fit residuals is thus

[TABLE]

Since $\langle R_{i}R_{j}\rangle\equiv N_{jk}$ for GLS fitting, we get

[TABLE]

where we used the symmetry of $N_{ij}$ throughout. Finally, using (114) for $N_{ij}$, we recover (115) with

[TABLE]

We thus obtain the same $R$-matrix-dependent transmission function ${\cal T}_{R}(f)$ found in (15), with the $R$-matrix given by the expression in parentheses, $R_{ij}\equiv\delta_{ij}-\psi_{ia}\left(\psi_{bl}^{T}N_{lm}^{-1}\psi_{ma}\right)^{-1}\psi_{bk}^{T}N^{-1}_{kj}$.