https://orcid.org/0000-0003-4453-3776
ABSTRACT
Using data from the Large European Array for Pulsars, and the Effelsberg telescope, we study the scintillation parameters of the millisecond pulsar PSR J0613−0200 over a 7 yr timespan. The ‘secondary spectrum’ – the 2D power spectrum of scintillation – presents the scattered power as a function of time delay, and contains the relative velocities of the pulsar, observer, and scattering material. We detect a persistent parabolic scintillation arc, suggesting scattering is dominated by a thin, anisotropic region. The scattering is poorly described by a simple exponential tail, with excess power at high delays; we measure significant, detectable scattered power at times out to |${\sim}5 \, \mu {\rm s}$|, and measure the bulk scattering delay to be between 50 to 200 ns with particularly strong scattering throughout 2013. These delays are too small to detect a change of the pulse profile shape, yet they would change the times of arrival as measured through pulsar timing. The arc curvature varies annually, and is well fitted by a one-dimensional scattering screen |${\sim}40{{\ \rm per\ cent}}$| of the way towards the pulsar, with a changing orientation during the increased scattering in 2013. Effects of uncorrected scattering will introduce time delays correlated over time in individual pulsars, and may need to be considered in gravitational wave analyses. Pulsar timing programmes would benefit from simultaneously recording in a way that scintillation can be resolved, in order to monitor the variable time delays caused by multipath propagation.
1 INTRODUCTION
Radio emission from pulsars experiences several propagation effects from the ionized interstellar medium (ISM), as the index of refraction varies with electron density and frequency. The signal acquires a group delay Δt, known as dispersion, scaling as |$\Delta t \propto \mathrm{DM}\, \nu ^{-2}$|, where DM is the integrated column density of free electrons, and ν is the observing frequency. Spatial variations in the electron density result in multipath propagation, with deflected paths acquiring a geometric time delay from the path-length difference compared to the direct line of sight. When these delays are large (compared to the pulse duration), it is observed as scattering, a one-sided broadening of pulses often resembling an exponential tail. When these delays are small, we observe it as scintillation, the constructive and destructive interference of different deflected images at the observer, resulting in a time and frequency dependence of the observed flux. These delays are steeper in frequency than dispersion, and are expected to scale roughly as τ ∼ ν−4.
One of the central goals of pulsar timing is to directly detect gravitational waves, in a so-called pulsar timing array (PTA, Hobbs 2013; Desvignes et al. 2016; Verbiest et al. 2016; Arzoumanian et al. 2018). The most stable pulsars are observed on weekly to monthly cadence over many years, and ∼nHz gravitational waves could be observed in timing residuals correlated in time and position on the sky (Hellings & Downs 1983). This effect is expected to be tiny, with a fractional change of the arrival time compared to the gravitational wavelength of order 10−15, so it requires careful understanding of all other effects that would change the arrival times of pulses. While PTA pulsars are selected for their stability, they all experience variable dispersion and scattering to some degree due to the relative motion of the pulsar and observer with respect to the ISM. Variable dispersion measures have been measured and corrected using multifrequency data (e.g. Keith et al. 2013), while changes in scattering time are often estimated using the statistical relation between the scintillation bandwidth (the frequency width of scintillation) and scattering time (e.g. Levin et al. 2016; Shapiro-Albert et al. 2020, see Verbiest & Shaifullah 2018 for a review of how these effects limit precision pulsar timing). Dispersion and scattering both scale strongly with frequency, and are often covariant. One can look for variable delays following a ν−2 and ν−4 scaling (Lam et al. 2019); the technique of wide-band template matching has recently been developed as a way to jointly fit for these effects (Liu et al. 2014; Pennucci, Demorest & Ransom 2014; Pennucci 2019; Alam et al. 2020b).
In this paper, we begin to apply the methods of Hemberger & Stinebring (2008) to PTA pulsars, in which scintillation arcs are used to estimate time delays from multipath propagation. We analyse PSR J0613−0200 over 7 yr, in roughly monthly cadence, using data from the Large European Array for Pulsars (LEAP) (Stappers & Kramer 2011; Bassa et al. 2016), and a 3-month bi-weekly observing campaign using the 100-m Effelsberg radio telescope. This pulsar is of particular interest; it shows the strongest evidence of a 15 nHz strain, but since the signal appears most strongly in this pulsar, it is believed to arise from an unmodelled non-GW signal (Aggarwal et al. 2019). In Section 2, we give an overview of some necessary background of scintillation, and summarize the methods of Hemberger & Stinebring (2008). In Section 3, we describe our observations with the LEAP telescope, our short-term observing campaign with the Effelsberg telescope. In Section 4, we outline our methods, in Section 5 we present our results, and we discuss the ramifications and future prospects in Section 6.
2 BACKGROUND ON THEORY OF SCINTILLATION
2.1 Thin screen theory and stationary phase approximation
The theory of scattering in thin screens is outlined in detail in Walker et al. (2004) and Cordes et al. (2006), and we summarize some of the pertinent relations here.
2.2 The interstellar response
In this section, we summarize and expand upon the method of Hemberger & Stinebring (2008), to use the secondary spectrum to estimate the total time delays from multipath propagation.
2.3 Estimating time delays from the secondary spectrum
Now we address how to estimate the time delays of the intensity response function in practice. As we are only concerned with measuring time delays, in this section we drop the time dependence for simplicity, focusing on the imprint of the impulse response function on the spectrum.
We see that in the limit of a strong central image, we can recover the total time delays from the secondary spectrum. More generically, how well the time delays can be computed from the secondary spectrum is dependent on the unknown distribution of images [or the functional form of gI(τ)]. In the case of strong scattering, there is no reason to expect a single undeflected line-of-sight image, but rather there may be many bright, scattered images at small angular separation. In this case, the time delay computed from the above formula will be overestimated, due to the cross-terms of bright central images interfering. This will bias the result high by a factor of ∼2m/(m + 1), where m is the number of bright images, leading to a difference as large as a factor of 2. In the case of a discretized secondary spectrum, this will only begin to matter if the image separations are larger than one pixel in τ, otherwise it will approximate the case of one bright central image. Additionally, we are still limited by the fact that the secondary spectrum measures differences in time delays, rather than absolute time delays; if there is a time-shift applied to all images, it would not be captured by our estimate.
With the above caveats mentioned, we use equations (20) and (23) as our basis to measure time delays throughout the paper. These include the assumption of a strong central image, which we believe gives a reasonable estimate for our purposes. We describe how to compute time delays in practice from our data in Section 4.3, after detailing our data reduction and secondary spectra creation.
3 OBSERVATIONS
3.1 LEAP
The LEAP is a phased array of five large radio telescopes in Europe; the Effelsberg telescope, the Lovell telescope at Jodrell Bank Observatory, the Westerbork Synthesis Radio Telescope, the Nançay Radio Telescope, and the Sardinia Radio Telescope (Stappers & Kramer 2011). The coherent addition of radio signals from all these telescopes results in an effective 195 m diameter dish. The overview of LEAP is given in Bassa et al. (2016). Observations have been made monthly since 2012, with whichever subset of these telescopes was available.
The voltage data from each site are shipped or transferred to Jodrell Bank Observatory to be correlated and coherently added on a designated CPU cluster, using a specifically designed software correlator (details in Smits et al. 2017). Correlation involves a polarization calibration based on an observation of PSR J1022+1001 or PSR B1933+16 from the same epoch, correlation on a calibrator to find an initial phasing solution, then self-calibration on the pulsar to determine the time delays and fringe drift rates for each telescope throughout the observation, using Effelsberg as a time and position reference. The coherently added voltages are stored on tape, allowing us to re-reduce the data with arbitrary time or frequency resolution. The high sensitivity, and the flexibility offered by storing the baseband data has enabled LEAP to do single pulse studies of MSPs (Liu et al. 2016; McKee et al. 2019); for these same reasons, it is an ideal telescope for the scintillation work presented in this paper. Typical observing lengths are 30−60 minutes, with bandwidths of 80−128 MHz (comprised of 16 MHz subbands), depending on the subset of telescopes used for a given observation. As we will show in Section 5.2, the angular extent of the scattering screen is unresolved by LEAP, so we can safely treat it as a single-dish instrument for our purposes.
3.2 Effelsberg 100-m Telescope
From March to June 2020, we had a roughly bi-weekly monitoring campaign using the Effelsberg telescope. Baseband data were recorded as 8-bit ‘dada’2 files using the PSRIX backend (described in Lazarus et al. 2016), using the central feed of the 7-beam receiver (‘P217mm’). The data were recorded in 25 MHz subbands, with a usable bandwidth of 1250–1450 MHz, and typical observation lengths of |$90\,$| min. While Effelsberg alone is less sensitive than LEAP, this is compensated through the larger exposure times and bandwidth.
4 METHODS
4.1 Creating dynamic and secondary spectra
We created folded archives from the baseband data using dspsr (van Straten & Bailes 2011), coherently de-dispersing and folding with 10 s bins, 128 phase bins, and sufficient channels to fully resolve scintillation −62.5 and 50.0 kHz channels for LEAP and Effelsberg, respectively. The subbands were combined in frequency using the psrchive tool psradd (Hotan, van Straten & Manchester 2004) to form one combined archive per observation. The following processing steps for data from either telescope are identical unless expressly stated otherwise.
After summing polarizations, each folded archive contains a data cube I(t, ν, phase). We use a fixed off-pulse region relative to the pulse, a contiguous |$50{{\ \rm per\ cent}}$| section with no apparent pulsed emission (in Fig. 2, phase 0.5–1.0) We divide by the time average of the off-pulse region across the full observation to approximately remove the bandpass, and in each time and frequency bin, we subtract the mean of the off-pulse region to remove variable background flux. Subintegrations with an off-pulse standard deviation >5× the mean rms value were masked, as were any time bins or frequency channels with |${\gt}30{{\ \rm per\ cent}}$| of flagged subintegrations.
The LEAP dada files are saved separately in each subband, in individual 10 s files; a small number of these files were missing, and were filled with zeros, and included in our mask. To reduce artefacts caused by Fourier transforming over a window function, masked pixels were iteratively in-painted using the mean of the nearest pixels. While more sophisticated methods of inpainting exist, this is sufficient for our analysis, as typically no more than |$5{{\ \rm per\ cent}}$| of data are flagged.
A time and frequency averaged profile was created, and zeroed everywhere the S/N was below 5|$\, \sigma$|. This profile was used to weight each phase bin, before summing over pulse phase to create the dynamic spectrum I(t, ν). Over a narrow band, it is sufficient to simply use a 2D FFT, which we used for this analysis (over a wider band, the ν−2 scaling of η causes arcs to smear in the secondary spectrum, summarized in Gwinn & Sosenko 2019). Before taking a FFT, we padded the edges by a factor of two with the mean value of the dynamic spectra, to mitigate artefacts caused by edge effects. A few representative LEAP dynamic and secondary spectra, at the same time of several years, are shown in Fig. 1.
Top: Dynamic spectra of five observations around the same time of year, to have comparable contributions from the Earth’s velocity. The colourbar extends from 2σ below the mean to 5σ above. Bottom: Corresponding secondary spectra, with a logarithmic colourbar extending three orders of magnitude. Clear arcs with noticeable localized clumps of power are seen, these correspond to prominent diagonal features in the above dynamic spectra. The observation from 2013 is anomalous, showing extremely fine stripes in the dynamic spectrum, corresponding to power at large time delays.
4.2 Measuring arc curvatures
The main power in the secondary spectrum of PSR J0613−0200 follows a parabolic arc, suggesting scattering dominated by a highly anisotropic, thin screen. As described in Section 5.2, the arc curvature encodes veff and deff; we wish to measure the arc curvature for each observation, to probe the changing velocities of the system, and to localize power for measuring time delays. One method often used to measure parabolic curvatures is the Hough transform, finding the peak of the power summed over different possible parabolic curvatures (Bhat et al. 2016). While this technique works very well for thin parabolic arcs, it leads to a broad maximum (and correspondingly large uncertainties) for broad parabolae, as seen in our data. We determine the curvature in a different way; in τ steps of |$0.125\, \mu$|s in the secondary spectrum (where |$\tau \gt 0.5\, \mu$| s, to avoid confusion in the bright centre), we find the peak value of I(fD) for both positive and negative fD. We keep only points where the peak is >4× the rms of the background, estimated from the region of |$I(|f_{D}| \gt 10\, \mu \text{s})$|. This set of points in fD and τ is fitted with a parabola, using an orthogonal minimization routine, to find the best-fitting curvature and error.
4.3 Integrating the secondary spectrum
Artefacts in the dynamic spectrum, such as radio-frequency interference (RFI), phasing imperfections, and the window function lead to correlated features in the secondary spectrum. As such, the noise properties are not always well behaved, and direct error propagation underestimates the error on 〈τ〉. We estimate our errors directly from the cumulative function in equation (24); at high enough T the integral plateaus, with residual variations caused by the effect of integrating noise in the secondary spectrum. We take the mean and standard deviation of equation (24) between |$T=4\!-\!8\, \mu$|s as our measurement and error of 〈τ〉, respectively.
4.4 ‘Timing’ a convolved template
To illustrate the effects of scattering on a profile, we can directly convolve our measure of the amplitude of gI(τ) into a template profile and measure the time offset using the standard Fourier template-matching algorithm outlined in the appendix of Taylor (1992). We create an analytic template using the standard psrchive tool paas (Hotan et al. 2004), fitting the profile with a series of von Mises functions, and interpolate the solution to have the equivalent 31.25 ns bins of our measured |I(τ)|2. We convolve the two, and measure the relative time delay between the convolved template against the original one. Fig. 2 shows this convolution applied to one of our observations. The measured time delay in this way agrees perfectly with the method in the previous section; timing recovers the shift correctly, even when the effects are not visibly noticeable.
The effects of scattering on a pulse profile. Top: Analytic profile, before and after convolving with the scattering tail measured from scintillation. The inset panel shows |$\tilde{I}(\tau)$| measured on 2017 February 17 in the way described in Section 4.3. The resulting convolved profile looks identical by eye to the original, as the time delays are largely subbin. Middle: The residual of the original and convolved profile, bottom: the residual between the template and shifted profile, after shifting them into closest alignment. The middle, uncorrected residuals are only as large as |$1{{\ \rm per\ cent}}$|, while after aligning the residuals are below the |$0.1{{\ \rm per\ cent}}$| level. This reinforces the fact that the main contribution of scattering is a bulk time shift, rather than a noticeable change in profile shape.
We note again that this is not precisely the intensity impulse response, but rather its autocorrelation, but it is close enough in amplitude to demonstrate that the convolved template is visually identical (with residuals at the |$0.1{{\ \rm per\ cent}}$| level after aligning the template), yet is measurably delayed. In addition, |I(τ)|2 is noticeably clumpy and poorly described by an exponential, even after being effectively smoothed by the autocorrelation.
4.5 Inferred time delay from the frequency ACF
A standard way to infer the time delays from scattering is to construct the autocorrelation functions (ACF) R(Δν) = (I*I)(Δν). Fitting the width (specifically, the HWHM) of the ACF in frequency gives the scintillation bandwidth νscint, which is inversely proportional to the bulk scattering delay as 〈τ〉 = C/2πνscint (C is commonly assumed to be 1, and depends on the assumptions of the scattering distribution). This method is often used when arcs cannot be resolved nicely, as Δν can typically be measured simply and robustly. However, if gE(t) is not smooth, then the ACF will be poorly described as a single Gaussian. Two such examples of an ACF, one well described, and one poorly described by a 1D Gaussian fit are shown in Fig. 3, along with their derived νscint and 〈τ〉. In the second case, a single Gaussian would preferentially fit the broad component, and result in an inferred time delay which is low, while power at large time delays results in the narrow peak smaller than 1 MHz.
5 RESULTS
5.1 Evolution of scattering time from LEAP: month to year time-scales
We present our measurements of the bulk time delay in the top panel of Fig. 4. We find significant persistent scattering at the |${\sim}80\,$| ns level, and a few cases of strong scattering variability on several month to year time-scales. The time-scales are set by the time it takes power to move through the secondary spectrum, as described in Section 2. The most striking feature is the strong excess scattering in 2013, where the bulk scattering is variable and extends above |$200\,$| ns. This is not captured very well by the ACF method, as the Gaussian fit latches on to the broad-scale scintillation rather than the narrow peak caused by the large time delays, as described in Section 4.5. As hinted at in Fig. 3, this could potentially be remedied by using a multicomponent model to the ACF, as in principle the ACF contains the equivalent information as the secondary spectrum, only differing by a Fourier transform.
Top: Variations in the bulk scattering time from every LEAP observation. Black points are estimated through integrating |$\tilde{I}(\tau)$|, as described in Section 4.3. The red points are delays inferred from the measurements of νscint from the ACF, described in Section 4.5. Middle: The measured arc curvatures for each observation, with a best-fitting model of a 1D scattering screen from 2014 onwards overplotted. Bottom: Measured DM values from NANOGrav 12.5 yr data release (Alam et al. 2020a).
In the bottom panel of Fig. 4, we plot the DM values of PSR J0613−0200 from NANOGrav’s 12.5 yr release (Alam et al. 2020a). The DM is steeply decreasing prior to 2013, and there is clear annual variation; this was studied in detail in Jones et al. (2017), with data spanning from 2006 until near the end of 2013. The authors fit the time variations of DM with a 1-yr period sinusoid and a linear trend, to capture the contribution of the pulsar’s observed trajectory through the ISM from the pulsar’s and Earth’s velocity, respectively. The residuals of the DM show a borderline significant dip of 2−3 × 10−3 pc cm−3 at MJDs 56300–56400, the beginning of 2013.
Although the ecliptic latitude of the pulsar is quite large (−25.4○), part of the annual variation that appears in the DM time series might be explained with the contribution of the Solar wind. The time of closest approach for PSR J0613−0200 is in mid-June each year, and by modelling the distribution of electrons in the Solar wind as spherical (Edwards, Hobbs & Manchester 2006) with an electron density of 7.9 cm−3 at the Earth orbit (Madison et al. 2019), we expect a DM displacement of the order of 2 × 10−4 pc cm−3 at the closest approach. By allowing the amplitude of the Solar wind approximation to vary year by year, we find a model that accounts for both the ISM and the Solar wind is preferred in 7 yr across the data set, while an ISM-only model favoured elsewhere. The years in which the complete model is preferred show a general compatibility with an amplitude of 7.9 cm−3, and after the subtraction of the time-dependent approximation of the Solar wind the most significant remaining feature is the steep gradient of the DM leading into 2013.
Similar events of increased scattering have been seen in PSR J1017−7156 and PSR J1603−7202 by the Parkes Pulsar Timing Array (PPTA), in which the scintillation bandwidth and time-scale decrease suddenly associated with a jump in DM of several 10−3 pc cm−3 , interpreted as an extreme scattering event (Coles et al. 2015). While we do not see an increase in DM of this order, the increased scattering we observe in PSR J0613−0200 may be of similar origin.
5.2 Location and nature of the scattering screen
As mentioned in Section 4.3, we fit the parabolic curvature of each observation to determine the masks for estimating the time delays. The arc curvatures contain the effective velocity, and vary throughout the year from the Earth’s motion. The existence of parabolic arcs suggests highly anisotropic scattering; for a one-dimensional screen, the arc curvature then depends only on the effective velocity parallel to the screen. By measuring the change in arc curvature over the year, one can measure the distance and orientation of the scattering screen.
We perform only a simple analysis here, currently ignoring the contribution from the pulsar’s orbital motion. We fit the observed curvature values beyond 2013 with a one-dimensional screen, using measured values of the pulsar’s distance of 780 ± 80 pc and proper motion of μra = 1.822 ± 0.008 mas yr−1, μdec = −10.355 ± 0.017 mas yr−1 from Desvignes et al. (2016). The three free parameters are the fractional screen distance s, the orientation of the screen ψ, and vism,||, the velocity of the scattering screen parallel to its axis of anisotropy. A 1D screen fits the data well, shown in the middle panel of Fig. 4, while the best-fitting values are in Table 1. Using the screen distance, and the largest detectable time delays of |$\tau \approx 5\, \mu$|s, the largest angular extent of the screen is |$\theta \approx \sqrt{2c \tau / d_{\rm eff}} \approx 3$| mas, smaller than the resolution of the longest baselines of LEAP.
Fit parameters of 1D screen to the arc curvatures, as defined in Section 5.2.
| s | ψscr (deg) | vism,|| (km s−1) | |
|---|---|---|---|
| 2014 onwards | 0.62 ± 0.06 | 16± 2 | −1.2 ± 2.5 |
| 2013 event | 0.58 ± 0.10 | -36 ± 9 | 12.8 ± 2.8 |
| s | ψscr (deg) | vism,|| (km s−1) | |
|---|---|---|---|
| 2014 onwards | 0.62 ± 0.06 | 16± 2 | −1.2 ± 2.5 |
| 2013 event | 0.58 ± 0.10 | -36 ± 9 | 12.8 ± 2.8 |
Fit parameters of 1D screen to the arc curvatures, as defined in Section 5.2.
| s | ψscr (deg) | vism,|| (km s−1) | |
|---|---|---|---|
| 2014 onwards | 0.62 ± 0.06 | 16± 2 | −1.2 ± 2.5 |
| 2013 event | 0.58 ± 0.10 | -36 ± 9 | 12.8 ± 2.8 |
| s | ψscr (deg) | vism,|| (km s−1) | |
|---|---|---|---|
| 2014 onwards | 0.62 ± 0.06 | 16± 2 | −1.2 ± 2.5 |
| 2013 event | 0.58 ± 0.10 | -36 ± 9 | 12.8 ± 2.8 |
During the increased scattering of 2013, the best-fitting model poorly matches the data. Here, we investigate this year separately. The secondary spectra spanning 2013 are shown sequentially in Fig. 5. Persistent clumps of power can be tracked throughout the year, and features can be seen to cross the line of fD = 0, indicating a changing sign of |$\boldsymbol {v}_{\rm eff}$|. From this, we can fit directly the signed value of |$1/\sqrt{\eta } \propto v_{\rm eff}$| (with the same free parameters s, ψ, and vism,||), where the locations of velocity zero-crossings are quite constraining. The measures of |$\boldsymbol {v}_{\rm eff}$|, and best-fitting model are shown in Fig. 6.
Secondary spectra of all observations covering the anomalous scattering in 2013, with a logarithmic colourbar extending four orders of magnitude. The power can be seen to cross the fD = 0 line, indicating the sign of |$\boldsymbol {v}_{\rm eff}$| is changing throughout the year, and power can be seen to travel from left to right along the parabola.
Effective velocity of 2013 scattering event, and best-fitting 1D model. The sign of veff was assigned to switch when the dominant power in the secondary spectra was clearly seen to cross the fD = 0 axis.
The best-fitting screen parameters for 2013, and for all data beyond 2013 are tabulated in Table 1. The distance of the screen is consistent between both fits, with the orientation and parallel screen velocity differing between the two. This implies that the strong scattering plausibly arises from the same physical region. In addition, the absolute velocity of the screen need not be changing, as the orientation differs and we are only sensitive to the component of the screen velocity parallel to ψ; the results of both fits are consistent with a screen velocity of |vscr| = 15 ± 2 km s−1 at an angle of ϕvel = -75 ○ ± 10○(East of North).
We note however that the models above are incomplete, as a proper treatment needs to include the binary motion of the pulsar, which we have neglected. The orbit is 1.2 d, and |$v_{\mathrm{ orb}} \sin (i) = 19.9\,$| km s−1. Each observation is much smaller in duration than the orbit, but is at an effectively random orbital phase. This will add scatter in the velocities, and thus the curvatures. To properly account for the orbital velocity, one would need to jointly fit for i and Ω - such a fit has been performed successfully on 16 years of arc curvature values of PSR J0437-4715 by Reardon et al. 2020, measuring the inclination with 0.3○ precision and the longitude of ascending node with 0.4○ precision. Regardless, a 1D screen is a good fit to the data beyond 2013 (where each year the curvature peaks around November and is minimal around May), where the annual variation is the strongest effect.
5.3 Arclet evolution with Effelsberg: scattering time on week to month time-scales
We investigate the variability of scattering on week to month time-scales with the Effelsberg observing campaign described in Section 3.2. The secondary spectra of all of our observations are shown in Fig. 7. A clear parabolic arc, with a hint of inverted arclets is seen, and can be seen to clearly move through the secondary spectrum from left to right. We show three examples of larger, zoomed in secondary spectra in Fig. 8 to emphasize these features. We estimate the total time delay from the secondary spectrum using the methods of Section 4.3, shown in the bottom panel of Fig. 9. The total time delays are consistent with what was found with LEAP’s monitoring, showing steady scattering at around 60–100 ns, decreasing slightly over 2 months.
Secondary spectra of our roughly bi-weekly monitoring campaign with Effelsberg, with a logarithmic colourbar extending three orders of magnitude. Power can be seen to travel from left to right along the parabola, most evident by following the power in the last six panels. The arc curvature does not abruptly change between observations, despite being at random orbital phases, suggesting that the scattering screen resulting in these arcs is not very sensitive to the pulsar’s orbit.
Zoomed in Effelsberg secondary spectra from Fig. 7 for three different epochs, panels 1, 2, and 4. The red dotted line shows our best-fitting parabolic curvature, while blue dotted inverted parabolae are plotted with the same curvature, with their apex on the main parabola – these are to guide the eye towards faint structures in the secondary spectra which may be inverted arclets, and to help visualize the general trend of power moving from left to right along the parabola. More sensitive, or longer observations will be required to reveal the possible structure of inverted arclets more clearly.
Arclet evolution, and time delays from bi-weekly observations with Effelsberg. Top: Magnification (fraction of total intensity) and time delay of one arclet seen moving through the secondary spectrum. Red crosses are expected positions of the undetected arclet. Middle: Contribution of the single arclet to the total time delay to the pulsar’s signal. Bottom: Estimate of the bulk scattering time inferred from the secondary spectra.
We also attempt to measure the effect of a single arclet, which is akin to contribution of a pulsar passing a single, compact point of scattering in the ISM, not unlike an echo. We track the position, and total fractional flux of the arclet seen travelling to the upper-right in the final seven panels of Fig. 7. This same feature first appears to be moving towards the origin in panels of MJD 58951−58958, although it is less prominent. We fit a flux centroid in an ellipse around the arc, to track its motion in fD and τ. We measure the fraction of the flux in the arc, compared to that integrated over the full parabola, which are plotted in the top panel of Fig. 9. The motion of the arclet unsurprisingly traces out a parabola over time, as seen in PSR B0834+06 (Hill et al. 2005), and is similar to echoes seen in the Crab pulsar (e.g. Backer, Wong & Valanju 2000; Lyne, Pritchard & Graham-Smith 2001). The strength of the arclet is quite asymmetric about the origin, and at its peak contains |${\approx}4{{\ \rm per\ cent}}$| of the total pulsar flux. The total time shift arising from this arclet can be estimated as 〈τ〉arclet ≈ τarclet(Iarclet/I), and is shown in the middle panel of Fig. 9. The contribution from a single arclet as it passes in front of the pulsar contributes a variable scattering of |${\sim}20\,$|ns over a 2-month period.
5.4 Comparison to earlier results
Previous analysis using the ACF in Levin et al. (2016) measures the time delay from scintillation to be 〈τ〉 = 11.7 ± 4.9 ns, monitoring this pulsar up until 2013 October. Additionally, Shapiro-Albert et al. (2020) estimate a time delay of 〈τ〉 = 43.6 ± 2.3 ns in a similar manner. The frequency channels used in these papers were 1.5625 MHz wide, a common standard in timing archives, and would have averaged over the fine scintillation structures due to power at high delays. Keith et al. (2013) estimate a scintillation bandwidth of 1.64 MHz, for which one would infer a time delay of 97 ns. This measurement is from before we have data, so we cannot directly compare, but this value is much closer to the order of the time delays we measure.
5.5 Effects of uncorrected scattering
In this section, we aim to make from our measurements a simple estimate of the single-pulsar gravitational wave signal arising from unaccounted scattering. We stress that this is a conservative upper limit, as we do not know the extent to which variable scattering is absorbed in red-noise modelling, or in measurements of DM. The GW strain h at a given periodicity P is related to the amplitude TOA variations δt as roughly h ∼ 2πδt/P. Since Aggarwal et al. (2019) find excess signal at |$15\,$| nHz, which is ∼2.1 yr, we estimate this using the long-term variable time delays from LEAP shown in Fig. 4. We perform a Lomb–Scargle periodogram on the measured values of 〈τ〉, and convert to a measure of h while taking into account the proper normalizations, shown in Fig. 10. The measured value at |$15\,$| nHz is ∼10−15, still an order of magnitude lower than the single pulsar limit of h = 9.7 × 10−15 from the EPTA (95 per cent upper limit, from table 1 in Lentati et al. 2015). As PTA upper limits are improved, scattering variations, if unaccounted for, may begin to limit the timing precision.
Estimated strain signal which would arise from uncorrected scattering, obtained from a Lomb–Scargle periodogram of the bulk time delays 〈τ〉 from Fig. 4.
6 CONCLUSIONS
We have measured variable time delays on a PTA millisecond pulsar, using similar methods to those laid out in Hemberger & Stinebring (2008). The next logical step will be to perform timing with scattering time-scales subtracted from TOAs, to see if this improves the timing residuals. One can apply this approach to study the variable scattering in many PTA pulsars. With LEAP we can re-reduce the data to whatever time and frequency resolution we like, but regular timing observations would benefit with a second reduction with fine frequency channels at the expense of phase bins. Going further, methods to obtain the interstellar response directly may be important, including holography (Walker et al. 2008), cyclic spectroscopy (Demorest 2011; Walker, Demorest & van Straten 2013; Palliyaguru et al. 2015; Dolch et al. 2020), or directly by using bright giant pulses in special cases (Main et al. 2017). Only these methods, in which the interstellar delays are measured directly (as opposed to measuring delay differences, as we do in the secondary spectrum) can retrieve overall delays that are not related to the characteristic time-scale of the scattering tail. New analysis techniques such as the θ–θ diagram (Sprenger et al. 2020), which expresses the secondary spectrum in terms of the angular coordinates on the scattering screen, may be useful as well. This technique can be used to precisely measure arc curvatures, and could be used to efficiently perform holography of 1D screens (Baker et al., in preparation).
Scattering is statistically expected to follow ∼λ4, yet scattering arising from discrete structures (observed as arclets, or localized clumps of power in secondary spectra) will be localized at a fixed τ as a function of wavelength. Observations over a wide frequency range will help to inform the amplitude scaling of arclets, and thus the contribution of discrete arclets to the total scattering time at different frequencies. In addition, scattering occurs from density gradients in the ISM, so the link between variable DM and scattering should be explored in more detail. In the case where the DM and scattering variations occur in the same scattering screen, they could potentially be inferred from the other quantity; a predictive model of scattering from DM (or vice versa) would be a great step towards removing these effects from timing observations.
Knowing the screen distance and orientation, the location of power in the secondary spectrum is predetermined, but the amplitudes are dependent on the physics of scattering and lens models. In recent years, some predictive models of scintillation properties have been developed (e.g. Simard & Pen 2018; Gwinn & Sosenko 2019), which can be tested using measurements tracking arclets in time and frequency.
We used annual variations in the arc curvature to determine properties of the scattering screen, while orbital variations were ignored. Orbital variations can give an additional orbital constraint, such as the inclination (including the ‘sense’; Rickett et al. 2014; Reardon et al. 2019; Reardon et al. 2020), and could possibly lead to precise pulsar distances (Boyle & Pen 2012). Such analysis would be greatly improved with precise, quantitative measurements of the arc curvature. Additionally, one could instantaneously measure the scattering screen’s properties using the multiple telescopes of LEAP, either using the visibilities (Brisken et al. 2010) or more simply using the inter-station time delays of the dynamic spectra (Simard et al. 2019a, while the method of combining of visibilities and intensities is outlined in Simard et al. 2019b). This is being investigated, and will be the subject of future work.
As scattering screens are likely much smaller than the angular separation between pulsars, scattering variations are very unlikely to directly correlate between pulsars in a way that mimics a Hellings & Downs curve (Hellings & Downs 1983). But while no direct correlation is expected between pulsars, it is possible that scattering is variable on similar time-scales if pulsar proper motions and distances are comparable, and if the screen distance is not at an extreme (i.e. not too close to the pulsar or to the Earth). Several of the EPTA pulsars show variable scintillation arcs, similar to those shown in this paper, and will be subject of future work. As PTAs become more sensitive, any PTA result relying on a small number of pulsars may need to consider the effects of variable scattering when interpreting the significance of a potential gravitational wave signal.
ACKNOWLEDGEMENTS
RAM thanks Daniel Reardon for useful advice on modelling variable arc curvatures, and Tim Sprenger for useful comments on the manuscript.
This work is supported by the ERC Advanced Grant ‘LEAP’, Grant Agreement Number 227947 (PI M. Kramer). The European Pulsar Timing Array (EPTA) is a collaboration between European Institutes, namely ASTRON (NL), INAF/Osservatorio di Cagliari (IT), the Max-Planck-Institut für Radioastronomie (GER), Nançay/Paris Observatory (FRA), The University of Manchester (UK), The University of Birmingham (UK), The University of Cambridge (UK), and The University of Bielefeld (GER), with an aim to provide high-precision pulsar timing to work towards the direct detection of low-frequency gravitational waves. The Effelsberg 100-m telescope is operated by the Max-Planck-Institut für Radioastronomie. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope is supported by a consolidated grant from the STFC in the UK. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Foundation for Radio Astronomy, ASTRON, with support from NWO. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique.
The Sardinia Radio Telescope (SRT) is funded by the Department of Universities and Research (MIUR), the Italian Space Agency (ASI), and the Autonomous Region of Sardinia (RAS), and is operated as a National Facility by the National Institute for Astrophysics (INAF). From Mar 2014 - Jan 2016, the SRT data were acquired as part of the Astronomical Validation of the SRT. We thus thank the SRT Astronomical Validation Team, and in particular: M. Burgay, S. Casu, R. Concu, A. Corongiu, E. Egron, N. Iacolina, A. Melis, A. Pellizzoni, and A. Trois. We also thank M. Pilia and A. Possenti for helping with observations.
SC acknowledges the support by the ANR Programme d’Investissement d’Avenir (PIA) under the FIRST-TF network (ANR-10-LABX-48-01) project. SC and IC acknowledge financial support from Programme National de Cosmologie and Galaxies (PNCG) and Programme National Hautes Energies (PNHE) funded by CNRS, CEA and CNES, France.
KL is supported by the European Research Council for the ERC Synergy Grant BlackHoleCam under contract no. 610058. KJL gratefully acknowledge support from National Basic Research Program of China, 973 Program, 2015CB857101, and NSFC 11373011. WWZ was supported by the National Natural Science Foundation of China Grant No. 11873067, the CAS-MPG LEGACY project and the Strategic Priority Research Program of the Chinese Academy of Sciences Grant No. XDB23000000.
DATA AVAILABILITY
The timing data used this article shall be shared on reasonable request to the corresponding author.
Footnotes
Equivalently, one must Fourier transform over a long enough timespan of E(t) fully encompassing g(t) (by a factor of 2, due to the Nyquist Theorem) – the longest time-scales correspond to the finest frequencies.
