Modelling annual and orbital variations in the scintillation of the relativistic binary PSR J1141$-$6545

TL;DR

This study uses long-term scintillation observations of pulsar PSR J1141$-$6545 to measure orbital, spatial, and interstellar medium parameters, improving pulsar astrometry and understanding of interstellar turbulence.

Contribution

It models annual and orbital scintillation variations to derive pulsar system parameters and interstellar medium characteristics, resolving previous ambiguities and providing new measurements.

Findings

Measured orbital inclination and longitude of ascending node.

Estimated pulsar distance and proper motion.

Characterized interstellar scattering properties and anisotropy.

Abstract

We have observed the relativistic binary pulsar PSR J1141$-$6545 over a period of $\sim$6 years using the Parkes 64 m radio telescope, with a focus on modelling the diffractive intensity scintillations to improve the accuracy of the astrometric timing model. The long-term scintillation, which shows orbital and annual variations, allows us to measure parameters that are difficult to measure with pulsar timing alone. These include: the orbital inclination $i$; the longitude of the ascending node $\Omega$; and the pulsar system transverse velocity. We use the annual variations to resolve the previous ambiguity in the sense of the inclination angle. Using the correct sense, and a prior probability distribution given by a constraint from pulsar timing ($i=73\pm3^\circ$), we find $\Omega=24.8\pm1.8^\circ$ and we estimate the pulsar distance to be $D=10^{+4}_{-3}$ kpc. This then gives us an estimate of this pulsar's proper motion of $\mu_{\alpha}\cos{\delta}=2.9\pm1.0$ mas yr$^{-1}$ in right ascension and $\mu_{\delta}=1.8\pm0.6$ mas yr$^{-1}$ in declination. Finally, we obtain measurements of the spatial structure of the interstellar electron density fluctuations, including: the spatial scale and anisotropy of the diffraction pattern; the distribution of scattering material along the line of sight; and spatial variation in the strength of turbulence from epoch to epoch. We find that the scattering is dominated by a thin screen at a distance of $(0.724\pm0.008)D$, with an anisotropy axial ratio $A_{\rm r} = 2.14\pm0.11$.

Tables

Table 1. Table 1: Model comparison with Bayes factors ( log ⁡ K 𝐾 \log{K} , Equation 8 ) relative to the standard model of an isotropic thin screen.

Geometry	$\log K$
Uniform $\cos i$ prior
Isotropic screen	0
Isotropic uniform	$-$292.0
Anisotropic screen	82.9
Anisotropic uniform	54.6
Gaussian $i$ prior: $73\pm 3^{\circ }$
Isotropic screen	$-$4.6
Isotropic uniform	$-$426.5
Anisotropic screen	82.1
Anisotropic uniform	$-$64.0

Table 2. Table 2: Measured and derived parameters for the long-term scintillation model of an anisotropic thin screen, using two different priors for the inclination angle i 𝑖 i . The Gaussian prior had a mean and standard deviation of 73 ∘ superscript 73 73^{\circ} and 3 ∘ superscript 3 3^{\circ} respectively, corresponding to the 2 σ 𝜎 \sigma confidence region of (Bhat et al., 2008 ) . The value given in brackets are the uncertainty on the last quoted decimal place. The uncertainties for our measurements of v μ , ∥ v_{\mu,\parallel} and v μ , ⟂ subscript 𝑣 𝜇 perpendicular-to v_{\mu,\perp} have been inflated with 10 10 10  km s -1 added in quadrature to account for possible contributions from an unknown V IISM subscript 𝑉 IISM V_{\rm{IISM}} of approximately this magnitude.

Measured parameters
	Inclination angle prior
	Uniform $\cos i$	Gaussian $73\pm 3^{\circ }$
$\kappa$	$0.432\left(17\right)$	$0.496\left(14\right)$
$s$	$0.299\left(9\right)$	$0.276\left(8\right)$
$\Omega$ (∘)	$27.1\left(18\right)$	$24.8\left(18\right)$
$i$ (∘)	$61.3\left(20\right)$	$68.6\left(1.1\right)$
$v_{\mu ,\parallel }$ (km s-1)	$21\left(11\right)$	$19\left(11\right)$
$v_{\mu ,⟂}$ (km s-1)	$242\left(23\right)$	$173\left(15\right)$
$\psi$ (∘)	$87.8\left(3\right)$	$86.8\left(4\right)$
$R$	$0.77\left(3\right)$	$0.64\left(3\right)$
Derived parameters
$A_r$	$2.80\left(18\right)$	$2.14\left(11\right)$
$M_p$ ($M_{\odot }$)	$1.18\left(2\right)$	$1.242\left(8\right)$
$M_c$ ($M_{\odot }$)	$1.11\left(3\right)$	$1.047\left(8\right)$
$D$ (kpc)	${12}_{-3}^{+4}$	${10}_{-3}^{+4}$
${\mu }_{\alpha }\cos \delta$ (mas yr-1)	$3.3\left(10\right)$	$2.9\left(10\right)$
${\mu }_{\delta }$ (mas yr-1)	$2.1\left(6\right)$	$1.8\left(6\right)$

Equations

$$C(\tau ,0)$$

$$C(0,0)$$

$$C(0,\delta \nu )$$

$$C(0,0)$$

$$V_{\mathrm{ISS}}=A_{\mathrm{ISS}}\frac{\sqrt{D\Delta {\nu }_{\mathrm{d}}}}{f{\tau }_{\mathrm{d}}},$$

$$V_{\mathrm{eff}}(s)=(1-s)(V_{\mathrm{p}}+V_{\mu })+s{V}_{\mathrm{E}}-V_{\mathrm{IISM}}(s),$$

$$V_{\mathrm{model}}={\kappa }_u\sqrt{v_{\mathrm{p},\parallel }^2+v_{\mathrm{p},\perp }^2},$$

$$\frac{1}{{\tau }_{\mathrm{d}}(\varphi {)}^2}=K_0+K_S\sin \varphi +K_C\cos \varphi +K_{S2}\sin (2\varphi )+K_{C2}\cos (2\varphi ),$$

$$P(\theta |d)=\frac{P(d|\theta )P(\theta )}{Z},$$

$$\log K=\log Z_h-\log Z_0,$$

$$V_{\mathrm{ISS}}={\left[\frac{\underset{0}{\overset{1}{\int }}\eta (x)|V_{\mathrm{eff}}(x){|}^{\alpha }dx}{\underset{0}{\overset{1}{\int }}\eta (x)x^{\alpha }dx}\right]}^{1/\alpha },$$

$$E(V_{\mathrm{ISS}})=V_{\mathrm{los}}=V_{\mathrm{eff}}/s,$$

$$v_{\mathrm{p},\parallel }$$

$$v_{\mathrm{p},\perp }$$

$$\theta =2\arctan \left[\sqrt{\frac{1+e}{1-e}}\tan \frac{E}{2}\right].$$

$$\left(\begin{array}{c}{v}_{\mathrm{E},\parallel }\\ v_{\mathrm{E},\perp }\end{array}\right)=\left(\begin{array}{cc}\sin \Omega & \cos \Omega \\ \cos \Omega & -\sin \Omega \end{array}\right)\left(\begin{array}{c}{v}_{\mathrm{E},\alpha }\\ v_{\mathrm{E},\delta }\end{array}\right)$$

$$v_{\parallel }(s)$$

$$v_{\perp }(s)$$

$$V_{\mathrm{eff}}(s)=\sqrt{v_{\parallel }(s{)}^2+v_{\perp }(s{)}^2}.$$

$$a$$

$$b$$

$$c$$

$$V_{\mathrm{eff}}(s)=\kappa \sqrt{a{v}_{\parallel }(s{)}^2+b{v}_{\perp }(s{)}^2+c{v}_{\parallel }(s)v_{\perp }(s)},$$

$$K_0=$$

$$+a(v_{\mathrm{C},\parallel }-V_{0}e\sin \omega {)}^2$$

$$+b(v_{\mathrm{C},\perp }+V_{0}e\cos \omega \cos i{)}^2$$

$$+c(v_{\mathrm{C},\parallel }-V_{0}e\sin \omega )(v_{\mathrm{C},\perp }+V_{0}e\cos \omega \cos i)]/s_{\mathrm{d}}^2$$

$$K_S=$$

$$+c(v_{\mathrm{C},\perp }+V_{0}e\cos i\cos \omega )]/s_{\mathrm{d}}^2$$

$$K_C=$$

$$+2b(v_{\mathrm{C},\perp }+V_{0}e\cos i\cos \omega )]/s_{\mathrm{d}}^2$$

$$K_{S2}=$$

$$K_{C2}=$$

Full text

††pubyear: 2018††pagerange: Modelling annual and orbital variations in the scintillation of the relativistic binary PSR J1141$-$6545–Appendix B: Reproducing our results

1 Introduction

PSR J1141$-$6545 is a relativistic binary pulsar in a $\sim$4.7 hour eccentric orbit with a young, white dwarf companion (discovered by Kaspi et al., 2000). This system has proved to be a unique and valuable laboratory for testing general relativity (GR) in such asymmetrical mass systems (Bhat et al., 2008; Manchester et al., 2010), and its runaway velocity is interesting for investigating their formation and binary evolution mechanisms (e.g. Tauris & Sennels, 2000; Davies et al., 2002; Church et al., 2006). The short orbital period and large centripetal acceleration of this pulsar makes it ideal for measuring and modelling changes in its scintillation timescale (Ord et al., 2002b) for the purpose of measuring astrometric parameters. Parameters including the inclination angle $i$ (and its sense), the longitude of ascending node $\Omega$, and the transverse velocity are useful for tests of GR and formation scenarios. We are able to uniquely determine these parameters and estimate the pulsar distance and proper motion by modelling the long-term scintillation of PSR J1141$-$6545 and accounting for anisotropy in the scattering. This is only possible by exploiting the annual variations observed with long-term monitoring of the pulsar.

Intensity variations caused by interstellar scintillation are seen in all observations of radio pulsars at centimetre to metre wavelengths. They are caused by transverse fluctuations in the electron density of the turbulent ionised interstellar medium (IISM). The primary mechanism is interference by waves scattered through diffraction (Rickett, 1969). The diffractive scintillations are modulated by refractive scintillations on larger spatial scales (Rickett et al., 1984; Romani et al., 1986).

The observed scintillations, which appear as variations in the source flux with time and observing frequency, are caused by the spatial diffraction pattern drifting across the line-of-sight. Thus the time scale of the observed scintillations $\tau_{\rm d}$ (typically of order minutes) is inversely proportional to this drift velocity $V_{\rm los}$, which is the velocity at which the line-of-sight crosses the IISM (Section 3). The diffraction pattern is frequency-dependent and becomes decorrelated over a bandwidth $\Delta\nu_{\rm d}$. The angular scattering broadens each pulse into a quasi-exponential pulse with timescale $\tau_{\rm s}$, which is related to the bandwidth by $2\pi\Delta\nu_{\rm d}\tau_{\rm s}\approx 1$ (Rickett, 1977). Strong diffractive scintillations have a narrow bandwidth (typically of order MHz), which is a useful measure of the strength of scattering and can be used to estimate the diffractive spatial scale $s_{\rm d}$ (Cordes & Rickett, 1998). Detailed overviews of pulsar scintillations are given by Rickett (1990) and Narayan (1992).

Measurements of a dynamic spectrum of intensity scintillation with time and frequency can therefore provide information on the spatial structure of the IISM, the transverse velocity of the pulsar, and the strength of scattering. Although scattering occurs throughout the line-of-sight, it is often dominated by one, or a few, local regions of more intense scattering.

For solitary pulsars $V_{\rm los}$ is often dominated by the pulsar proper motion and is relatively simple to model, depending only on the distance from the scattering region to the Earth (Lyne & Smith, 1982). If the pulsar has a binary companion, orbital dynamics can also be studied from the transverse orbital motion. This was first used by Lyne (1984) to measure the orbital inclination angle of PSR B0655$+$64 for the first time. Then Ord et al. (2002b) analysed two consecutive orbits of PSR J1141$-$6545 and found that the time scale varied smoothly over the orbit, but since they had to assume that the scattering was isotropic the parameter estimates likely suffered a bias. The technique was extended to deal with anisotropic scattering to analyse scintillations of the double pulsar PSR J0737$-$3039A by Coles et al. (2005). It was further extended by Rickett et al. (2014) to analyze multiple observations of PSR J0737$-$3039A over several years, which requires including the variation of the Earth’s velocity into the model. We have been provided access to the latter data and used it to calibrate our analysis against that of Rickett et al. (2014).

For this work we use the dynamic spectrum as our basic observable and model the diffractive scintillations of PSR J1141$-$6545. The physical parameters measured in this way are important for testing a variety of theories such as GR and other boost-invariant theories of gravity (Damour & Taylor, 1992). However, poor timing precision currently limits the measurement of the Shapiro delay, which would otherwise provide a measurement of the inclination angle and companion mass to further constrain the tests of GR (Bhat et al., 2008). An independent measurement of the inclination angle from scintillation will provide complementary tests of GR. In addition, the distance to the pulsar is poorly constrained at present, so the contamination from kinematic effects (in this case the Shklovskii effect; Shklovskii, 1970) in the measured relativistic orbital parameters is unknown. A lower-bound distance estimate of 3.7 kpc was found by Ord et al. (2002a) using HI absorption spectra and this was used in the scintillation modelling Ord et al. (2002b). However Verbiest et al. (2012) showed that distances derived in this way may be overestimated because of a luminosity bias, and presented a revised distance estimate of $3\pm 2$ kpc (we take this as an initial value for our models). At this large distance, the proper motion is too small to be measured with current timing precision. Further, the white dwarf companion was not identified in a recent targeted search for binary pulsar companions (Jennings et al., 2018) in the second Gaia data release (Gaia Collaboration et al., 2018), which could have provided independent estimates of distance and proper motion. For these reasons, both the measured transverse velocity (to estimate the Shklovskii effect) and inclination angle (to constrain the pulsar and companion masses) from scintillation modelling are particularly useful for improving the tests of GR with this system.

In this paper we present new short-term and long-term scintillation models for different scattering geometries, which are described in Section 3. We compare the evidence for each of these geometries using Bayesian methods in Section 4. Using the short-term models, in Section 5.1 we demonstrate that the Earth’s velocity is detectable from scintillation timescale modulation over a year, and that the relativistic advance of periastron can be measured from scintillation alone. Our best long-term model, presented in Section 5.2 is a thin scattering screen and includes anisotropy in the IISM. The annual variation in this model allows us to measure the orientation of the pulsar’s orbit in celestial coordinates, resolve the sense of the inclination angle, constrain the scattering anisotropy, determine the distance to the scattering region, and estimate the proper motion of the pulsar. We compare these results with previous measurements from scintillation and pulsar timing, provide a revised distance estimate to the pulsar in Section 6.1, and predict the contamination in the orbital period-derivative from the Shklovskii effect in Section 6.2.

2 Dataset

2.1 Observations and dynamic spectra

For this work, we use a selection of archival PSR J1141$-$6545 observations from the Parkes 64 m radio telescope, spanning $\sim$6 years from June 2009 to June 2015. The observations were part of the P361 observing project, which was a long-term campaign to improve the stringent tests of gravitational theories. The data were received with the central beam of the Parkes 20 cm multibeam receiver, and recorded either with a digital polyphase filterbank system with a 256 MHz bandwidth and 0.25 MHz channel width, or a coherent dedispersion machine with a 400 MHz bandwidth and 0.78 MHz channel width. We selected only the observations which spanned at least 142 minutes, so that each observation covered at least half of an orbit. This yielded 126 individual observations which we treated as 23 distinct “epochs" separated by at least 60 days.

A dynamic spectrum for each observation was produced using the data processing pipeline designed for the upcoming data release 2 (DR2) of the Parkes Pulsar Timing Array (PPTA; Manchester et al., 2013) project. In brief, observations of a pulsed noise diode that excites both X and Y polarisations in phase are performed before each of the observations described above, to allow correction of the complex gain. The noise diode is itself calibrated to absolute flux density using on- and off-source observations of the bright radio galaxy Hydra A. Polarisation calibration is done using the noise diode observations combined with regular observations of the highly polarised pulsar PSR J0437$-$4715. To compute dynamic spectra, we perform a least-squares fit of an analytic model of the pulse profile in total intensity (Stokes I) to the observed pulse profile for each sub-integration and frequency channel. This fit provides the amplitude and its uncertainty. Because the observed pulsar profiles are already absolutely calibrated, the amplitude measurement yields the pulsar flux density directly. In addition, using this analytic pulse profile simultaneously optimises the signal-to-noise ratio for both the pulse amplitude and time of arrival. The calibrations and measurements described above are performed with the psrchive (Hotan et al., 2004) package.

Figure 1 shows the dynamic spectra of four consecutive observations of PSR J1141$-$6545 on MJD 56391. The modulation of the scintillation timescale $\tau_{\rm d}$ caused by the orbital motion of the pulsar is visible by eye. We characterise the statistics of this spectrum by a two dimensional autocovariance function (ACF). Following convention (Cordes & Rickett, 1998), the half-width at half maximum of this ACF in frequency is the decorrelation bandwidth $\Delta\nu_{\rm d}$ and the half-width at $1/e$ in time is $\tau_{\rm d}$. The scale of spatial variations in the diffraction pattern $s_{\rm d}$ is related to the scintillation timescale by $\tau_{\rm d}=s_{\rm d}/V_{\rm los}$. The spatial scale is determined by the strength of scattering, which can vary with time but can be estimated from $\Delta\nu_{\rm d}$. If the spatial structure is anisotropic, $\tau_{\rm d}$ will depend on the angle between the semi-major axis of the anisotropy and the velocity vector, which our long-term data is sensitive to.

Fortunately, since $\Delta\nu_{\rm d}$ is a direct measure of the strength of scattering, this can be used to correct for temporal variations in the spatial scale caused by changes in the strength of scattering, provided that the anisotropy does not also change with time. A formalism for using $\Delta\nu_{\rm d}$ to correct for changes in $s_{\rm d}$ has been provided by Cordes & Rickett (1998). The near-constant $\Delta\nu_{\rm d}$ for the dynamic spectrum in Figure 1 shows that the strength of scattering does not change significantly over the orbital period. This is because the projected size of the orbit is smaller than the scattering disk. However we do observe small strength of scattering changes from epoch to epoch (discussed in Section 2.3) and we use the scheme of Cordes & Rickett (1998) to account for this.

We cut each observation into segments $<12$ minutes in length (on average they are $\sim 11$ minutes) and measure the scintillation timescales and bandwidths for each segment (as described in the following section). In this way we can measure the modulation of $\tau_{\rm d}$ across orbital phase for each observation. The decorrelation bandwidth $\Delta\nu_{\rm d}$ measurements and the grouping of observations into the 23 epochs are shown in Figure 2. The epoch to epoch variation represents real changes in the turbulence level of the IISM and is not related to the motion of the pulsar or the Earth. By comparison with the 20 millisecond pulsars monitored by the Parkes Pulsar Timing Array (Keith et al., 2013) this variation is quite modest.

The measurements in Figure 2 are from all available archival observations (with $t_{\rm obs}>24$ mins) during our selected observing span, not just the 126 observations from P361 that we use for modelling. This is because many of the P361 observations have channel bandwidths $B_{\rm c}=0.78$ MHz, which is greater than $\Delta\nu_{\rm d}$. For these observations the ACF measurement is not useful but we can use the rms flux over the spectrum to estimate the true $\Delta\nu_{\rm d}$. This method is described in Section 2.3. The estimates are calibrated by a scaling factor using observations in the same epoch for which $\Delta\nu_{\rm d}>B_{c}$.

2.2 Measurement of $\tau_{d}$ and $\Delta\nu_{d}$

We measure the scintillation parameters from $\sim$11 minute segments of the dynamic spectrum because the scintillation timescale varies rapidly for the relativistic orbit. We characterise the statistics of these segments using the estimated autocovariance function (ACF), $C(\tau,\delta\nu)$. To calculate this we remove the mean, then pad each segment with an equal length of zeroes in both dimensions, perform a 2-D FFT on the zero-padded segment, take the squared magnitude of the result, and perform an inverse 2-D FFT. We then perform a least squares fit of analytical models to $C(\tau,0)$ and $C(0,\delta\nu)$ to obtain $\tau_{\rm d}$ and $\Delta\nu_{\rm d}$ respectively. First we fit $C(\tau,0)$ with

[TABLE]

where $T_{\rm obs}$ is the length of the segment ($\sim$11 minutes), $\Lambda(\tau,T_{\rm obs})$ is the triangle function of length $\pm T_{\rm obs}$, and $W$ is the variance noise spike. This function is a slight modification to the previous standard of a Gaussian function (Cordes & Rickett, 1998), where the exponent of $5/3$ gives better fit to the shape of Kolmogorov scintillations (e.g. Coles et al., 2005, 2010; Rickett et al., 2014). This does not bias the estimate of $\tau_{\rm d}$, but slightly reduces its uncertainty. After obtaining $A$, $W$, and $\tau_{\rm d}$, we keep $A$ and $W$ constant because we typically have marginal resolution in frequency and it is easier to isolate the noise spike using the time cut, $C(\tau,0)$. We then fit $C(0,\delta\nu)$ with

[TABLE]

to obtain $\Delta\nu_{\rm d}$, where B is the receiver bandwidth.

This fitting provides straightforward estimators of $\tau_{\rm d}$ and $\Delta\nu_{\rm d}$, which are easily checked manually, but the least squares fit is not optimal because samples of the observed ACF are heavily correlated. We therefore repeat the fit using the same analytical models, but we perform the fit in the Fourier transform domain where we simply transform the autocovariance function and the models.

The 2-D FFT of $C(\tau,\delta\nu)$ is the power spectrum (or secondary spectrum) $P(f_{\rm dop},t_{\rm del})$, where the dimensions are the differential time delay $t_{\rm del}$ and the differential Doppler shift $f_{\rm dop}$ of the interfering waves. To obtain $\tau_{\rm d}$ from this data we first sum $P(f_{\rm dop},t_{\rm del})$ over the $t_{\rm del}$ dimension and divide by the number of samples $N_{\rm del}$ and then we fit the transform of Equation 1. As before we hold $A$ and $N$ fixed from this fit before obtaining $\Delta\nu_{\rm d}$ by summing $P(f_{\rm dop},t_{\rm del})$ over the $f_{\rm dop}$ dimension, dividing by the number of samples $N_{\rm dop}$, and fitting the resulting power spectrum with the transform of Equation 2.

The errors on these average power spectra are independent, but not equal. In fact they are proportional to the average power spectra itself, so we use a weighted least squares fit with the models providing the weights. This approach provides a second (but not independent) estimator for both $\tau_{\rm d}$ and $\Delta\nu_{\rm d}$ for which we believe the measurements and uncertainty estimates are more reliable. We check that the two methods agree within the uncertainty of the second method and review the dynamic spectra if they do not. For this work we use the $\tau_{\rm d}$ and $\Delta\nu_{\rm d}$ measurements from this Fourier-domain method.

2.3 The effects of inhomogeneity in the IISM

The primary physical mechanism underlying scintillation is angular scattering by small scale irregularities in electron density that is diffractive in nature. Observations that are sensitive to the spectral exponent of the power-law of density fluctuations have shown a Kolmogorov spectrum truncated at an inner scale (e.g. Armstrong et al., 1995; Rickett et al., 2009). Thus it has been assumed that the microstructure in the IISM is turbulent in origin and by default they have been assumed homogeneous. However more recent observations of phenomena such as extreme scattering events (ESEs) (e.g. Fiedler et al., 1987; Coles et al., 2015) suggest that the turbulence is often inhomogeneous or that inhomogeneous structures that dominate the scattering are often present at some place on the line of sight from the source to the observer.

Furthermore the spatial structure is now often found to be localised (along the line-of-sight) and anisotropic, for example through analysis of parabolic arcs in the two-dimensional fourier transform of the dynamic spectrum that were first discovered by Stinebring et al. (2001). This has been called the “secondary spectrum" or the “delay-Doppler distribution" because its axes are the differential Doppler shift and the differential time delay of the interfering waves that cause the intensity variations (e.g. Walker et al., 2004; Cordes et al., 2006; Brisken et al., 2010). The parabolic arcs occur when the scattering is dominated by a compact local region on the line of sight. We have analysed the secondary spectra for PSR J1141$-$6545 and found no clear parabolic arcs that could be used provide information on the scattering geometry, likely because of unfavourable sampling. Arcs will be displayed only when dynamic spectra are significantly oversampled in both time and frequency. Thus we have no a priori evidence from our observations for a particular scattering geometry to use in our models of the scintillation velocity.

The spatial scale of the diffractive scattering, $s_{\rm d}$ is defined as the transverse separation where incident waves have a 1 radian rms difference in phase. The width of the angular scattering is then $\theta_{\rm d}\approx 1/\left(ks_{\rm d}\right)$, for incident radiation with wavenumber $k$ (Rickett, 1990). The radiation received by the observer arrives from a scattering region of diameter $s_{\rm r}=\theta_{\rm d}D_{\rm e}$, where $D_{\rm e}$ is the distance to the scattering screen from the observer. The intensity will also show variations from refraction, with a spatial scale equal to the diameter of this scattering region (the refractive scale). We estimate the size of the scattering region $s_{\rm r}$ from $s_{\rm r}/r_{F}=\sqrt{\nu/\Delta\nu_{\rm d}}$, where $r_{F}=\sqrt{D_{\rm e}/k}$ is the “Fresnel scale," and find $s_{\rm r}\approx 10^{8}$ km (AU-scale) for scintillation bandwidth $\Delta\nu_{\rm d}=0.36$ MHz (mean from Figure 2) and approximate distance $D_{e}\sim D=3$ kpc (Verbiest et al., 2012). Because the projected size of the orbit is $5.6\times 10^{5}$ km at this distance, the scattering disk is much larger than the orbital diameter regardless of the location of the scattering screen or a significant error in the pulsar distance. So for a single binary orbit, the line of sight to PSR J1141$-$6545 does not travel outside of the scattering region. Therefore we do not expect the strength of scattering or the anisotropy of the IISM to change during a single orbit. However epoch-to-epoch variation must be expected if the line-of-sight velocity through the scattering region is $\gtrsim 10$ km s*-1*, since in the $\sim 100$ days between epochs it will traverse the refractive scale $s_{\rm r}$.

We compute the weighted mean and the rms of $\Delta\nu_{\rm d}$ at each epoch, but we find that the apparent value of $\Delta\nu_{\rm d}$ depends on the channel bandwidth $B_{\rm c}$. In a few epochs we have observations with $B_{\rm c}=0.78$ MHz and also with $B_{\rm c}=0.25$ MHz. In these cases $\Delta\nu_{\rm d}$ should be the same, but the estimator determined from the ACF saturates near $B_{\rm c}$ and typically is close to 0.3 MHz. Unfortunately many of our observations are made with $B_{\rm c}>\Delta\nu_{\rm d}$. In this case fitting the ACF does not provide a useful estimate, but we know that the intensity variance $V_{I}$ will be reduced by a factor of $\Delta\nu_{\rm d}/B_{\rm c}$. When $B_{\rm c}<\Delta\nu_{\rm d}$ we know that $V_{I}=M_{I}^{2}$ where $M_{I}$ is the mean intensity. So when $B_{\rm c}>\Delta\nu_{\rm d}$ we can use $\Delta\nu_{\rm d}=F_{\rm c}B_{\rm c}V_{I}/M_{I}^{2}$, where $F_{\rm c}$ is a calibration factor that is needed because the method depends on the actual shape of the ACF. We use the epochs for which we have observations with two different values of $B_{\rm c}$ to determine $F_{\rm c}=0.82$. This is shown in Figure 2, where the direct $\Delta\nu_{\rm d}$ measurements are in dark grey, the scaled estimates from $V_{I}$ are in light grey, and the resulting weighted mean $\Delta\nu_{\rm d}$ for each epoch is in black. This technique was also used briefly by Kerr et al. (2018), but they had insufficient data to determine the calibration factor $F_{\rm c}$. We show here for the first time that the method works well enough to detect small changes in bandwidth and generally agrees well with the direct measurement of $\Delta\nu_{\rm d}$ from the ACF.

3 The models

In this section we describe how we use the procedure of Cordes & Rickett (1998) and a model of the line-of-sight velocity to describe the time variation of $\tau_{\rm d}$ and $\Delta\nu_{\rm d}$. The equations describing the variation of $\tau_{\rm d}$ over a binary orbit are very non-linear in the physical parameters and it is not obvious a priori how many parameters can be estimated from one orbit. However Rickett et al. (2014) showed that the orbital variation of $\tau_{\rm d}$ when the scattering is anisotropic, can be described by 5 harmonic coefficients. In addition one can measure $\Delta\nu_{\rm d}$ which provides an independent degree of freedom. So at best one can model only six parameters which must include: $s_{\rm d}$, the axial ratio of anisotropy $A_{\rm r}$, the angle of anisotropy $\psi$, the distance to screen $s$, and two constant components of the velocity $V_{x}$, and $V_{y}$. If $A_{\rm r}$ is small, then $\omega$ and $i$ can also be measured at each epoch. Rickett et al. (2014) also showed that it was possible to use the shape of the ACF to independently measure the anisotropy, thus providing two more degrees of freedom. However our ACF estimates are not sufficiently accurate to take advantage of this option.

For a long-term analysis we can also use the variation of the Earth’s velocity over the year, and the known relativistic advance of the longitude of periastron (from precise pulsar timing; Bhat et al., 2008), to provide additional constraints. Of course one must re-estimate $s_{\rm d}$ from $\Delta\nu_{\rm d}$ at each epoch because the strength of scattering is likely to change on a spatial scale of AU; the estimated size of the scattering region.

The model of Cordes & Rickett (1998) is designed to account for variations in the strength of scattering, which will change the $\Delta\nu_{\rm d}$ and consequently $\tau_{\rm d}$. They define a “scintillation velocity" $V_{\rm ISS}$, which is the ratio of the spatial scale of the diffraction pattern at the observer $s_{\rm d}$, to the temporal scale $\tau_{\rm d}$. This $V_{\rm ISS}$ can be modelled if one knows the distribution and velocity of interstellar plasma along the line-of-sight. We considered two simple models for the plasma distribution:

Uniform medium: A continuous, uniform distribution of plasma with Kolmogorov turbulence along the line-of-sight.

Thin screen: A single compact “blob" of plasma with Kolmogorov turbulence, the “scattering screen", at some position $s$ between the pulsar at $s=0$ and the Earth at $s=1$.

We also assume that the density irregularities originate from turbulence and are described by a Kolmogorov spectrum. Cordes & Rickett (1998) derive the spatial scale $s_{\rm d}$ given the observed decorrelation bandwidth $\Delta\nu_{\rm d}$ for both of these models, assuming that the scattering is isotropic. In general the scintillation velocity is given by

[TABLE]

where $D$ is the distance to the pulsar in kiloparsecs, $f$ is the observing frequency in GHz, $\Delta\nu_{\rm d}$ is in MHz, and $\tau_{\rm d}$ is in seconds. The factor $A_{\rm ISS}$ depends on the assumed geometry of the scattering medium and on the exponent of the density spectrum. For a uniform medium, $A_{\rm ISS}=2.53\times 10^{4}$ km s*-1*, while for a thin screen, $A_{\rm ISS}=2.78\times 10^{4}\sqrt{2(1-s)/s}$ km s*-1* (Cordes & Rickett, 1998). We also include the extension to anisotropic scattering that was first presented by Coles et al. (2005) in both of our models.

$V_{\rm ISS}$ (Equation 3) is related to the effective transverse line-of-sight velocity $V_{\rm eff}(s)$ through the scattering medium at position $s$, which is a linear combination of the pulsar, Earth, and IISM velocities:

[TABLE]

where $V_{\rm p}$, $V_{\mu}$, $V_{\rm E}$, and $V_{\rm IISM}$ are the velocities from the pulsar’s orbit, the pulsar proper motion, the Earth, and the IISM respectively. Each is relative to the Solar system barycentre, and the IISM velocity can vary as a function of distance along the line-of-sight, $s$, for extended scattering. In Appendix A we describe the model for $V_{\rm eff}$, which includes the definitions for the pulsar orbital velocity, the orientation of the orbit in celestial coordinates, and the extension to anisotropic scattering developed by Coles et al. (2005).

In addition to these models that we use for the long-term scintillation, we also consider the short-term model of Ord et al. (2002b) for each of the 23 epochs of observations in Figure 3. This physical model necessarily assumes isotropic scattering, but we also fit the harmonic coefficient model of Rickett et al. (2014), which is useful for understanding sources of noise in our long-term analysis. Both of these models are summarised in Section 3.1 with results given in Section 5.1.

Each scintillation velocity model includes a single constant component, which (from Equation 4) is $V_{C}=(1-s)V_{\mu}-V_{\rm{IISM}}$ for a thin screen (for a uniform medium there is an appropriate integration over $s$ with $V_{\rm{IISM}}(s)$; Appendix A, Equation 9). The proper motion of PSR J1141$-$6545 is not currently known from pulsar timing, however we expect that the transverse velocity of the pulsar system ($V_{\mu}$) is larger than any IISM velocity ($V_{\rm{IISM}}$), since the latter may be of order $\sim 10$ km s*-1*. Therefore if we parameterise $V_{C}$ with only $V_{\mu}$ in our models (i.e. by setting $V_{\rm{IISM}}=0$ in Equation 4), the measurement will include a contamination from any non-zero $V_{\rm IISM}$. So when interpreting our measurement of $V_{\mu}$ as the proper motion, we would implicitly make the assumption that the velocity of the scattering screen is small compared with $V_{\mu}$, i.e. $V_{\rm{IISM}}<<(1-s)V_{\mu}$. To account for this, we include a conservative $10$ km s*-1* factor added in quadrature to the measurement uncertainty of each component of $V_{\mu}$.

The parameters in our long-term $V_{\rm{eff}}$ model that are unknown from pulsar timing of PSR J1141$-$6545 are $s$, $i$, $\Omega$, $v_{\mu,\parallel}$, $v_{\mu,\perp}$, $R$, and $\psi$. The subscripts $\parallel$ and $\perp$ denote coordinates parallel and perpendicular to the line of nodes of the pulsar’s orbit (Appendix A). We also require a scaling factor $\kappa$ to account for any systematic errors in the calculation of $s_{\rm d}$ from $\Delta\nu_{\rm d}$. The constraints from the known Earth’s velocity and $\dot{\omega}$ provide the additional degrees of freedom required to uniquely determine these parameters in a long-term analysis. For the uniform medium model we have one fewer parameter because we integrate over $s$ (Equation 9 in Appendix A). For the case of isotropic scattering, which we consider for both the uniform medium model and the thin screen model, we have $R=0$ (reducing the model to Equation 15 in Appendix A, scaled by the factor $\kappa$). However, we also found in Section 2.3 that it is necessary to re-estimate the spatial scale from $\Delta\nu_{\rm d}$ at each epoch because of AU-scale inhomogeneities in the IISM.

The calculation of the spatial scale from these $\Delta\nu_{\rm d}$ measurements is imperfect, and Rickett et al. (2014) found some disagreement between $\Delta\nu_{\rm d}$ and the spatial scale (by analysing the harmonic coefficients described in the following section). To account for this, and some other correlated time-variability in properties such as the anisotropy and IISM velocity, we chose to use a scaling factor at each epoch. To do this, we scaled the data for each epoch to the mean $V_{\rm ISS}$, then fitted for the mean scaling factor $\kappa$ between the model and scaled data. Ideally $\kappa$ would then represent the major systematic biases in the model, such as from an inaccurate assumed pulsar distance, while the scaling at each epoch accounts for time variations in $s_{\rm d}$ and perhaps some of the variation due to changing anisotropy and IISM velocity.

3.1 Short-term models of earlier work

For the first scintillation analysis of PSR J1141$-$6545, Ord et al. (2002b) modelled the dynamic spectrum for a single 10 hr observation. From this dynamic spectrum they derived the $V_{\rm ISS}$ of Cordes & Rickett (1998) (Equation 3), and modelled this using only the pulsar’s orbital velocity components ($v_{\mu,\parallel}$ and $v_{\mu,\perp}$; Appendix A, Equation 11), with a scaling factor $\kappa_{u}$. They assumed that the scattering was isotropic and the velocity model was then given by

[TABLE]

with $\kappa_{u}$, $i$, $\omega$, $v_{\rm C,\parallel}$, and $v_{\rm C,\perp}$ being the five fitted parameters. We use $v_{\rm C,\parallel}$ and $v_{\rm C,\perp}$ here instead of $v_{\mu,\parallel}$ and $v_{\mu,\perp}$ because they will include a contribution from the Earth’s velocity, and as a result will show annual variations with time. The parameter $\kappa_{u}$ can be used to absorb any errors in the calculation of $V_{\rm ISS}$, such as from an incorrect pulsar distance, which is what we intend for $\kappa$ of the long-term models (Equation 17). However, since the pulsar velocity in this model is not scaled by a screen distance $s$ because the Earth’s velocity is not included to constrain it, one could also interpret $\kappa_{u}$ in this case as a scale factor for $A_{\rm ISS}$. In this interpretation, $\kappa_{u}$ is also related to the screen distance $s$ and can account for different scattering geometries provided they are isotropic.

We repeat this analysis of Ord et al. (2002b) as described above for each of the 23 epochs of observations and we refer to this as the single-epoch “physical model". We present the results from this model in Section 5.1 and show that this approach can also be used to measure the advance of periastron $\dot{\omega}$.

Like PSR J1141$-$6545, which we analyse here, the double pulsar PSR J0737$-$3039A shows orbital modulation of the scintillation timescale, which was first measured and modelled by Ransom et al. (2004). The analysis was extended to include anisotropy in the IISM by Coles et al. (2005) and further extended to include the Earth’s velocity by Rickett et al. (2014), who also showed that the modulation of $1/\tau_{\rm d}^{2}$ with orbital phase can be modelled as the sum of five harmonics

[TABLE]

where the harmonic coefficients contain all of the information on the diffractive interstellar scintillation available in the data. The relationship between these coefficients and the physical parameters of the scattering and velocity are given in Equation 10 of Rickett et al. (2014) and are reproduced in Appendix A (Equation 18)

It is useful to note that each of the coefficients are inversely proportional to $s_{\rm d}^{2}$, and that $K_{S2}$ and $K_{C2}$ are constant with time for a constant anisotropy. Rickett et al. (2014) used these facts to normalise the coefficients by $K_{C2}$, which corrects each one for changes to the spatial scale $s_{\rm d}$ with time and with observing frequency. This allowed them to model observations at multiple observing frequencies simultaneously.

We have measured these normalised harmonic coefficients ($k_{0}=K_{0}/K_{C2}$, $k_{S}=K_{S}/K_{C2}$, $k_{C}=K_{C}/K_{C2}$, and $k_{S2}=K_{S2}/K_{C2}$) for each epoch of PSR J1141$-$6545 observations. The results are shown in Figure 4. It is worth noting that in the case of isotropic scattering (with $a=b=1$ and $c=0$), $k_{S}$ only changes with time because of $V_{\rm C,\parallel}$ (ignoring $\dot{\omega}$), $k_{C}$ changes with $V_{\rm C,\perp}$, and $k_{S2}=0$. While we do find that $k_{S2}$ is consistent with zero on average (Figure 4), there is some variation that is correlated with variations in the other normalised parameters. Since $k_{S2}$ is constant with orbital phase (for a constant anisotropy over the orbit), this indicates epoch-to-epoch variation to the anisotropy, which is an unmodelled source of noise in our data (as well as most scintillation observations). Because $k_{S2}=0$ on average, we can also see that our long-term model must be consistent with $c=0$, and thus should indicate that the scattering is isotropic, with $R=0$, or anisotropic with orientation near $\psi=0^{\circ}$ or $\psi=90^{\circ}$.

The harmonic coefficient model for each epoch was found with a linear least squares fit, while the physical model of Ord et al. (2002b) was fitted with nonlinear regression (using the lsqnonlin function in Matlab) because this was a convenient way to propagate uncertainties to the physical parameters. For the long-term models we use Bayesian methods for determining parameter uncertainties and performing model selection, as described in the following Section.

4 Bayesian inference and model selection

Bayesian methods have recently become popular for analysing astrophysical datasets because they can provide a robust means for parameter estimation, model selection, and visualisation of parameter correlations. An introduction to Bayesian inference, with examples taken from gravitational-wave astronomy, is given in Thrane & Talbot (2018). In brief, we aim to estimate the posterior probability distribution $P(\theta|d)$ for our set of model parameters $\theta$, given our data $d$, using Bayes theorem:

[TABLE]

where $P(d|\theta)$ is the likelihood function of the data given the set of parameters of a model, $P(\theta)$ is the prior probability distribution for the parameters, and $Z\equiv P(d)=\int P(d|\theta)P(\theta)d\theta$ is the fully-marginalised likelihood function or “evidence". For a single model, the evidence $Z$ is a normalisation factor that is often ignored, but it becomes meaningful when comparing evidences for multiple models. For each model of scattering geometry, we sample $P(\theta|d)$ and estimate $Z$ with the sampler of Veitch & Vecchio (2010), implemented in MATLAB by Pitkin & Romano (2013), with 1000 initial sample points and a multivariate Gaussian likelihood function.

We initially used a uniform prior for each parameter across a physically-motivated bound region: $\Omega\in[0,2\pi)$, $\cos{i}\in[-1,1]$, $v_{\mu,\parallel},v_{\mu,\perp}\in[-2000,2000]$ km s*-1*, $R\in[0,1)$, $\psi\in[0,\pi)$, $\kappa\in(0,3]$, and $s\in(0,1)$. However, for the anisotropic uniform medium model, we found that the resulting solution was non-physical (because the measured inclination angle $i=48.6\pm 1.9^{\circ}$ is ruled out by timing analysis), suggesting that the model is incorrect. To test this, we also consider each of these models with a Gaussian prior probability distribution on $i$ with mean and standard deviation $73\pm 3^{\circ}$111Although each model with a uniform $\cos{i}$ prior returned a unique solution with $i<90^{\circ}$, we also considered the opposite sense for $i$, with a Gaussian prior of $107\pm 3^{\circ}$ to quantify the strength of evidence for the measurement of the sense.. This prior corresponds to the 2$\sigma$ range derived using GR in Bhat et al. (2008). We choose this 2$\sigma$ range because the estimate in Bhat et al. (2008) did not account for correlations between $i$ and other parameters in the timing model, because it was not a true measurement of the Shapiro delay.

We also found that there are significant correlations between $i$, $v_{\mu,\perp}$, and $R$ (Figure 8, which has also been reported previously, e.g. Coles et al., 2005; Rickett et al., 2014). Using the Gaussian prior distribution on $i$ allows us to constrain the models to something physical, given that the anisotropy is completely unconstrained a priori. In addition $v_{\mu,\perp}$ and $R$ likely have some unmodelled time-dependence, which may contribute to a bias in these parameters and correlated parameters such as $i$. Using this Gaussian prior will therefore give a more reliable estimate of the anisotropy and pulsar velocity.

4.1 Bayes factors

For the purpose of model selection, we choose the null hypothesis to be isotropic scattering by a thin screen with a uniform prior for $\cos{i}$. This can be considered as the standard model for scintillation velocity analysis, but this is often out of necessity because of degrees-of-freedom constraints. We use the Bayes factor $K$ to quantify the evidence for (or against) our test models. This Bayes factor is the ratio of model evidences, which we choose to quote as the difference between log-evidences

[TABLE]

where $Z_{0}$ and $Z_{h}$ are the evidence for the null hypothesis and test model respectively.

There are numerous suggestions for how to interpret the strength of evidence for a model, given a Bayes factor (e.g. Jeffreys, 1998; Kass & Raftery, 1995). $|\log{K}|>8$ is commonly used as a significance threshold, while $|\log{K}|>20$ can be considered “strong" evidence for one model over the other, with the sign indicating which model is favoured (Kass & Raftery, 1995). However, for this work we find that the magnitude of Bayes factors are much greater than these contentious thresholds.

In Table LABEL:tab:compmodel we show the $\log{K}$ for each long-term model with respect to the isotropic thin screen. In the top half of this table we show the values obtained with a uniform prior for $\cos{i}$. This shows that the isotropic uniform medium is decisively ruled out, while the anisotropic thin screen and anisotropic uniform medium are strongly preferred by the data. Although there is strong evidence for the anisotropic thin screen over the anisotropic uniform medium, with $\log{K}=28.3$, we also reconsidered these models with a Gaussian prior on $i$ to help constrain these models to physical solutions. The Bayes factors for these models are in the lower half of Table LABEL:tab:compmodel and show that only an anisotropic thin screen is favoured over the isotropic screen.

Finally, to quantify the evidence for our measurement of the sense of the inclination angle, we also considered Gaussian priors of $107\pm 3^{\circ}$ with each model. Each model strongly favours $i<90^{\circ}$, for example comparing the two priors for the anisotropic thin screen we have $\log{K}=178.6$ in favour of $73\pm 3^{\circ}$.

The decisive evidence for anisotropic scattering reflects the fact that a long-term analysis such as this is highly sensitive to changes in $\tau_{\rm d}$ as a function of the position angle of $V_{\rm eff}$. This is particularly true because the annual variation is completely known except for a scale factor, and this is also the reason for such strong evidence for the thin screen model and inclination $<90^{\circ}$. This highlights the importance of annual variations for modelling scintillation velocities.

5 Results

Here we present and compare the various scintillation models for PSR J1141$-$6545. We first fitted each of the 23 epochs separately with a sum of five harmonics of orbital phase, which are described in Rickett et al. (2014) and in Section 3.1. We then fit a physical model for isotropic scattering to each epoch (Section 3.1), which provided measurements of $\omega$ and $V_{\rm C}$ over time to clearly show that the relativistic advance of periastron and the modulation from the Earth’s velocity can be recovered in the data. These models are presented below in Section 5.1.

We then held the periastron advance fixed at the precisely-measured value from pulsar timing, $\dot{\omega}=5.3096$ ${}^{\circ}{\rm yr}^{-1}$ (Bhat et al., 2008), and calculated the components of the Earth’s velocity ($V_{\rm E}$) transverse to the line-of-sight of PSR J1141$-$6545 for each observation. We combined the 23 epochs of data and fit several long-term scintillation models, using the additional degrees of freedom provided by $\dot{\omega}$ and $V_{\rm E}$ to help constrain additional parameters including the scattering anisotropy ($R$ and $\psi$) and the longitude of the ascending node $\Omega$. We fit these long-term velocity models for two distributions of Kolmogorov turbulent plasma along the line-of sight: A uniform distribution, and a thin screen. The uniform medium models are ruled out with high significance, as is purely isotropic scattering (see previous Section). We present our long-term anisotropic thin screen model in Section 5.2 for two prior probability distributions on the inclination angle.

Finally, we have included in Appendix B the links to access the data (raw observations and processed dynamic spectra) used for this work, as well as the MATLAB code used for the analysis for the purpose of reproducibility, including the scripts used to measure scintillation parameters from the dynamic spectra and fit the scintillation velocity as described in the previous section.

5.1 Individual epochs

For each epoch of observations shown in Figure 3, we fit the $V_{\rm ISS}$ of Equation 3 as a function of orbital phase $\phi$, with a sum of five harmonics (Section 3.1). The resulting normalised harmonic coefficients are given in Figure 4 as a function of time. In the case of PSR J0737$-$3039A, the inclination of the orbit is essentially edge-on, which means that $k_{C}$ and $k_{S2}$ are almost zero (Rickett et al., 2014). This is not true for PSR J1141$-$6545, but we do still find that the mean of $k_{S2}$ is close to zero. This suggests that the scattering may be nearly isotropic since $k_{S2}$ is independent of the changing pulsar and Earth velocities. There is however some time variability in $k_{S2}$, which is correlated with some other normalised harmonic coefficients. This may be due to time variability in the anisotropy, which one would expect even from random realisations of a truly isotropic medium. This time variability is a source of noise for the long-term scintillation models. This may be the primary reason for a large reduced chi-squared value ($\chi_{r}^{2}\sim 2.3$) for the long-term models (e.g. Figure 7).

Considering the evidence for isotropic scattering is important because the single-epoch physical model (Ord et al., 2002b) depends on it. With only five degrees of freedom available from a single epoch of observations, accurate measurement of the pulsar proper motion and inclination angle can only be made for isotropic scattering (Coles et al., 2005). We therefore fit this physical model to each of the epochs under the assumption of isotropic scattering to obtain a time series for each of the physical parameters described in Section 3.1.

There are four possible solutions for each epoch, arising from a known degeneracy between the proper motion and inclination angle (Lyne, 1984). This produces a more edge-on solution with higher proper motion and a more face-on solution with lower proper motion, and their corresponding pairs with an opposite “sense" of inclination about $i=90^{\circ}$. Ord et al. (2002b) considered only the $i<90^{\circ}$ solutions and determined the more edge-on solution to be favorable using the implied pulsar mass. For each of our fits we also take the $i<90^{\circ}$ solutions, and we show in Sections 4.1 and 5.2 that this is physical. The measurement of $\omega$ is unaffected by these degeneracies.

The time series of our fitted $\omega$ values is shown in Figure 5. The clear gradient is the relativistic advance of periastron, $\dot{\omega}=5.6\pm 0.3$ ${}^{\circ}/{\rm yr}$, and is close to the measurement of pulsar timing $\dot{\omega}=5.3096\pm 0.0004$ ${}^{\circ}/{\rm yr}$ (Bhat et al., 2008). We find that $\kappa_{u}$ and $i$ are constant with time with the exception of random variations similar to those seen in the harmonic coefficients (Figure 4). The locus of the two solutions in inclination angle from these degenerate solutions are $i=80.1\pm 1.1^{\circ}$ and $i=71\pm 3^{\circ}$, and the magnitudes of their corresponding velocities are $76\pm 3$ km s*-1* and $44\pm 2$ km s*-1* respectively. These inclination angles are both consistent with the previous scintillation measurement of $76\pm 2.5^{\circ}$, to within about 1$\sigma$ of both measurements.

The time-series of $V_{C}$ components for the more edge-on solution is shown in Figure 6 to demonstrate the clear annual variation, however we do not comment on which solution is physical because these parameters are highly correlated with $R$ and are therefore both biassed under the assumption of isotropic scattering. The reduced chi-squared value for these solutions (across all epochs, with total number of parameters $m=115$; five per epoch), is $\sim\chi^{2}_{r}=1.5$, suggesting that either the measurement errors for scintillation parameters are underestimated or there is excess noise in the data. The $\chi^{2}_{r}$ for individual epochs with observations spanning several days is generally higher than the few epochs with a single day, which is usually close to unity. It is therefore likely that the cause is small random changes in the scattering on a timescale of $\sim$days, such as variations to the anisotropy (which we would expect even from isotropic scattering).

5.2 Long-term model: Anisotropic thin screen

The annual variations of $V_{\rm ISS}$, which are clear in the single-epoch models shown in Figure 6, are due to the known $V_{\rm E}$ and can be included in our long-term model. This provides additional constraints that break the usual degeneracies between $i$ and $v_{\mu,\perp}$, and in the sense of $i$. It also provides the additional degrees of freedom required to measure the scattering anisotropy.

Any anisotropy cannot be accounted for by measuring a single epoch of $V_{\rm ISS}$ alone because there are too few degrees of freedom available for the additional two parameters. Instead, Coles et al. (2005) was able to use the correlated scintillations of both pulsars in the double pulsar system at their apparent closest approach (when the magnetosphere of pulsar B eclipses pulsar A) to produce a spatial correlation pattern that revealed the anisotropy. Later, Rickett et al. (2014) was able to use the annual variations in the harmonic coefficients for PSR J0737$-$3039A to measure the anisotropy and determine the sense of the inclination angle by comparing the two models with inclinations fixed by the $\sin{i}$ measurement from pulsar timing. They also showed that the anisotropy and IISM velocity (because the proper motion velocity was known from timing) could be measured from a model of the individual ACFs alone, but that the measurements differed from the average value obtained with long-term fitting, suggesting some weakness in the model and difficultly in confidently determining the anisotropy. The annual variation approach is equivalent to the long-term physical models we use here, but we do not have sufficient frequency resolution to use the ACFs of PSR J1141$-$6545 for an independent anisotropy estimation. Fortunately though, our data are sensitive to the changing spatial scale for different orientations of the rotating velocity vector, and we have been able to rule-out a purely isotropic model with high significance.

Our best-fit long-term model, which has an anisotropic thin screen and includes a Gaussian prior on inclination angle, is shown as $V_{\rm ISS}(\phi)$ in Figure 7. The model parameters are given in the second column of Table LABEL:tab:bestparams. The full posterior probability distributions for these parameters are shown in Figure 8. Each observation epoch has been independently scaled to the mean $V_{\rm ISS}$ in Figure 7 as described in Section 3, and the apparent white noise in the model is actually the out-of-phase variation due to the Earth’s velocity. For this model we find $i=68.6\pm 1.1^{\circ}$, which is more than $1\sigma$ inconsistent with the prior, suggesting a strong preference in the data for a lower inclination angle. The anisotropy for this model is $R=0.64\pm 0.03$, which corresponds to an axial ratio $A_{r}=2.14\pm 0.11$. This is comparable in magnitude to what has been observed for scattering towards other objects (e.g. Coles et al., 2005; Frail et al., 1994).

If the scattering were truly isotropic, we would expect that each observation samples a single realisation of this, and thus would randomly appear slightly anisotropic, with an rms of $R_{\rm rms}\sim 0.7(s_{\rm d}/s_{\rm r})^{-1/6}$ (Romani et al., 1986; Rickett et al., 2014). In our case $R_{\rm rms}\sim 0.18$, and from 23 epochs we would expect to observe $R_{\rm rms}/\sqrt{23}\sim 0.04$. The measurement is significantly larger than this, which suggests that there is significant anisotropy in the scattering. We find $\psi=86.8\pm 0.4^{\circ}$, which is a rather close alignment with the pulsar orbit. Given that the IISM is independent of the pulsar, $\psi$ must be uniformly distributed over 180∘. So the probability of such alignment is only 3.6%. However including the anisotropy greatly improves the fit and the accuracy of all the derived parameters, so we have to assign the alignment to chance.

Anisotropy of this magnitude, in the regime of strong scattering, is expected to produce a set of reversed sub-arclets distributed along a primary arc in the secondary spectra (e.g. Walker et al., 2004; Cordes et al., 2006; Brisken et al., 2010). However, without fine sampling in the dynamic spectra, this collection of arclets would together appear as broad unfocussed power and may not have a clear curvature. Indeed our analysis of the secondary spectra failed to show any clear scintillation arcs, which may be due to a combination of the sampling characteristics of the dynamic spectra and the broadening of the arcs through reverse sub-arclets. Future observations and studies of PSR J1141$-$6545 could be designed to maximise the detail in the ACFs and secondary spectra to allow for an independent measurements of the anisotropy and its time variability, which could then be used in these long-term $V_{\rm ISS}$ models to improve the $i$ and $v_{\mu,\perp}$ measurements.

The reduced chi-squared value for the best models is $\chi_{r}^{2}\sim 2.3$, and (as mentioned for the single-epoch models as well) this is likely high because of model errors, such as our assumption of a time-stationary anisotropy and IISM component in $v_{\mu}$. Indeed the normalised harmonic coefficients in Figure 4 show some time variability that would not be explained by the Earth’s motion or by changes to the strength of scattering that would be taken into account by scaling of each epoch. This random time-variability is likely to be our largest source of noise and the major reason for a large $\chi_{r}^{2}$, but there may also be some contribution from measurement errors on the scintillation parameters for example. Obtaining independent constraints on the anisotropy and IISM velocity from modelling high-resolution ACFs directly (e.g. Rickett et al., 2014) is the best chance for future work to improve the quality of the fit.

6 Discussion

Our modelling of the long-term scintillation of PSR J1141$-$6545 shows that the dominant scattering region is centred at a fractional distance of $s=0.276\pm 0.008$ from the pulsar, and that the scattering is slightly anisotropic with an axial ratio of $A_{r}=2.14\pm 0.11$. We were able to provide an independent estimate of the orbital inclination angle $i=61.3\pm 2^{\circ}$ from the anisotropic thin screen model using uniform priors on $\cos{i}$, however this was inconsistent with the value inferred from pulsar timing. We repeated the fit with a Gaussian prior on the inclination angle, to obtain improved measurements of the other correlated scintillation parameters. The measured parameters from both of these models, the derived anisotropy axial ratio $A_{r}$, the mass of pulsar $M_{\rm p}$, and the mass of the companion $M_{\rm c}$ are given in Table LABEL:tab:bestparams. We also derive a new estimate of the pulsar distance $D$ and proper motion in celestial coordinates for the first time, as described below in Section 6.1. We then discuss the implications of these derived measurements for future tests of GR in Section 6.2.

6.1 Pulsar distance and proper motion estimates

The distance to PSR J1141$-$6545 is currently poorly constrained, with the best estimate of $D=3\pm 2$ kpc from Verbiest et al. (2012) originating from a luminosity bias correction to an earlier lower-limit of $3.7\pm 1.7$ kpc that was derived from neutral hydrogen absorption spectrum (Ord et al., 2002a). There is also a lower estimate of 1.7 kpc from the most recent Galactic electron-density model (Yao et al., 2017). This poor measurement precision is problematic for our data because our derivation of $V_{\rm ISS}$ from the scintillation parameters is proportional to $\sqrt{D}$. We accordingly fitted for a scaling factor $\kappa$ (see Section 3), which would absorb any errors in the calculation of $V_{\rm ISS}$, for example from an error in the assumed pulsar distance of $D=3$ kpc. However, this scaling factor would also include systematic errors in the measurements of $\Delta\nu_{\rm d}$ and $\tau_{\rm d}$, as well as any error in the numerical relationship between $\Delta\nu_{\rm d}$, the strength of scattering, and $s_{\rm d}$ derived by Cordes & Rickett (1998). This numerical relationship is summarised by the $A_{\rm ISS}$ coefficient, $A_{\rm ISS}=2.78\times 10^{4}\sqrt{2(1-s)/s}$ km s*-1* for the thin screen model.

We have also chosen to scale each epoch of observations to the mean $V_{\rm ISS}$, which will primarily take into account any epoch-to-epoch errors in the relationship between the strength of scattering and $\Delta\nu_{\rm d}$. In Section 2.3 we showed that it was necessary to re-estimate $\Delta\nu_{\rm d}$ at each epoch because of inhomogeneities in the IISM on an AU scale. However we also know that the derivation of $s_{\rm d}$ from $\Delta\nu_{\rm d}$ is imperfect and can depend on the shape of the ACF, among other factors. This was also noted by Rickett et al. (2014), who found discrepancies between their $\Delta\nu_{\rm d}$ measurements and the spatial scale inferred from measurements of the harmonic coefficient $K_{C2}$ with time. For our analysis, we intend for our scaling of each epoch to take into account such discrepancies with time, which may also come partially from time-evolution in other parameters, such as the anisotropy. However, we assume that the scaling factor for the mean $V_{\rm ISS}$, $\kappa$, represents the largest systematic error in the calculation of $V_{\rm ISS}$, and we assume that this is from an error in the pulsar distance $D$.

For the two best models, we show this scaling factor $\kappa$ and its error in Table LABEL:tab:bestparams. For our preferred model with a Gaussian inclination angle prior, $\kappa=0.496\pm 0.014$, suggesting a significant systematic error in the $V_{\rm ISS}$. If this is indeed a measurement of the error in pulsar distance, we can use this to provide a new estimate. To do so, we first assume a 20% uncertainty on $A_{\rm ISS}$, originating from a discrepancy of approximately this magnitude between the spatial scale $s_{\rm d}$ derived from the measured $\Delta\nu_{\rm d}$, and that derived from $K_{C2}$ in this work and that of Rickett et al. (2014). We then calculate the distribution for pulsar distance given that $D\propto 1/\kappa^{2}$ and assume a Gaussian distribution for $\kappa$ (which is reasonable, from Figure 8) with mean and standard deviation given by the measurement in Table LABEL:tab:bestparams. We do this for both models (for comparison), and find that the distance is estimated to be larger than the assumed $3$ kpc, but that it is poorly constrained by this method, $D=10^{+4}_{-3}$ kpc.

If the pulsar is this distant, it may be surprising that we observe significant anisotropy and localised scattering. But while most pulsars that have been studied with scintillation are relatively nearby, there are examples of distant thin screens observed out to of order few kpc (e.g. Putney & Stinebring, 2006; Popov et al., 2016). Indeed, in this work we have precisely measured the fractional distance to a scattering screen, which would be placed at $\sim 2$ kpc even for the lower suggested pulsar distance (Verbiest et al., 2012; Yao et al., 2017). We find no clearly associated HII regions for this line-of-sight in a search of the WISE catalogue (Anderson et al., 2014) (the nearest likely region has a radius 8.6 arcmin centred $\sim$30 arcmin from the line-of-sight). Our estimate shows that in principle scintillation can be used to estimate pulsar distances, but in practice it is complicated by the $D\propto 1/\kappa^{2}$ relationship and by the large uncertainty for $A_{\rm ISS}$. With improved scintillation modelling and understanding of the IISM along the line-of-sight, we expect that the precision on distance estimates with this method can be improved in the future.

Using our own distance estimate (for self-consistency) with our measurements of the pulsar’s transverse velocity $v_{\mu,\parallel}$ and $v_{\mu,\perp}$, and the longitude of the ascending node $\Omega$, we have been able to derive the pulsar proper motion in celestial coordinates (right ascension $\alpha$ and declination $\delta$). We used the same distance distribution derived from $\kappa$ above, and calculated the distribution of proper motions assuming Gaussian distributions for $v_{\mu,\parallel}$, $v_{\mu,\perp}$, and $\Omega$ given by the measurements and their uncertainties listed in Table LABEL:tab:bestparams. The derived proper motions in $\alpha$ and $\delta$ are shown in the table for both models, and are consistent because of the large uncertainty. For the Gaussian $i$ prior model, we have $\mu_{\alpha}\cos{\delta}=2.9\pm 1.0$ mas yr*-1* and $\mu_{\delta}=1.8\pm 0.6$ mas yr*-1*. We believe that this is a fairly conservative estimate (because of the large uncertainty on $D$), but may be under-estimated if other existing distance estimates are correct. For example, taking the current lowest distance estimate of $D=1.7$ kpc from Yao et al. (2017) and assuming a 20% uncertainty, we would have $\mu_{\alpha}\cos{\delta}=15\pm 4$ mas yr*-1* and $\mu_{\delta}=9\pm 3$ mas yr*-1*. However a proper motion of this magnitude may soon be ruled out with improved pulsar timing sensitivity.

The estimated proper motion is highly uncertain, primarily because of uncertainty in the pulsar distance. For the following section we will assume the values derived from our new distance estimate. However, if in the future the pulsar distance is constrained with higher confidence, the proper motion should be re-derived from our measurements of $v_{\mu,\parallel}$, $v_{\mu,\perp}$, and $\Omega$ in Table LABEL:tab:bestparams. Alternatively, given a measurement of the pulsar proper motion in the near future from improved timing precision, our independent measurements of the velocity can be used to derived the distance.

6.2 Implications for timing and tests of gravitational theories

PSR J1141$-$6545 is a highly relativistic pulsar in an eccentric, asymmetrical mass system, which makes it an ideal laboratory for testing GR. Bhat et al. (2008) analysed the gravitational radiation losses from this system through pulsar timing and noted that with increasing precision of the orbital period-derivative, contamination from kinematic effects (e.g. the Shklovskii effect) and Galactic acceleration would start to dominate the uncertainty in the near future. The transverse velocity of the pulsar system is accompanied by a radial acceleration, which produces a time-dependent Doppler-shift to the pulsar spin frequency and orbital period. This is the Shklovskii effect (Shklovskii, 1970), and it results in an apparent orbital period-derivative $\dot{P}_{b}^{\rm kin}$ that is considered a contamination to the orbital period-derivative measurement from gravitational radiation losses, $\dot{P}_{b}^{\rm GR}$. However, the exact level of this effect was unknown because the proper motion was not measured via pulsar timing. We are now able to determine the contribution from the Shklovskii effect for the first time because our scintillation work is sensitive to the transverse motion of the pulsar instead of the radial motion probed by pulsar timing.

This Shklovskii effect can be calculated from $\dot{P}_{b}^{\rm kin}=DP_{b}\mu^{2}/c=P_{b}V^{2}_{\mu}/cD$ (Bell & Bailes, 1996). With the distance and transverse velocity provided by our model, we calculate $P_{b}^{\rm kin}=3.5\times 10^{-15}$, which is $\sim$1% of the $\dot{P}_{b}$ measurement of Bhat et al. (2008) and well below the 6% measurement precision. However, if we again assume the distance from Yao et al. (2017), then $\dot{P}_{b}^{\rm kin}=1.9\times 10^{-14}$. This is at the level of the current expected timing precision of $\sim$2% for $\dot{P}_{b}$. It is therefore important to further constrain the system transverse velocity and/or distance from improved scintillation modelling and/or pulsar timing. In addition to the Shklovskii effect on $\dot{P}_{b}$, the proper motion changes the projected geometry of the binary orbit, resulting in apparent $\dot{x}$ and $\dot{\omega}$. However, these kinematic contaminations are far below the measurement precision for this pulsar and are typically only observed in precisely-timed millisecond pulsars (Arzoumanian et al., 1996; Kopeikin, 1996).

The ongoing long-term timing campaign on this system will also benefit from the robust measurements of the sense of the inclination angle and $\Omega$ obtained from this work. These measurements will help understand relative contributions to the secular variations of the orbit from phenomena such as relativistic aberration, spin-orbit coupling, and proper motion. Measurements of $\Omega$ will also be helpful in understanding the pulsar emission beam geometry with the precessional evolution of its polarisation (Kramer & Wex, 2009). Measurements of $\Omega$ combined with $\dot{\omega}$ from timing, can be used to test preferred-frame effects predicted by a variety of alternative theories of gravity (Wex & Kramer, 2007; Shao & Wex, 2012).

7 Conclusion

We have presented new scintillation models for PSR J1141$-$6545 using six years of data from the Parkes 64 m radio telescope. We found that like many pulsars, the scattering shows some anisotropy, and is dominated by a single scattering region centred at $s=0.276\pm 0.008$. By accounting for anisotropy in the scattering, we measured the system inclination angle of $i=61.3\pm 2^{\circ}$, which is significantly lower than the constraint of $73\pm 3^{\circ}$ from pulsar timing (inferred from the pulsar and companion masses derived with GR), suggesting that a weakness in the model, such as time-stationary IISM velocity and anisotropy may cause some systematic errors. However, by using this pulsar timing constraint as a prior probability distribution, we have been able to measure several astrometric parameters for the first time.

Using the significant annual and relativistic variations observed in the scintillation velocity to constrain long-term models, we have been able to resolve the “sense" of the inclination angle, and we find that $i<90^{\circ}$. This in turn resolved the ambiguity in the direction of the proper motion velocity in pulsar coordinates. With our new measurement of the orientation of the orbit in celestial coordinates $\Omega=24.8\pm 1.8^{\circ}$, and estimate of the pulsar distance $D=10^{+4}_{-3}$ kpc, we have been able to estimate the proper motion for the first time. We determine the proper motion in right ascension $\mu_{\alpha}\cos{\delta}=2.9\pm 1.0$ mas yr*-1* and in declination $\mu_{\delta}=1.8\pm 0.6$ mas yr*-1*, and we use these numbers to calculate the contribution of the Shklovskii effect to $\dot{P}_{b}$. This effect is the most significant source of contamination for tests of GR, but our low proper motion suggests that it exists only at the $\sim$1% level. Our improved accuracy and precision for the pulsar’s transverse velocity is also important for understanding the formation of this system.

Our distance measurement is model-dependent, with the accuracy determined partially by the relationship between the measured scintillation bandwidth and the spatial scale, for which we have used the method of Cordes & Rickett (1998). In the near future, wide-bandwidth observing systems, such as that of the MeerKAT radio telescope or the ultra-wideband low-frequency (UWL) receiver for the Parkes 64 m radio telescope, could provide an experimentally-derived relationship through analysis of the frequency-dependence on the scintillation bandwidth. In this way scintillation studies may be used to give distance measurements to pulsars with predictable modulation of the scintillation timescale (e.g. relativistic binaries).

We may also soon see an independent proper motion measurement from pulsar timing with improved observation span and techniques, which would provide another method for determining the pulsar distance in combination with our velocity measurements. Future high-quality observations, with the sensitivity and resolution to estimate the scattering anisotropy directly from the autocovariance functions of individual dynamic spectra (as in Rickett et al., 2014) will be valuable for improving the accuracy and precision of future inclination angle measurements of scintillating binary pulsars.

Acknowledgements

We thank Barney Rickett for useful discussion and for allowing us access to code and data used for the earlier work in Rickett et al. (2014). 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 CSIRO. Work at NRL is supported by NASA. MB is supported by the ARC Centre of Excellence OzGrav under grant CE170100004.

Appendix A: Scintillation velocity models

From the dynamic spectrum of intensity scintillations as a function of frequency and time, one can derive a scintillation velocity $V_{\rm ISS}$ (Equation 3) that relates to the structure and velocity of the diffraction pattern drifting across the observer. A line-of-sight velocity model for $V_{\rm ISS}$ is found by integrating $V_{\rm eff}$ (Equation 4) along the line-of-sight from $x=0$ at the pulsar and $x=1$ at the observer, with a weight that corresponds to the geometry of the scattering medium (Cordes & Rickett, 1998)

[TABLE]

where $\alpha=5/3$ for a Kolmogorov medium, and $\eta(x)$ is the mean-square scattering angle per unit distance and functions as the weight for the integral to describe different geometries. For uniform scattering along the line-of-sight $\eta(x)=1$, while for a thin screen $\eta(x)$ is a delta function at $x=s$ and the expectation of $V_{\rm ISS}$ reduces to the line-of-sight velocity with respect to the diffraction pattern $V_{\rm los}$, defined at the location of the observer,

[TABLE]

where $E(V_{\rm ISS})$ denotes the expectation of the observationally-derived $V_{\rm ISS}$, and $V_{\rm{eff}}$ is given in Equation 4. The transverse orbital velocity of the pulsar $V_{\rm p}$ has components $v_{{\rm p},\parallel}$ along the line of nodes, and $v_{{\rm p},\perp}$ perpendicular to this in the plane of the sky. These velocities are defined as a function of orbital phase from the line of nodes $\phi=\theta+\omega$, where $\theta$ is the true anomaly and $\omega$ is the longitude of periastron,

[TABLE]

where $V_{0}=2\pi xc/(\sin{i}P_{b}(1-e^{2})^{1/2})$ is the mean orbital velocity, $x$ is the projected semi-major axis in seconds, $P_{b}$ is the binary orbital period, $e$ is the eccentricity, $i$ is the inclination angle. The true anomaly is first calculated by numerically computing the eccentric anomaly, $E$ from Kepler’s equation $E-e\sin{E}=M$, with mean anomaly $M=(2\pi/P_{b})(t-T_{0})$, where $T_{0}$ is the epoch of periastron. The true anomaly $\theta$ is then given by

[TABLE]

The known contribution of the Earth’s velocity is then rotated from celestial coordinates to the pulsar-frame coordinates according to the longitude of the ascending node $\Omega$ for the pulsar’s orbit

[TABLE]

where $v_{\rm E,\parallel}$ and $v_{\rm E,\perp}$ are the components of the Earth’s velocity aligned with $v_{{\rm p},\parallel}$ and $v_{{\rm p},\perp}$ respectively, and $\Omega$ is defined as a rotation East of North. This was also used by Rickett et al. (2014) (but with the direction of the perpendicular axis reversed in their definition) to include the Earth’s velocity and proper motion of PSR J0737$-$3039A into their scintillation model. We then combine the pulsar and Earth velocities, include the pulsar system transverse velocity ($v_{{\mu},\parallel}$ and $v_{{\mu},\perp}$), and scale them appropriately by the distance to the scattering region $s$

[TABLE]

In the case of isotropic scattering where the angular size of the pulsar orbit is compact enough to remain in the scattering disk, the spatial scale $s_{\rm d}$ is constant with orbital phase (but it may change on longer time scales if the IISM is inhomogeneous; see Section 2.3). Consequently, the decorrelation bandwidth $\Delta\nu_{\rm d}$ is also constant with orbital phase and the effective velocity is simply given by

[TABLE]

However, for anisotropic scattering, $s_{\rm d}$ depends on the direction of $V_{\rm ISS}$. To account for such scattering in PSR J0737$-$3039A, Coles et al. (2005) considered the spatial diffraction pattern as an ellipse. The pattern is then described by a quadratic form $Q(\mathbf{s_{\rm d}})=as_{\rm d,\parallel}^{2}+bs_{\rm d,\perp}^{2}+cs_{\rm d,\parallel}s_{\rm d,\perp}$, where the coefficients $a$, $b$, and $c$ are parametrised by the axial ratio $A_{r}$ of the ellipse and its orientation $\psi$ with respect to the coordinates of the pulsar orbit as defined above. Rickett et al. (2014) used this anisotropy model for PSR J0737$-$3039A but parametrised the quadratic coefficients in terms of $R=(A_{r}^{2}-1)/(A_{r}^{2}+1)$, which is bound between 0 and 1. If the orientation angle $\psi$ is defined clockwise from the line of nodes, then from Rickett et al. (2014) the coefficients are

[TABLE]

Finally, we introduce a scaling factor $\kappa$ to the model, which will account for any errors (for example an error in the pulsar distance $D$) in the calculation of $V_{\rm ISS}$ from the dynamic spectrum (Equation 3). Our final model for the effective velocity is then

[TABLE]

which we use in Equation 9 with $\eta(x)=1$ for the uniform medium model and Equation 10 for the thin screen model.

This anisotropic scattering model can alternatively be represented as a sum of five harmonics, which we used to model each of the individual epochs. This model is described in Section 3.1, but the full dependence on physical parameters for each of the harmonic coefficients in Equation 6 is given below (Rickett et al., 2014)

[TABLE]

where $v_{\rm C,\parallel}$ and $v_{\rm C,\perp}$ are the constant (for a given epoch) components of the line-of-sight velocity, which will be dominated by the pulsar’s proper motion, but also includes the Earth’s velocity and any IISM velocity.

Appendix B: Reproducing our results

The raw data from the Parkes radio telescope that were used for this work are available from the CSIRO data access portal (DAP; https://data.csiro.au), and were processed with a pipeline developed for the second PPTA data release. This pipeline is used to produce TOAs and dynamic spectra, and is briefly summarised in Section 2. The dynamic spectra files that were produced from this pipeline, and the collection of MATLAB codes used to analyse these dynamic spectra for the results presented in this paper, are available from the DAP at https://doi.org/10.4225/08/5ae7b7b1b65c8.

The code is presented for reproducibility and in general is not intended to be used for other applications, however a brief description of each script and function is included in README.txt, and some may be useful for other scintillation studies. For example, we have included a function (getDynspecParams.m) that may be used to measure scintillation parameters for any dynamic spectrum given as a two-dimensional matrix of intensity versus observing time and frequency, using the methods described in Section 2.2. We have also included the scintillation velocity model (vissmodel.m), which is described in Section 3, and a script that numerically calculates the true anomaly from MJD and then uses nonlinear regression to fit the model to $V_{\rm ISS}$ derived from scintillation parameters (modelDynspecParams.m).

The MATLAB code was executed in version 2017b, and requires jd2date.m from the Astromatlab library (https://webhome.weizmann.ac.il/home/eofek/matlab/; Ofek, 2014), MATLAB Multinest (Pitkin & Romano, 2013), and the “Image Processing", “Statistics and Machine Learning", “Signal Processing", “Optimization", “Curve Fitting", “Symbolic Math", and (optionally) “Parallel Computing" MATLAB toolboxes.