Robust Estimation of Scattering in Pulsar Timing Analysis

TL;DR

This paper introduces a robust, profile-domain method for modeling time-variable scattering in pulsar timing, improving the accuracy of dispersion measure and timing estimates, especially at low frequencies, by jointly estimating scattering variations.

Contribution

It presents a new approach and GPU-accelerated code for simultaneously modeling scattering, DM, and pulse profile evolution in pulsar timing analysis.

Findings

Ignoring scattering variability biases timing parameters.

Including scattering models reduces uncertainties in DM and timing.

Joint modeling improves parameter estimates at low frequencies.

Abstract

We present a robust approach to incorporating models for the time-variable broadening of the pulse profile due to scattering in the ionized interstellar medium into profile-domain pulsar timing analysis. We use this approach to simultaneously estimate temporal variations in both the dispersion measure (DM) and scattering, together with a model for the pulse profile that includes smooth evolution as a function of frequency, and the pulsar's timing model. We show that fixing the scattering timescales when forming time-of-arrival estimates, as has been suggested in the context of traditional pulsar timing analysis, can significantly underestimate the uncertainties in both DM, and the arrival time of the pulse, leading to bias in the timing parameters. We apply our method using a new, publicly available, GPU accelerated code, both to simulations, and observations of the millisecond pulsar PSR J1643$-$1224. This pulsar is known to exhibit significant scattering variability compared to typical millisecond pulsars, and we find including low-frequency ($< 1$ GHz) data without a model for these scattering variations leads to significant periodic structure in the DM, and also biases the astrometric parameters at the $4\sigma$ level, for example, changing proper motion in right ascension by $0.50 \pm 0.12$. If low frequency observations are to be included when significant scattering variations are present, we conclude it is necessary to not just model those variations, but also to sample the parameters that describe the variations simultaneously with all other parameters in the model, a task for which profile domain pulsar timing is ideally suited.

Tables

Table 1. Table 1: Timing Parameters for PSR J1643 − - 1224

Measured Quantities
Scenario	S1	S2	S3
Right ascension, (hh:mm:ss)	16:43:38.161512(21)	16:43:38.161572(11)	16:43:38.161510(12)
Declination, (dd:mm:ss)	$-$12:24:58.6724(14)	$-$12:24:58.6725(8)	$-$12:24:58.6734(9)
Pulse frequency, (s-1)	216.373337142644(5)	216.373337142639 (3)	216.373337142639(4)
First derivative of pulse frequency, (s-2)	$-$8.6450(5)$\times {10}^{-16}$	$-$8.6448(3)$\times {10}^{-16}$	$-$8.6450(4)$\times {10}^{-16}$
Proper motion in right ascension, (mas yr-1)	6.00(10)	5.50(6)	5.85(8)
Proper motion in declination, (mas yr-1)	3.8(5)	3.5(3)	3.9(4)
Parallax, $\pi$ (mas)	1.3(4)	1.3(3)	1.0(3)
Orbital period, $P_b$ (d)	147.01739774(6)	147.01739775(4)	147.01739772(5)
Epoch of periastron, $T_0$ (MJD)	49283.9336(5)	49283.9337(4)	49283.9336(5)
Projected semi-major axis of orbit, $x$ (lt-s)	25.0726165(29)	25.0726182(19)	25.0726165(23)
Longitude of periastron, ${\omega }_0$ (deg)	321.8487(13)	321.8489(10)	321.8485(11)
Orbital eccentricity, $e$	5.05755(12)$\times {10}^{-4}$	5.05757(10)$\times {10}^{-4}$	5.05752(11)$\times {10}^{-4}$
First derivative of $x$, $\dot{x}$ (${10}^{-12}$)	$-$5.2(5)$\times {10}^{-14}$	$-$5.5(3)$\times {10}^{-14}$	$-$5.2(4)$\times {10}^{-14}$
Set Quantities
Epoch of frequency determination (MJD)	55000	55000	55000
Epoch of position determination (MJD)	55000	55000	55000
Epoch of DM determination (MJD)	55000	55000	55000

Equations

$$B_n(t;\Lambda )\equiv {\Lambda }^{-1/2}{\Phi }_n({\Lambda }^{-1}t),$$

$${\Phi }_n(t)\equiv {\left[2^{n}n!\sqrt{\pi }\right]}^{-1/2}{H}_n\left(t\right)\exp \left(-\frac{t^2}{2}\right),$$

$${\tilde{\Phi }}_n(f)=i^n{\Phi }_n(f)$$

$${\tilde{B}}_n(f;\Lambda )=i^n{B}_n(f;{\Lambda }^{-1})$$

$$\tilde{s}(f,\zeta ,\Lambda )=\sum _{n=0}^{n_{\max }}{\zeta }_n(\nu ){\tilde{B}}_n(f;{\Lambda }^{-1}),$$

$${\zeta }_n(\nu )=\sum _{k=0}^p(\nu -{\nu }_c{)}^k{\zeta }_{n,k}(\nu ),$$

$$\tilde{s}(f,A,\zeta ,\Lambda )=A\sum _{n=0}^{n_{\max }}{\zeta }_n(\nu ){\tilde{B}}_n(f;\Lambda ),$$

$$R(f,\delta \varphi )=\exp \left(2\pi if\delta \varphi \right),$$

$$\mathrm{PBF}(t,\bar{\tau },\nu ,\alpha )=H(t)\exp \left(-\frac{t}{{\nu }^{-\alpha }10^{\bar{\tau }}}\right).$$

$$\mathrm{PBF}(f,\bar{\tau },\nu ,\alpha )=\frac{1}{(2\pi f{\nu }^{-\alpha }10^{\bar{\tau }}{)}^2+1}+\frac{-2\pi f{\nu }^{-\alpha }10^{\bar{\tau }}}{(2\pi f{\nu }^{-\alpha }10^{\bar{\tau }}{)}^2+1}i,$$

$$\frac{d\mathrm{PBF}}{d\bar{\tau }}$$

$$\frac{d\mathrm{PBF}}{d\alpha }$$

$$\mathrm{Pr}(\tilde{\mathbf{d}}|\theta )$$

Full text

1 Introduction

The eventual detection of gravitational waves in the nanoHertz window using a pulsar timing array (Foster & Backer, 1990) will require a thorough understanding of the myriad mechanisms that can impact either the shape, or the time of arrival (ToA) of the pulses of light from those pulsars. The Ionised Interstellar Medium (IISM) is known to be the dominant such mechanism in pulsar timing experiments (e.g., Lam et al. 2016), and introduces both changes in the shape of the profile, and delays in the arrival times.

The delays in the arrival times are primarily the result of interstellar dispersion. As the pulse propagates through the ionised plasma that makes up the IISM, it interacts with free electrons causing a frequency dependent delay. This delay is proportional the the integrated column density of electrons along our line of sight to the pulsar, called the dispersion measure (DM), and scales as $\nu^{-2}$, with $\nu$ the observing frequency. As our line of sight to the pulsar changes with time, so too does the observed column density, introducing a time-variable delay in arrival times.

Over the last several years significant progress has been made in modeling these variations in DM, and in propagating our uncertainties in that model through to the timing parameters, with multiple different methods in use by groups around the world (e.g. Keith et al. 2013; Demorest et al. 2013; Lentati et al. 2014).

A more subtle effect due to the IISM is that of scattering, where inhomogeneities in the IISM cause both intensity variations (called diffractive scintillation, e.g. Narayan 1992) and also a broadening of the pulse profile. The impact of scattering, as for DM variations, will also change with time as our line of sight to the pulsar changes, and in this paper we will be concerned with how to robustly incorporate scattering variations into pulsar timing analysis.

Given the complexity of the IISM, and the steep frequency dependence of the impact it has on pulsar timing, one can rightly ask why we should not just migrate to high-frequency observations alone. Indeed, this approach has already resulted in the most sensitive limit on an isotropic gravitational wave background using a pulsar timing array to date (Shannon et al., 2015). The simple answer is that pulsars are known to have steep spectral indices ($S_{\nu}\propto\nu^{-1.8}$ on average, Maron et al. 2000), and so are much brighter at lower observing frequencies. Given the strong dependence of the detection probability of gravitational waves using a pulsar timing array on the number of pulsars in the array (e.g., Taylor et al. 2016), being able to include more pulsars by virtue of observing at lower frequencies where they are brightest has clear benefits.

Previous attempts to model scattering have either made significant assumptions about a priori unknown quantities, such as the intrinsic pulse shape, or only consider scattering independent of other aspects of the model such as the intrinsic variability in the pulse shape, or DM variations. In both cases the result is a failure to propagate the uncertainties from these other parameters into the scattering measurements.

For example, Löhmer et al. (2001) use profile data from 4.9 GHz observations to construct a model for the intrinsic pulse profile, and then assume this model when determining the scattering time scales and scaling with frequency. However, significant profile evolution is known to occur across wide frequency ranges which can bias these parameter estimates. In this case the uncertainties on the measured parameters were multiplied by a factor of three, however this is clearly unsatifactory, and a more statistically robust approach would be preferred. In Bhat, Cordes & Chatterjee (2003), no assumptions are made about the intrinsic shape of the pulse profile, however the DM is assumed fixed, and the CLEAN algorithm is sub-optimal when in the scattering timescale is small.

Further complications arise when wanting to include scattering variations in timing analysis. Existing approaches have suggested obtaining an estimate for the scattering time scale, and then using that to ‘correct’ the ToAs (e.g., Levin et al. 2016). While cyclic spectroscopy has shown significant promise, the robustness of the method breaks down for low signal-to-noise profiles (S/N), where the S/N is less than $\sim$ 100 (Palliyaguru et al., 2015). Even in the high S/N simulations, however, the approach advocated is to ‘correct’ the ToAs using estimates of the scattering delay obtained from the cyclic spectroscopy. In the PSR J1643$-$1224 data set analysed in Section 4 the mean S/N per epoch is $\sim$ 90, with a minimum of only 15. We therefore require an approach that is robust across all observed profile S/N. Further, neither approach accounts for the significant covariances that exist between the scattering timescales and other parameters of interest.

We illustrate two of the key challenges associated with incorporating scattering measurements into timing analysis in Fig. 1. We simulate a pulsar with a 4.6 ms period, using a Gaussian pulse shape with a full-width at half-maximum of 5 % pulse phase. The scattering time scale is chosen to be $10^{-4.5}$ s at a reference frequency of 1 GHz, consistent with the values observed in our analysis of PSR J1643$-$1224 in Section 4. We use a frequency range of 1.25-1.5 GHz, separated into eight channels, again chosen to be consistent with a typical 20 cm observation used in Section 4. In the top left panel we show the simulated pulse profile at the lowest-frequency channel (red lines, 1263 MHz), and the highest-frequency channel (black lines, 1475 MHz). We have aligned the leading edge of the profiles (shown as a zoomed in region in the bottom sub-panel), and can see that the trailing edge has been broadened in the low-frequency channel (shown as a zoomed in region in the top sub-panel). In the top-right panel we show the same thing, however after adding noise to the simulation. By eye the scattering is undetectable, however we will now show that it is still sufficient to significantly bias parameter estimates when ignored in the analysis.

In the bottom-left panel, we show in black the one- and two-dimensional marginalised posterior distributions from a three-dimensional analysis of this simulated observation. The parameters included in the analysis are:

• (i)

the ToA, measured in units of phase,

• (ii)

the DM in units of the Tempo2 uncertainty,

• (iii)

and $\log_{10}$ of the scattering timescale, $\tau$, measured in seconds at a reference frequency of 1GHz.

One can clearly see that the scattering timescale is correlated in a non-linear way with both the pulse ToA, and the DM. The ToA in particular has a bi-modal distribution, with one peak associated with ‘large’ scattering timescales ($\log_{10}\tau\ga-5$), and one associated with small scattering time scales ($\log_{10}\tau\la-5$). Note that the peak of the one-dimensional probability distribution for the scattering time scale is consistent with the simulated value, indicating that despite being practically invisible by eye, it is still marginally detected in our analysis.

In the same panel we show the posterior distributions for the phase and DM when fixing the scattering timescale at $10^{-10}$ (red lines), and $10^{-4.4}$ (blue lines). These correspond to models with no scattering, or using the maximum likelihood scattering timescale obtained from the three-dimensional analysis. In both cases this results in significant bias in both the measured values, and the uncertainties of the other two parameters. In particular, the measured ToA is completely inconsistent between the two models, with a difference of approximately 10 $\mu$s, two to three orders of magnitude more than the shift expected from an isotropic gravitational wave background in pulsar timing observations given the most stringent current upper limits (Shannon et al., 2015). While the frequency dependence of scattering variations will help to decrease their covariance with gravitational waves, and so will limit how much this bias truly impacts our sensitivity, clearly leaving such variations unmodelled in the analysis must be sub-optimal.

In Fig. 1 (bottom-right panel) we then compare the posterior distribution for the ToA from the three-dimensional analysis (black line) to a Gaussian approximation of the posterior (red line). Typically when forming ToAs, a Gaussian approximation is made of the true probability density function of the arrival time. While methods such as those presented in Liu et al. (2014) and Pennucci, Demorest & Ransom (2014), which perform a simultaneous fit to broadband profile data in order to fit for both phase and DM can be extended to also fit for scattering, the fundamental approximation of Gaussianity implicit in the ToA forming process is invalid in regimes such as that shown in Fig. 1. While the Gaussian approximation is conservative, it is sub-optimal compared to incorporating the full probability density function of the arrival time, and while in principle the full posterior distribution from this kind of analysis could be used in subsequent pulsar timing analysis, such an approach will always be an approximation to simply performing the analysis in the profile domain.

In this paper we present a solution to these problems by extending the profile domain pulsar timing framework (Lentati et al. 2016a and references therein, henceforth L16) to incorporate scattering variations. This means that the parameters that describe the scattering as a function of time can be estimated simultaneously with all other parameters in the model. This includes DM variations, the pulsar timing model, timing noise, pulse jitter, and also any models for the profile, profile evolution as a function of frequency, and pulse shape instability.

In order to incorporate scattering into the existing framework we perform the full analysis in the Fourier domain, as opposed to the time domain as in L16. In Section 2 we describe this modified framework, before applying it to simulated, and then real data in Sections 3-4, before ending with some concluding remarks in Section 5. All the analysis presented subsequently is performed using a new, publicly available GPU accelerated code111https://github.com/LindleyLentati/TempoNest2, and we have included the simulated data sets used in Section 3 in the repository with instructions for how to replicate the results presented in this paper.

2 A Fourier-domain model

In the following section we will build on the methods presented in L16. In that previous work, the profile-domain analysis was performed in the time-domain. Here, we will define our models in the Fourier-domain which allows us to optimize the analysis in two ways. Firstly, the Fourier representation makes it trivial to truncate our model for the profile at a particular harmonic, decreasing the size of the linear-algebra operations required to evaluate it. Secondly, when including scattering in the analysis, we only need to multiply the pulse broadening function with our profile model in the Fourier domain, as opposed to computing a convolution in the time domain, and the gradients can likewise be obtained at less computational expense for all parameters.

For a full description of the general profile domain framework we refer the reader to L16. Below we will give details of how the methodology has been changed to allow analysis in the Fourier domain, and of our implementation of models for scattering.

2.1 The Profile Model

As in L16 we construct our profile model using the shapelet basis (Refregier, 2003), for which the Fourier-transform can be obtained analytically. In the time domain, shapelets are described by a position $t$, a scale factor $\Lambda$, and a set of $n_{\mathrm{max}}$ amplitude parameters, with which we can construct the set of basis functions:

[TABLE]

with $\Phi_{n}(t)$ given by:

[TABLE]

where $H_{n}$ is the $n$th Hermite polynomial. In the Fourier domain, the dimensionless basis functions $\Phi_{n}(t)$ become:

[TABLE]

and the basis functions $B_{n}(t;\Lambda)$ become:

[TABLE]

for Fourier frequencies $f$. We can then write our profile model in the Fourier domain, $\tilde{s}(f,\mathbf{\zeta},\Lambda)$, as the sum:

[TABLE]

where $\zeta_{n}$ are the shapelet amplitudes, and $n_{\mathrm{max}}$ is the number of shapelet basis vectors included in the model.

As in L16 we have explicitly written the shapelet amplitudes as a function of the observing frequency $\nu$, and use a general polynomial expansion of the shapelet amplitudes with frequency in order to model any potential smooth profile evolution. This model is defined such that for the $p$ terms in the polynomial we can write the $n$th shapelet amplitude $\zeta_{n}(\nu)$ as:

[TABLE]

where $\nu_{c}$ is an arbitrary reference frequency, and $\zeta_{n,k}$ is the amplitude parameter for the $k$th polynomial of the $n$th term in the shapelet model.

Lastly, as in L16, we use the shapelet basis to describe the overall profile shape, and then scale this in amplitude for each epoch. We therefore use the zeroth-order term as a reference, taking $\zeta_{0}=1$, leaving only $n_{\mathrm{max}}-1$ free parameters $\zeta_{n}$ which are the amplitudes for the shapelet components with $n>0$. Written in this way Eq. 5 becomes:

[TABLE]

with $A$ the overall scaling factor for a particular epoch.

2.2 Shapelet interpolation

As in L16 we do not re-evaluate our shapelet model for every likelihood calculation, but adopt an interpolation scheme, where the shapelet basis vectors are precomputed on a grid from $t=0$ up to the duration of the longest phase bin in the data set. In principal, as we are performing our analysis in the Fourier domain we can rotate the shapelet model exactly by muliplying our complex shapelet model by a rotation vector:

[TABLE]

where $\delta\phi$ is the amount of phase by which to rotate the model. However, when performing the sampling with Hamiltonian Monte Carlo (HMC) methods, we must compute the gradients for all the shapelet amplitudes. As such, if we rotate the profile model to match the data, we must also rotate all the basis vectors for each profile in the data set. We find it is more efficient to pre-compute sub-bin shifts in the profile model, and then to rotate the data by an integer number of bins, again using Eq.  8 so that it aligns with the nearest interpolated set of basis functions. In this way we need only perform a single rotation per profile, which is a significant computational saving.

As in L16, we use an interpolation interval of 1 ns, chosen to be sufficiently small that no bias enters our analysis as a result of the interpolation process.

2.3 Truncating the Shapelet Model

In addition to the interpolation scheme described in Section 2.2, the Fourier representation makes it straightforward to truncate the shapelet model at a particular harmonic. This means that we need only include the subset of the harmonics in the profile model that contribute above some threshold in the analysis. While the interpolated shapelet model used in L16 significantly reduced the time taken to evaluate the mean profile compared to numerical evaluation, it still required a $\mathrm{N}_{\mathrm{b}}\times\mathrm{N}_{\mathrm{c}}$ matrix-vector product, where $N_{c}$ is the number of amplitudes in the shapelet model, and $N_{b}$ is the number of bins in the profile data. By only including Nh harmonics in the profile model, we can reduce this to a $2\mathrm{N}_{\mathrm{h}}\times\mathrm{N}_{\mathrm{c}}$ product, providing an immediate decrease in the computation time of the profile model by a factor of $\mathrm{N}_{\mathrm{b}}/2\mathrm{N}_{\mathrm{h}}$.

In Fig. 2 we compare the time domain (left-panel), and Fourier domain (right-panel) representations of the time averaged profile for PSR J1643$-$1224, using the 1200 MHz to 1500 MHz data that we analyse in Section 4. While the time-domain profile extends across $\sim 300$ phase bins, the Fourier representation contains over 99% of the total signal in only 40 harmonics, representing a significant reduction in the size of the matrix-vector multiplication required to evaluate the model. In Section 4 we include up to the 80th harmonic in our model. Beyond this the relative contribution of higher harmonics in the profile model is less than one part in $10^{10}$.

In principle, variations in the shape of the pulse profile may introduce fluctuations at harmonics higher than that required for the mean profile model. In the analysis presented here, such fluctuations will be absorbed by the white noise component of our model. In principle, however, one could compare models that incorporate these higher frequency terms and determine the optimal maximum harmonic to include in the analysis at the expense of increased computation time.

2.4 Scattering

The primary addition to the model used in our analysis in this work, as opposed to L16, is broadening of the pulse profile due to scattering as the pulse passes through the IISM. In the time-domain approach, pulse broadening would have to be included as a convolution of the pulse profile with a pulse broadening function (PBF). Typically this is performed by performing a Fast Fourier Transform (FFT) of the model and the PBF, multiplying them in the Fourier-domain, and then transforming back to the time domain. When calculating the gradient of the likelihood for each parameter in the model, these too require Fourier-transforms to be performed, and so the complexity of the likelihood calculation grows rapidly. In the Fourier domain we need only perform the multiplication of the Fourier transform of the PBF with our profile model, and the gradients can likewise be obtained trivially for all parameters.

In principal, any model for scattering that allows for the calculation of a gradient can be incorporated into our analysis. Even if the analytic Fourier transform of the PBF is not known, if the gradient can be computed in the time domain it can then be transformed via FFT for use in our analysis framework, although this will be less computationally efficient. In our analysis in this work we simply assume a thin screen model for the PBF (Williamson, 1972), using a single parameter for the time scale, $\tau$, which in the time domain will be given by:

[TABLE]

where $\nu$ is the observing frequency, and $H(t)$ is the Heavyside step function. Note that, as we always require the scattering timescale to be positive, we fit for $\log_{10}$ of the scattering timescale, $\bar{\tau}$. The analytic Fourier transform of Eq. 9 is then given by:

[TABLE]

with gradients with respect to $\bar{\tau}$ and $\alpha$:

[TABLE]

[TABLE]

where as before $i$ indicates a complex number and $\omega=2\pi f$. We can then multiply the Fourier representation of our shapelet model, $\mathbf{\tilde{s}}$, by Eq. 10 to get the scattered profile model.

In our analysis of both simulated data, and the PSR J1643 data set, we use a piecewise-constant $\tau$(t) model for the scattering. This therefore makes the same assumptions made when modelling DM variations using the DMX parameterisation (Demorest et al., 2013). In principal, one could also model the power spectrum of the scattering, under the assumption that it is a smooth, wide-sense stationary process in the same manner as the smooth model for DM variations in L16. Deviations from this simple model, such as those predicted from simulations (Coles et al., 2010) could then be encapsulated as additional, non-stationary features in the model, in the same manner as the extreme scattering events observed in PSR J1713+0747 (e.g. Desvignes et al. 2016; Lentati et al. 2016b). However, for this demonstration of the profile domain approach, we do not yet consider such a model.

2.5 Evaluating the profile domain model

The remainder of the timing framework described in L16 requires only minimal changes to operate on the Fourier representation of the profile data. As in L16, we consider the data in terms of a set of $N_{\mathrm{e}}$ epochs. Each epoch $i$ then has $N_{\mathrm{c},i}$ channels, such that the total number of profiles $N_{\mathrm{p}}=\sum_{i=1}^{N_{\mathrm{e}}}N_{\mathrm{c},i}$. The profile in the $j$th channel of the $i$th epoch then consists of a set of $N_{i,j}$ values representing the amplitude of the profile as measured at a set of times $\mathbf{t_{i,j}}$ which we denote $\mathbf{d}_{i,j}$. The Fourier representation of the profile data, $\mathbf{\tilde{d}}_{i,j}$ will therefore be a set of $N_{i,j}/2+1$ complex values. While in L16 a baseline offset for each profile was required as part of the model, in the Fourier representation we simply discard both the DC offset term and the nyquist term, so that our final data vector for each profile is of length $N_{i,j}/2-1$.

All further parameters, such as those describing pulse jitter or shape variation, can then be incorporated into our Fourier-domain shapelet model $\tilde{s}$ in exactly the same manner as in L16. We can then write the final likelihood that the data is described by the parameters in our model, which we collectively refer to as $\theta$, as:

[TABLE]

where $\textbf{{N}}_{\mathbf{{i,j}}}$ is the white noise covariance matrix for the Fourier domain profile corresponding to the $j$th channel in the $i$th epoch, with elements $(N_{i,j})_{mn}=\sigma_{i,j}\delta_{mn}$, with $\sigma_{i,j}$ the root-mean-square (RMS) deviation of the uncorrelated radiometer noise in the profile. When performing the analysis with the Guided Hamiltonian Sampler (GHS), $\sigma_{i,j}$ is a free parameter in our analysis for every profile.

3 Application to Simulated Data

In order to test the efficacy of the analysis method described in the preceding sections, we first apply it to a simulated data set. In particular we simulate 100 observational epochs, each assuming 256 MHz of bandwidth split into 8 channels, with a central frequency of 1369 MHz and covering a total time span of 4.5 years. Each epoch has an integrated S/N across the observing band of $\sim 500$. Both the MJD and the bandwidth of each simulated epoch are taken from the 1400 MHz data used in Section 4. We use a non-evolving, Gaussian profile, and a timing model consistent with that given in Table 1 for PSR J1643$-$1224.

In addition to this timing model, each epoch is subject to a random change in the DM with a standard deviation for the variations of $0.007$cm*-3*pc, and changes in the scattering of $0.81\pm 0.16$ ms at a reference frequency of 1.4 GHz. The scattering is simulated using the thin screen model described in Section 2.4, with the frequency scaling parameter, $\alpha$, equal to four. We then perform a simultaneous analysis of all epochs, fitting for:

• i)

The DM and scattering time scale at each epoch,

• ii)

the frequency scaling of the scattering, $\alpha$,

• iii)

the timing model parameters, and

• iv)

the pulse amplitude and instrumental noise level in each profile.

In Fig. 3 we show both the injected values (red crosses) and mean parameter estimates with 1 $\sigma$ confidence intervals obtained from this analysis (black points with error bars) for the scattering (top-left panel) and DM variations (top-right panel). We find the parameter estimates are consistent in both cases, and that the histogram of the residuals after subtracting the maximum likelihood values are consistent with a Gaussian distribution with unit variance (middle panels, red and black lines respectively). The Pearson product-moment correlation coefficient between the simulated and recovered values is $0.999\pm 0.005$. In Fig. 3 (bottom-left panel) we then show the one dimensional marginalised posterior distribution for the frequency dependence of the scattering model, $\alpha$. We indicate the simulated value with a vertical red line, and find our inferred values are consistent with the simulation.

We also use this simulated data set to check the impact of not modeling the scattering variations on the parameter estimates for the DM variations. In Fig. 3 (bottom-right panel) we plot the difference between the measured DM variations from the simulated data set when including scattering variations in the model, or when assuming only a mean scattering timescale (red points). Not modeling the variations in the scattering leads to significant bias in the measured DM variations. Note the different scales on the y-axis of the bottom-right, and top-right panel, the bias introduced in this case is a factor two larger than the DM variations themselves. We compare this bias with the structure present in the scattering variations (a shifted and scaled version of the scattering variations is overplotted in black). The Pearson product-moment correlation coefficient between the two sets of measurements is $0.92\pm 0.04$, indicating significant correlation between the unmodeled scattering, and the resultant bias in the DM variations.

4 Application to Real Data

We perform our analysis using observations of PSR J1643$-$1224 made with the 64-m Parkes radio telescope. We choose to analyse this particular MSP because it is known to display significant scattering variability, which has been attributed to its location behind an HII region (Gaustad et al., 2001; Villamariz & Herrero, 2005). As such we stress that the results presented will likely not be typical for a standard MSP. However, the benefit for more typical systems will naturally depend on combinations of factors such as the observing frequency of the data set, and the brightness of the pulsar at low observing frequencies. Further application of the techniques we have discussed in the preceding sections to a broader, and more typical population of pulsars will occur in subsequent work. Data were collected using two receiver packages, a co-axial system at 10 cm and 40 cm, and the centre pixel of a multi-beam 20 cm system. Data at 3100 MHz (co-axial, hereafter ‘10 cm’) and 1369 MHz (multi-beam, hereafter ‘20 cm’) were recorded using a digital polyphase filterbank (PDFB4) with a typical resolution of 1024 channels and 1024 bins and respective bandwidths of 1024 MHz and 256 MHz. Data from the lower co-axial band, centred at 732 MHz (hereafter ‘40 cm’), were recorded over a 64 MHz bandwidth with a similar polyphase system, PDFB3. More details about the observing systems and data reduction process are given in Manchester et al. (2013). For all three systems we average the profile data into eight channels per band.

As there will be discrete time offsets (known as jumps) between the different observing systems we include these as free parameters in our analysis together with the rest of the timing model parameters given in Table 1. In all the analysis that follows we include a model for the profile and smooth profile evolution as a function of frequency using a quadratic polynomial in the shapelet amplitudes, which we found was sufficient to describe the observed variation. In addition, we then include DM variations using the DMX parameterisation (Demorest et al., 2013) in which the variations in DM are modelled as a piecewise constant function, and we use an epoch length of 30 days, chosen so that the median number of 40 cm observation per epoch is one. Below we summarise three scenarios used to investigate the impact scattering variations have on the parameter estimates and uncertainties for the DM variations and timing model parameters.

• (S1)

We include only the 10 cm and 20 cm data. No scattering variations are included in the model.

• (S2)

As (S1) however we also include the 40 cm data.

• (S3)

As (S2), however we include scattering variations with the model described in Section 2.4, with the same cadence as the DMX parameters.

As in our analysis of simulated data in Section 3, we also attempt to recover the scaling of the scattering timescale with frequency. However, we find that the relatively narrow 40 cm band, combined with the fact that the scattering variations are already sub-dominant to the DM variations in the 20 cm band means we are unable to constrain the quantity in our analysis. We find that the only quantity that varies significantly when changing the scaling factor across a range of 3-5 is the mean scattering timescale, however the variations in scattering, and their impact on the other parameter estimates are consistent throughout. We have therefore set $\alpha=4$ in all the results discussed below, and will refer only to the magnitude of variations in the scattering timescale, rather than the absolute value.

4.1 DM Variations

In Fig. 4 (top panel) we compare parameter estimates for the DM variations in our PSR J1643$-$1224 data set for scenarios S1-S3 (red, blue and black points respectively).

Excess time-correlated noise has been observed in this pulsar previously at low observing frequencies ($<$ 1 GHz) (e.g., Keith et al. 2013; Lentati et al. 2016b henceforth L16a). In particular, L16a analysed the first International Pulsar Timing Array data set for PSR J1643$-$1224 and found that this ‘band-noise’ had a steep dependence on observing frequency, such that above 1 GHz the timing fluctuations were consistent with DM variations only. We can therefore expect that in scenario (S1), we will also be dominated by DM variations, rather than scattering variations. We find that the changes in the DM are consistent with the smooth variations observed in L16a when including band-noise terms.

When including the 40 cm data in our analysis in scenario (S2), significant additional structure can be seen in the DM that is inconsistent with the variations observed from scenario (S1). This structure has been reported previously in the literature as being the result of yearly variations in the DM (e.g., Arzoumanian et al. 2015; Jones et al. 2016)

In scenario (S3) the additional structure seen in scenario (S2) is no longer present, and we find the parameter estimates are consistent with those obtained from scenario (S1). That the parameter estimates for the DM variations are inconsistent between Scenarios (S1) and (S2) is highly suggestive that unmodelled effects present in the 40 cm data are introducing bias into our analysis, and are subsequently mitigated when including scattering variations into the model.

4.2 Scattering Variations

In the middle panel of Fig. 4 we show the parameter estimates for the variations in the scattering timescale, $\tau$, from scenario (S3). In the bottom panel of Fig. 4 we compare this structure with the additional structure seen in the DM variations in scenario (S2) compared to scenario (S1). In particular, we overplot the difference in the DM variations from scenarios (S1) and (S2) (red points) with a shifted and scaled version of the scattering variations (black points). The two clearly mirror one another, implying that the additional DM variations observed when adding low frequency data in scenario (S2) are the result of unmodelled scattering variations, consistent with the result of not modeling the scattering variations observed in the analysis of simulated data in Section 3.

It is clear, then, that simply adding low-frequency data without considering scattering can significantly bias estimates of the DM variations. With LOFAR (Stappers et al., 2011) already performing pulsar timing at frequencies of 100 MHz to 200 MHz, and upcoming telescopes such as the Square Kilometre Array providing unprecedented sensitivity in the frequency range of 500 MHz to 800 MHz, this bias introduced into the analysis when failing to appropriately model scattering variations will become extremely important, especially in the context of gravitational wave astronomy using a pulsar timing array. Looking forward this will become increasingly true even if we only consider more ’typical’ pulsars compared to PSR J1643$-$1224. This is both because in the SKA era we can expect to add many distant pulsars to PTAs, for which the presence of detectable variations in scattering will be increasingly probable, and also simply because by observing for longer the impact of such effects will grow.

4.3 Timing Bias

Existing methods for incorporating scattering variations into pulsar timing have suggested using maximum-likelihood estimates of scattering delays and ‘correcting’ the measured ToA by fixing that delay in the subsequent timing analysis. This can be done either by obtaining estimates of the delay from measurements of scintillation, as in Levin et al. (2016), or using cyclic spectroscopy (e.g., Palliyaguru et al. 2015). With the profile domain framework we can directly compare the impact of fixing the scattering timescales in the analysis at their maximum-likelihood values, which is equivalent to simply fixing the time delay induced by the scattering when performing the analysis using ToAs. We find that fixing the scattering timescales results in significant bias in the uncertainties for the DM variations of a factor of $\sim 1.8$, with some epochs being much larger when less multi-frequency information was available.

In Table 1 we list the timing model parameter estimates for all three scenarios, and in Fig. 5 we then compare the parameter estimates between scenarios (S1) and (S2), and scenarios (S3) and (S2). The bias in the parameter estimates introduced as a result of not modeling the scattering variations is substantial in several of the astrometric parameters. From the middle panel of Fig. 4 there is a clear annual term in the scattering signal which we would expect to correlate strongly with both position and proper motion. This is indeed the case, with both position in RA, and proper motion in RA suffering $\sim 4~{}\sigma$ changes in the parameter estimates. For example, the proper motion in RA changes from $6.0\pm 0.1$ to $5.50\pm 0.06$ when first adding in the 40 cm data without a model for scattering variations, and then changes to $5.85\pm 0.08$ when including scattering variations in the model.

A natural question to ask is whether it is worth including low-frequency pulsar timing data into the analysis, if one is interested only in the timing model parameters and not the behavior of the ISM. From Table 1 we can see that the uncertainties on the parameters from Scenario (3) are still improved compared to Scenario (1), in which only the 10 cm and 20 cm data were included in the analysis. Therefore, if the scattering variations are modeled appropriately such that no bias is introduced in the parameter estimates, the addition of the 40 cm data does improve the precision of the results.

However, this implies we must also ask whether the model used for scattering is sufficient. If the thin screen assumption is not valid, then bias will still enter into our parameter estimates. In Fig. 6 we show the timing residuals after subtracting the timing model given in Table 1 obtained from our analysis of scenario (S3). In the middle panel we then show the timing residuals obtained from the maximum likelihood model profile data that includes variation in the pulse width due to scattering. In both cases significant additional structure in the 40 cm data (magenta points) can clearly be seen as a result of the scattering. In the bottom panel of Fig. 6 we then show the difference between the first two panels. The residuals, after subtracting out delays that result from scattering variations, are consistent with white noise. We stress that this process is performed purely for visualization of the profile domain model. As discussed, simply subtracting the delay induced by scattering is not an appropriate means of correcting for scattering variations.

5 Conclusions

Over the last several years different groups worldwide have sought to improve the precision with which pulsar timing can be performed by incorporating data across wide band widths, extending to low radio frequencies ($<$ 1 GHz). The principal aim was to better correct for DM variations, which are known to dominate the noise budget for the majority of millisecond pulsars being used as part of an ongoing effort to detect gravitational waves in the nHz band.

As better models for these DM variations were developed, excess noise in this low frequency data has been observed in an increasing number of pulsars. In Lentati et al. (2016b), four pulsars in the International Pulsar Timing Array data release were found to have ‘band noise’, including PSR J1643$-$1224. Even the most stable millisecond pulsar known to date, PSR J1909$-$3744, has been observed to suffer from band noise in the low frequency data, which, unless modeled appropriately, decreased the sensitivity of the data set to gravitational waves by a factor of two (Shannon et al., 2015).

These band-noise models described the excess noise as time-correlated shifts in the arrival time of the pulses. Ultimately, however, if the origin of this noise is the result of time-variable broadening of the pulse profile due to scattering in the IISM, describing it in terms of shifts in the ToAs will be sub-optimal compared to appropriately modeling the changes in the pulse shape.

We have described a statistically robust approach to incorporating models for scattering into pulsar timing analysis by extending the existing profile-domain timing framework. This approach makes it possible to simultaneously estimate temporal variations in both the DM and scattering, together with a model for the pulse profile that includes smooth evolution as a function of frequency, the pulsars timing model, and models for pulse jitter.

We have shown that unless both the DM variations and scattering are sampled together in the timing analysis one can significantly underestimate the uncertainties in both the DM, and the arrival time of the pulse, leading to bias in the timing parameters. This means that fixing the scattering timescales and simply ‘correcting’ the ToAs, as has been suggested as a means of handling scattering variations in the literature, is not a viable approach to performing a robust analysis of low frequency pulsar timing data. We also showed, however, that the true probability density function (PDF) of the arrival time of the pulse can be highly non-Gaussian when scattering is only marginally detected. In this case, the typical procedure for forming a ToA, which assumes a Gaussian approximation to the true PDF, is inappropriate. However, by performing a profile domain analysis the non-Gaussianity of the PDF is automaticlaly incorporated into the final parameter estimates.

We have applied our method using a new, publicly available, GPU accelerated code both to simulations, and observations of the millisecond pulsar PSR J1643$-$1224 that cover a frequency range from 700 MHz up to 3800 MHz over a period of 4.7 yr. This particular MSP was chosen because it is known to display significant scattering variability compared to typical MSPs. The techniques we have introduced will be applied to a broader, and more typical population of pulsars in subsequent work. When including the low frequency data in the PSR J1643$-$1224 analysis the parameter estimates for both the DM variations and timing model are significantly biased by the variations in the pulse width due to scattering. Further, simply fixing the scattering timescales to maximum likelihood values - equivalent to fixing the time offset when making ToAs - results in almost a factor two underestimation in the uncertainties of the DM variations compared to the simultaneous analysis.

If low-frequency observations from telescopes such as LOFAR, or the Square Kilometre Array are to be incorporated into pulsar timing analysis then the bias that results from failing to appropriately model scattering variations will become extremely important. The increased sensitivity afforded by these instruments will likely require more complex models than simple thin-screen approximations for scattering. Any differences between the true pulse-broadening function and the models used will manifest as additional band-dependent noise, down-weighting the low-frequency data. Given finite observing time, one must ask how much of that time to spend observing at low-frequencies. If we are able to completely model the variations in pulse shape, then low-frequency observations will naturally be of great benefit. The actual improvement, however, will vary from pulsar to pulsar depending on the strength and complexity of the scattering variations.

For pulsars that are bright at high frequencies ($\ga 3$GHz) the optimal approach may well be to spend all available observing time in a regime where IISM effects are negligible. Such an approach has already resulted in the most sensitive limit on an isotropic gravitational wave background using a pulsar timing array to date (Shannon et al., 2015). However, pulsars are known to have steep spectral indices ($S_{\nu}\propto\nu^{-1.8}$ on average, Maron et al. 2000), and so in many cases observations at lower frequencies are necessary. In such instances our ability to propagate our uncertainties in the scattering variability through to all other parameters will be critical, especially in the context of gravitational wave astronomy. This is all handled automatically as part of a profile domain analysis, which will inevitably form the basis of modern, broad-band pulsar timing analysis moving forward.

6 Acknowledgements

The Parkes radio telescope is part of the Australia Telescope National Facility which is funded by the Commonwealth of Australia for operation as a National Facility managed by Commonwealth Science and Industrial Research Organization (CSIRO). We thank all of the observers, engineers, and Parkes observatory staff members who have assisted with the observations reported in this paper.

LL was supported by a Junior Research Fellowship at Trinity Hall College, Cambridge University. SO was supported by the Alexander von Humboldt Foundation.

We thank Stephen Taylor for his assistance in building the python software used in this paper.