State-space analysis of a continuous gravitational wave source with a pulsar timing array: inclusion of the pulsar terms

We assume that the rest frame spin frequency of the $n$-th pulsar evolves according to a mean-reverting Ornstein-Uhlenbeck process, described by a Langevin equation with a time-dependent drift term (Vargas & Melatos, 2023),

where $f_{\rm em}^{(n)}$ is the deterministic component of the evolution, an overdot denotes a derivative with respect to $t$, $\gamma^{(n)}$ is a damping constant whose reciprocal specifies the mean-reversion timescale, and $\xi^{(n)}(t)$ is a white noise stochastic process which satisfies

In Equation (3), $[\sigma^{(n)}]^{2}$ parametrizes the noise amplitude and produces characteristic root mean square fluctuations $\approx\sigma^{(n)}/[\gamma^{(n)}]^{1/2}$ in $f_{\rm p}^{(n)}(t)$ (Gardiner, 2009). The deterministic evolution $f_{\rm em}^{(n)}(t)$ is attributed to magnetic dipole braking for the sake of definiteness, with braking index $n_{\rm em}=3$ (Goldreich & Julian, 1969). PTAs are typically composed of millisecond pulsars (MSPs), for which the quadratic correction due to $n_{\rm em}$ in $f_{\rm p}^{(n)}(t)$ is negligible over the observation time $T_{\rm obs}\sim 10\,{\rm yr}$. Consequently, $f_{\rm em}^{(n)}(t)$ can be approximated accurately by

where $t_{1}$ labels the first TOA.

Equations (1)–(4) are not unique. Rather, they offer one possible phenomenological description consistent with the main qualitative features of a typical PTA pulsar’s observed spin evolution, i.e. random, mean-reverting, small-amplitude excursions around a smooth, secular trend (Agazie et al., 2023b; Antoniadis et al., 2023b; Zic et al., 2023). A phenomenological approach is obligatory, because a predictive, first-principles theory of timing noise does not currently exist, c.f. the multiple theorized mechanisms referenced in Section 1 of K24. Langevin equations like (1)–(4) have been applied successfully to analyse anomalous braking indices (Vargas & Melatos, 2023) and in hidden Markov model glitch searches (Melatos et al., 2020; Lower et al., 2021; Dunn et al., 2022; Dunn et al., 2023). However, they are highly idealised (Meyers et al., 2021b; Meyers et al., 2021a; Antonelli et al., 2023; Vargas & Melatos, 2023). Idealisations include (i) the exclusion of physics, which is likely to be present in reality, e.g. the classic, two-component, crust-superfluid structure inferred from post-glitch recoveries (Baym et al., 1969; van Eysden & Melatos, 2010; Gügercinoǧlu & Alpar, 2017; Meyers et al., 2021a; Meyers et al., 2021b); (ii) the exclusion of non-Gaussian excursions such as Lévy flights (Sornette, 2004); (iii) the whiteness of $\xi^{(n)}(t)$ for MSPs in PTAs; and (iv) the formal interpretation of $d^{2}f_{\rm p}/dt^{2}$, noting that $\xi^{(n)}(t)$ in Equation (1) is not differentiable.

In the presence of a GW, the rest-frame spin frequency of the $n$-th pulsar is related to the radio pulse frequency measured by an observer on Earth via a measurement equation,

where $d^{(n)}$ labels the distance to the $n$-th pulsar, $f_{\rm p}^{(n)}$ is evaluated at the retarded time $t-d^{(n)}$, and $\varepsilon^{(n)}(t)$ is a Gaussian measurement noise which satisfies

In Equation (7), $[\sigma_{\rm m}^{(n)}]^{2}$ is the variance of the measurement noise at the telescope. In Equation (5) the measurement function $g^{(n)}(t)$ is given by (e.g. Maggiore, 2018)

where $[q^{(n)}]^{i}$ labels the $i$-th coordinate component of the $n$-th pulsar’s position vector $\bm{q}^{(n)}$, $\Omega$ is the constant angular frequency of the GW, $\bm{n}$ is a unit vector specifying the direction of propagation of the GW, $H_{ij}$ is the spatial part of the GW amplitude tensor, and $\Phi_{0}$ is the phase offset of the GW with respect to some reference time.

Regarding point (i) a pulsar’s sky position varies due to its own orbital motion (if it is located in a binary) as well as the rotation and revolution of the Earth. However, pulsar TOAs are defined relative to the Solar System barycentre, after correcting for the pulsar’s binary motion (if any). The barycentering correction is typically applied when generating TOAs, e.g. with tempo2 (Hobbs et al., 2006; Edwards et al., 2006) and related timing software, and is inherited by $f_{\rm m}^{(n)}(t)$. Some PTA pulsars do have non-negligible proper motions of order $10^{2}$ km s^-1 after barycentering (e.g. Jankowski et al., 2018), but we do not consider this effect in this introductory paper.

Regarding point (ii), there are two timescales to consider. The first is the timescale set by the observation period, $T_{\rm obs}$. The second is the timescale set by the light travel time between the pulsar and the Earth, $T_{\rm light}^{(n)}=d^{(n)}/c$. Regarding the former, studies of SMBHB inspirals in the PTA context show that the GW frequency $f_{\rm gw}$ ($=\Omega/2\pi$) evolves over a short time-scale by an amount (e.g. Sesana & Vecchio, 2010)

where $M_{\rm c}$ is the chirp mass of the SMBHB, $f_{\rm gw}(t=t_{1})$ is the GW frequency at the time of the first observation, and one has typically $T_{\rm obs}\sim 10$ years. A source can be considered monochromatic if $\Delta f_{\rm gw}$ is less than the PTA frequency resolution $1/T_{\rm obs}$ ¹¹1Strictly speaking, the frequency resolution is proportional to $1/T_{\rm obs}$ divided by the signal-to-noise ratio; brighter sources are resolved more accurately.. A vast majority of SMBHBs resolved by PTAs are expected to satisfy $\Delta f_{\rm gw}<1/T_{\rm obs}$ and we are therefore justified in treating the GW source as monochromatic over the $T_{\rm obs}$ timescale as a first approximation (Sesana et al., 2008; Sesana & Vecchio, 2010; Ellis et al., 2012).

Regarding the light travel time, for SMBHBs which are sufficiently massive and have sufficiently high orbital frequencies, the source may undergo non-negligible evolution during $T_{\rm light}^{(n)}$, such that the frequency of the GW which is incident upon the Earth does not equal the frequency of the GW which is incident upon the pulsar. Most SMBHBs detectable with PTAs are not expected to satisfy this condition, but some do; for a PTA composed of pulsars with a mean distance of 1.5 kpc, 78% of simulated SMBHBs detectable with the current IPTA undergo negligible evolution, whilst for the second phase of the Square Kilometre Array this fraction drops to 52%; see Figure 7 in Rosado et al. (2015). In this introductory paper, as a first pass, we focus exclusively on sources which undergo negligible evolution during $T_{\rm light}^{(n)}$, following previous investigations of pulsar-term biases in the context of standard PTA analyses (e.g. Zhu et al., 2016; Chen & Wang, 2022). However, the issue is an important one, so we perform some preliminary tests regarding how sensitively the results depend on the monochromatic assumption in Appendix A. We find that the results do not depend sensitively on the assumption that the GW source has a constant angular frequency, especially for systems with low signal-to-noise ratio (SNR). More extensive testing is deferred to future work, after the state-space method matures sufficiently (e.g. by tracking pulsar phase) to be applied to real, astronomical data.

The model described in Sections 2.1 and 2.2 comprises $5N$ static parameters, that are specific to the pulsars in the array, viz.

It also comprises seven parameters, that are specific to the GW source, viz.

where $h_{0}$ is the characteristic wave strain, $\iota$ is the orbital inclination, $\psi$ is the polarisation angle, $\delta$ is the declination and $\alpha$ is the right ascension. These parameters enter the model through Equation (8), with $H_{ij}=H_{ij}(h_{0},\iota,\psi,\delta,\alpha)$ and $\bm{n}=\bm{n}(\delta,\alpha)$. The complete set of $7+5N$ static parameters is denoted by $\bm{\theta}=\bm{\theta}_{\rm gw}\cup\bm{\theta}_{\rm psr}$. While we assume no prior information about $\bm{\theta}_{\rm gw}$, there are constraints on $\bm{\theta}_{\rm psr}$ from electromagnetic observations; for example estimates of $d^{(n)}$ are accurate to $\sim$ 10$\%$ typically (Cordes & Lazio, 2002; Verbiest et al., 2012; Desvignes et al., 2016; Yao et al., 2017).

The GW modulates radio pulses according to Equation (8). The modulation is composed of an Earth term, proportional to $\cos(-\Omega t+\Phi_{0})$, and a pulsar term, proportional to $\cos\left\{-\Omega t+\Phi_{0}+\Omega\left[1+\bm{n}\cdot\bm{q}^{(n)}\right]d^{(n)}\right\}$. The Earth term describes the GW phase at the observer on Earth, while the pulsar term describes the GW phase at the $n$-th pulsar. The Earth term depends only on the GW source parameters and is common across all pulsars. In contrast the pulsar term is a function of $d^{(n)}$ and $\bm{q}^{(n)}$ and varies between pulsars.

The phases of the pulsar terms are related unpredictably, which is why some published PTA analyses approximate the pulsar terms collectively as a source of self-noise. The pulsar term depends on $d^{(n)}$ which is generally poorly constrained by electromagnetic observations, with uncertainties greater than the typical GW wavelength. Consequently the pulsar term is often — although not always — dropped in standard PTA analyses (e.g. Sesana & Vecchio, 2010; Babak & Sesana, 2012a; Petiteau et al., 2013; Zhu et al., 2015; Taylor et al., 2016; Goldstein et al., 2018; Charisi et al., 2023). Dropping the pulsar term is also convenient computationally because it reduces the dimensionality of the parameter space; the $N$ values of $d^{(n)}$ are not inferred. Within the state-space model described in Section 2, dropping the pulsar terms equates to using a modified measurement equation,

with

However, dropping the pulsar term leads to well-known parameter estimation biases, especially in $\alpha$ and $\delta$, and reduces the detection probability by $\sim 5\%$ (Zhu et al., 2016; Chen & Wang, 2022; Kimpson et al., 2024).

In this paper, we generalize the Kalman filter analysis in K24 to include the pulsar terms. To do so, we define a new parameter

such that Equation (8) becomes

The reparametrisation in terms of $\chi^{(n)}$ treats the pulsar-dependent phase correction $\Omega\left[1+\bm{n}\cdot\bm{q}^{(n)}\right]d^{(n)}$ in the argument of the cosine of the pulsar term as a composite parameter to be inferred for each pulsar. In principle, it is possible to disentangle $\Omega$, $\bm{n}$, and $d^{(n)}$ and infer them individually, if $N$ is large enough; for example, $\bm{n}$ appears independently in the phase of the pulsar term and in the denominator of the first line of Equation (8). In practice, however, the likelihood surface is highly corrugated along the $d^{(n)}$-axis, and the prior on $d^{(n)}$ is broad compared to the wavelength $2\pi/\Omega$, so $d^{(n)}$ and hence $\Omega$ are not identifiable for reasonable $N\lesssim 10^{2}$. By replacing $d^{(n)}$ with $\chi^{(n)}$, we obtain a smooth likelihood function, whose maximum is located efficiently and accurately by the nested sampler. This point is discussed in more detail in Appendix E.

The new parametrisation does not increase the dimension of the parameter space, because we trade $d^{(n)}$ for $\chi^{(n)}$. Once $\Omega$, ${\bm{n}}$, and $\chi^{(n)}$ are inferred, it is possible to solve for $d^{(n)}$, albeit not uniquely; $d^{(n)}$ can be inferred up to an integer multiple of the Doppler-shifted wavelength, due to the ${\rm mod}\,2\pi$ operation in Equation (16). However, this ambiguity is unavoidable, whether we replace $d^{(n)}$ with $\chi^{(n)}$ or not, and occurs in every PTA analysis. The static parameters specific to the pulsars in the array (c.f. Equation (10)) are now

whilst $\bm{\theta}_{\rm gw}$ remains unchanged. We define the complete set of $7+5N$ static parameters to be inferred as $\bm{\theta^{\prime}}=\bm{\theta}_{\rm gw}\cup\bm{\theta^{\prime}}_{\rm psr}$.

Figure 1 displays the posterior distribution of $\bm{\theta}_{\rm gw}$ for ten noise realisations in the form of a traditional corner plot for the representative SMBHB source parametrized in the top portion of Table 1. The histograms are one-dimensional posteriors, marginalized over the six other parameters. Each coloured curve corresponds to a different noise realisation. The solid orange line marks the known injected value. The two-dimensional contours mark the (0.5, 1, 1.5, 2)-sigma level surfaces. All histograms and contours are consistent with a unimodal joint posterior, which peaks near the known, injected values. There is scant evidence of railing against the prior bounds. There is no strong evidence for correlations between parameter pairs, e.g. banana-shaped contours, although weak correlations are evident between $\Omega$ and $\Phi_{0}$, $\psi$ and $\alpha$, and $\iota$ and $h_{0}$. The injected value falls within the 90% credible interval of the one-dimensional marginalized posteriors in 60 out of the 70 possible combinations of the seven parameters and 10 realizations. Indeed, the injected value falls within the 90% credible interval for $\Omega,\Phi_{0},\psi$, and $\iota$ 40 out of 40 times.

There is appreciable natural dispersion in the one-dimensional posterior medians between realisations. We quantify the degree of dispersion using the coefficient of variation,

where $\mu_{\rm med}$ is the mean of the 10 posterior medians and $\sigma_{\rm med}$ is their standard deviation. The minimum and maximum $CV$ are 0.2% and 52% for $\Omega$ and $\Phi_{0}$ respectively. The remaining parameters typically satisfy $CV$ $\lesssim 10\%$. Analogous results for synthetic data with higher strain ($h_{0}=1\times 10^{-12}$) are presented in Appendix G (see Figure 12) for completeness. The dispersion between noise realisations decreases as $h_{0}$ increases, and the one-dimensional posteriors for $\iota$ and $h_{0}$ are more symmetric about the injected value.

It is important to check how well the pulsar-term phase $\chi^{(n)}$ is estimated, since the reparameterization involving $\chi^{(n)}$ is a key feature of our approach, as explained in Section 3.2. In short, it turns out that $\chi^{(n)}$ is estimated accurately for all $n$, except under special circumstances. Figure 2 displays the results as a corner plot for the representative subset $\chi^{(1)},\dots,\chi^{(5)}$ over the ten noise realisations, in the same style as Figure 1. We do not display the corner plot for $\chi^{(6)}$, $\dots$, $\chi^{(47)}$ because it is too big. With the exception of the results for $\chi^{(2)}$ (discussed below), all histograms and contours are consistent with a unimodal joint posterior, which peaks near the known, injected values. There is no evidence for correlations between parameter pairs. Most phases are recovered unambiguously across all noise realisations; for the 47 $\chi^{(n)}$ parameters (i.e. not just the five plotted in Figure 2), the injected value is contained within the 90% credible interval of the one-dimensional marginalized posteriors in 383 out of the 470 possible combinations of the 47 parameters and ten realizations. The next paragraph explains why 383 out of 470 is fewer than 90%.

Sometimes, albeit rarely, $\chi^{(n)}$ is not estimated consistently across multiple realizations for some $n$, e.g $\chi^{(3)}$ (third row, third column of Figure 2) corresponding to PSR J0340+4130. Out of all the pulsars in the array, this object matches most closely the direction to the synthetic GW source, with $\bm{n}\cdot\bm{q}^{(2)}=-0.96$ and $\cos\chi^{(2)}\approx 0.07$. The ability to accurately infer $\chi^{(2)}$ improves with the SNR. For example, for $h_{0}=1\times 10^{-12}$, $\chi^{(2)}$ is recovered unambiguously; see Figure 13 in Appendix G.

Similar results are obtained for the 4$N$ parameters in $\bm{\theta}_{\rm psr}^{\prime}$ besides $\chi^{(n)}$. The injected values are recovered unambiguously across the 10 noise realisations. For the sake of brevity we do not display the calculated posterior distributions, because inferring $\bm{\theta}_{\rm gw}$ is the main focus of this paper and most published PTA analyses. In short, the estimates of $f_{\rm em}(t_{1})$ and $\dot{f}_{\rm em}(t_{1})$ are guided into narrow ranges by the narrow priors. The one-dimensional posteriors inferred for $\sigma^{(n)}$ are generally broader due to the broader prior, but contain the injection within the 90% credible interval in a majority of cases. We remind the reader that $\gamma^{(n)}$ is not estimated in this paper; typically it satisfies $\gamma^{(n)}\sim 10^{-5}T_{\rm obs}$ (Price et al., 2012; Meyers et al., 2021a; Meyers et al., 2021b; Vargas & Melatos, 2023), so its influence is muted.

Section 5.1 focuses on a single representative source, summarised in Table 1. In this section we test the method for a range of sources, varying $\bm{\theta}_{\rm gw}$. The aim is to verify that the analysis scheme works across an astrophysically relevant domain and that the arbitrary choice of $\bm{\theta}_{\rm gw}$ in Table 1 is not advantageous by accident.

We analyse 200 injections constructed by fixing $h_{0}=5\times 10^{-15}$ and $\iota=1.0$ rad and drawing the remaining five elements of $\bm{\theta}_{\rm gw}$ randomly from the prior distributions defined in Table 1. We fix $h_{0}$ and $\iota$ in order to maintain an approximately constant SNR across the 200 injections. For each injection we compute the posterior distribution of ${\bm{\theta}}_{\rm gw}$. The static pulsar parameters $\bm{\theta^{\prime}}_{\rm psr}$ are specified in Section 4 and Table 1. It is prohibitive to display corner plots analogous to Figure 1 for 200 injections and seven elements of $\bm{\theta}_{\rm gw}$. Therefore, we summarise the results with the aid of a parameter-parameter (PP) plot (Cook et al., 2006). A PP plot displays the fraction of injections included within a given credible interval of the estimated posterior, as a function of the credible interval itself. In the ideal case of perfect recovery, the PP plot should be a diagonal line of unit slope.

The PP results are displayed in Figure 3(a). The shaded grey contours enclose the $1\sigma$, $2\sigma$, and $3\sigma$ significance levels for 200 injections. The curves for the five parameters, colour-coded in the legend, are approximately linear. Some parameters are consistently well estimated across the parameter domain, e.g. the blue curve for $\Omega$ lies everywhere within the 2$\sigma$ shaded region. Conversely some parameters, e.g. $\Phi_{0}$, stray outside the $1\sigma$ shaded region and show evidence of being slightly over-constrained; there are more injections contained within lower-value credible intervals than would be expected statistically, and fewer injections contained within higher-value credible intervals. This occurs due to the low SNR of the injected GW signal and the varying sensitivity of the specific PTA configuration as a function of sky position. More quantitatively, for the best estimated parameter $\Omega$ the injection is contained within the 90% credible interval in 91% of cases. For the worst estimated parameter $\Phi_{0}$, the injection is contained within the 90% credible interval in 81% of cases. The next paragraph explains why 81% is less than the theoretical ideal 90%.

Manually inspecting the individual corner plots for each of the 200 injections, we see that the nested sampler does not return a unimodal posterior in the rare event that the source is located unfavourably, with $\bm{n}\cdot\bm{q}^{(n)}\approx-1$, such that one cannot infer $\chi^{(n)}$ accurately (see Section 5.2). All injected sources are “observed” with the same settings such as $T_{\rm obs}$, $T_{\rm cad}$ and $n_{\rm live}$. In Figure 3(b) we display a second PP plot, arranged identically to Figure 3(a), but assuming that the injected values of $\chi^{(n)}$ are known exactly for the sake of testing, i.e. setting a delta-function prior on $\chi^{(n)}$. The excursions out of the shaded region diminish compared to Figure 3(a), confirming the importance of estimating $\chi^{(n)}$ accurately. For the best estimated parameter $\Omega$ the injection is contained within the 90% credible interval in 91% of cases, the same as in Figure 3(a). For the worst estimated parameters such as $\Phi_{0}$, the injection is contained within the 90% credible interval in 89% of cases, an improvement over the results of Figure 3(a).

In this section we apply the Kalman filter in conjunction with nested sampling in order to infer the joint posterior distribution for the static parameters with and without the pulsar terms. For $\mathcal{M}_{\rm Earth}$ we apply the Kalman filter using Equation (15) to return $\mathcal{L}(\bm{Y}|\bm{\theta}_{\rm Earth})$ and the nested sampler to estimate $p(\bm{\theta}_{\rm Earth}|\bm{Y})$. For $\mathcal{M}_{\rm psr\&Earth}$ we apply the Kalman filter using Equation (17) to return $\mathcal{L}(\bm{Y}|\bm{\theta}^{\prime})$ and the nested sampler to estimate $p(\bm{\theta^{\prime}}|\bm{Y})$. The settings of the nested sampler (for example the number of live points; see Appendix C) are identical for both models.

We consider two representative systems. The first is a “low-SNR” system with $h_{0}=5\times 10^{-15}$, i.e. the system described in Table 1. The second is a “high-SNR” system which has the same static parameters as in Table 1 except with $h_{0}=1\times 10^{-12}$. The “high-SNR” system is considered in order to quantify any systematic biases, as distinct from random errors caused by the measurement noise. In order to enable a clear comparison between the posterior probability distributions calculated using the two models, in this section we present a single noise realisation of the synthetic data $\bm{Y}$. The results quoted are consistent across different noise realisations.

The results for the seven parameters in $\bm{\theta}_{\rm gw}$ are shown in Figure 4 for the high-SNR system and in Figure 5 for the low-SNR system. The corner plot is arranged identically to Figure 1, except that the different coloured curves now correspond to different inference models, rather than different realisations of $\bm{Y}$. The blue curves are the results derived using $\mathcal{M}_{\rm Earth}$. The green curves are the results derived using $\mathcal{M}_{\rm psr\&Earth}$.

For the high-SNR results in Figure 4, the one-dimensional posteriors for $\mathcal{M}_{\rm Earth}$ are biased, as was observed by K24, as well as Zhu et al. (2016) and Chen & Wang (2022). Particular biases are observed in $\psi,\iota$ and $\alpha$, with the blue posterior displaced with respect to the orange injection line. For example, for $\iota$, the injected value is contained within the 90% credible interval, but the median of the posterior of the Earth-term model is shifted from the injected value by 0.35 radians. The inclusion of the pulsar terms corrects for this bias. The green one-dimensional posteriors exhibit no bias and are generally symmetric about the injected value. Moreover, for every pair of parameters in ${\bm{\theta}}_{\rm gw}$, the injected values are contained within the 2-sigma contours for $\mathcal{M}_{\rm psr\&Earth}$. This is not true for $\mathcal{M}_{\rm Earth}$, where the injected values fall outside the 2-sigma contour for 15 out of the 21 parameter pairs.

We define a relative error to quantify the accuracy of the one-dimensional posteriors with respect to the injected value. The relative error is defined to be the relative unsigned displacement of the mode of $p_{\rm M}\left(\theta|\bm{Y}\right)$ from the injection, viz.

In Equation (22), $p_{\rm M}\left(\theta|\bm{Y}\right)$ is the one-dimensional posterior for parameter $\theta$ returned by nested sampling (c.f. Equation (12)). The subscript $\rm M\in\{\text{Earth, psr\&Earth}\}$ indicates whether the posterior is estimated using Equation (15) or Equation (8) respectively. $\theta_{\rm inj}$ denotes the true injection value. We note that $\Delta_{\rm M}(\theta)$ quantifies the accuracy of the estimates returned by the two models, i.e. the closeness of the most probable estimate of $\theta$ relative to $\theta_{\rm inj}$. In contrast, $\Delta_{\rm M}(\theta)$ does not quantify the uncertainty in the estimates. The error $\Delta_{\rm M}(\theta)$ for each parameter in Figure 4 is summarised in Table 2. We find that the estimates are more accurate using the pulsar terms, i.e. $\Delta_{\rm Earth}(\theta)>\Delta_{\rm psr\&Earth}(\theta)$ for $\theta\in\bm{\theta}_{\rm gw}$. In some cases the difference is modest; $\Omega$ is recovered with high accuracy by both models, with $\Delta_{\rm Earth}(\Omega)-\Delta_{\rm psr\&Earth}(\Omega)=1.5\times 10^{-5}$. The improvement from including the pulsar terms is largest for $\iota$, with $\Delta_{\rm Earth}(\iota)-\Delta_{\rm psr\&Earth}(\iota)=0.14$. Whilst we present only a single noise realisation in this section, the improvements from including the pulsar terms are found to be comparable across different realisations.

For the low-SNR results in Figure 5, there is less improvement in the parameter estimates. Qualitatively, the green and blue contours and histograms mostly overlap, although the green contours are centred slightly better on the injected values. The relative error $\Delta_{\rm M}(\theta)$ for the low-SNR results is reported in the lower half of Table 2. We see that the error is larger for every element in $\bm{\theta}_{\rm gw}$ than in the high-SNR case, as expected. The inclusion of the pulsar terms improve the estimates for five out of the seven static parameters; we find $\Delta_{\rm Earth}(\theta)>\Delta_{\rm psr\&Earth}(\theta)$ for all parameters except $\psi$ and $h_{0}$. However it is hard to draw strong conclusions; the improvements from including the pulsar terms vary randomly across different noise realisations. We show in Section 6.2 that including the pulsar terms increases the detection probability.

In this section we compute the minimum detectable strain for the representative source in Table 1, using $\mathcal{M}_{\rm Earth}$ and $\mathcal{M}_{\rm psr\&Earth}$. We frame the detection problem in terms of a Bayesian model selection procedure, following the lead of other PTA analyses (e.g. Agazie et al., 2023a; Antoniadis et al., 2023a; Reardon et al., 2023; Xu et al., 2023). We define $\mathcal{M}_{\rm null}$ as the null model that assumes no GW exists in the data. This is equivalent to setting $g^{(n)}(t)=1$ in Equation (8). The evidence integral $\mathcal{Z}$ returned by nested sampling, Equation (13), is the probability of the data $\bm{Y}$ given a model $\mathcal{M}_{\rm M}$. The support in the data for the presence of a GW signal, described by model $\mathcal{M}_{\rm M}$, over the absence of a GW signal is quantified via the Bayes factor

In this paper we consider $\rm M\in\left\{\rm Earth,psr\&Earth\right\}$.

The Bayes factors $\beta_{\rm Earth}$ and $\beta_{\rm psr{\&}Earth}$ are plotted as functions of $h_{0}$ in Figure 6 for the representative source in Table 1. We vary the source amplitude from $h_{0}=10^{-15}$ (undetectable) to $h_{0}=10^{-12}$ (easily detectable). To sharpen the comparison between $\mathcal{M}_{\rm psr\&Earth}$ and $\mathcal{M}_{\rm Earth}$ we present only a single noise realisation of the synthetic data $\bm{Y}$, as in Section 6.1; the conclusions drawn below are consistent across different noise realisations. Moreover, the noise processes in the synthetic data are identical realisations for each value of $h_{0}$; the only change from one $h_{0}$ value to the next is $h_{0}$ itself, to smooth the curves in Figure 6.1.

Figure 6 reveals an approximate quadratic relationship $\beta\propto h_{0}^{2}$ for $h_{0}\gtrsim 10^{-14}$ for both $\mathcal{M}_{\rm Earth}$ and $\mathcal{M}_{\rm psr\&Earth}$. Moreover, we obtain $\beta_{\rm psr\&Earth}>\beta_{\rm Earth}$ for all $h_{0}$. That is, for a given $h_{0}$, $\mathcal{M}_{\rm psr\&Earth}$ provides greater evidence for a GW signal in the noisy data than $\mathcal{M}_{\rm Earth}$. The GW source is detectable with decisive evidence ($\beta\geq 10$) for $h_{0}\gtrsim 3.2\times 10^{-15}$ for $\mathcal{M}_{\rm Earth}$ and $h_{0}\gtrsim 2.8\times 10^{-15}$ for $\mathcal{M}_{\rm psr\&Earth}$, a relative improvement in the minimal detectable strain of $14\%$. The minimum detectable strain and the relative improvement through using $\mathcal{M}_{\rm psr\&Earth}$ are particular to the system in Table 1 and the realisation of $\bm{Y}$. They are influenced in general by $T_{\rm obs}$, $T_{\rm cad}$, ${\bm{\theta}}_{\rm gw}$, and ${\bm{\theta}}_{\rm psr}$. A full parameter sweep is postponed to future work, after the analysis scheme in this paper is upgraded from ingesting pulse frequencies to pulse TOAs to enable a like-for-like comparison with standard PTA analyses. In the interim, we note that a $14\%$ improvement in sensitivity when including the pulsar terms is comparable to the improvements of $5\%$ achieved by Zhu et al. (2016).