On Frequency-Dependent Dispersion Measures and Extreme Scattering Events

TL;DR

This paper investigates the causes of frequency-dependent dispersion measures in pulsar signals, demonstrating that turbulence in the interstellar medium explains observed variations without requiring extreme scattering events, using LOFAR data.

Contribution

The study introduces methods to differentiate between turbulence-induced and lensing-induced frequency-dependent DMs, applying them to LOFAR observations of PSR J2219+4754.

Findings

Frequency-dependent DM variations are consistent with ISM turbulence.

No evidence found for extreme scattering events or small-scale lensing.

Possible presence of an enhanced scattering screen near one end of the line of sight.

Abstract

Radio emission propagating over an Earth-pulsar line of sight provides a unique probe of the intervening ionized interstellar medium (ISM). Variations in the integrated electron column density along this line of sight, or dispersion measure (DM), have been observed since shortly after the discovery of pulsars. As early as 2006, frequency-dependent dispersion measures have been observed and attributed to several possible causes. Ray-path averaging over different effective light-cone volumes through the turbulent ISM contributes to this effect as will DM misestimation due to radio propagation across compact lensing structures such as those caused by "extreme scattering events". We present methods to assess the variations in frequency-dependent dispersion measures due to the turbulent ISM versus these compact lensing structures along the line of sight. We analyze recent Low-Frequency Array (LOFAR) observations of PSR J2219+4754 to test the underlying physical mechanism of the observed frequency-dependent DM. Previous analyses have indicated the presence of strong lensing due to compact overdensities halfway between the Earth and pulsar. Instead we find the frequency dependence of the DM timeseries for PSR J2219+4754 is consistent with being due solely to ISM turbulence and there is no evidence for any extreme scattering event or small-scale lensing structure. The data show possible deviations from a uniform turbulent medium, suggesting that there may be an enhanced scattering screen near one of the two ends of the line of sight. We present this analysis as an example of the power of low-frequency observations to distinguish the underlying mechanisms in frequency-dependent propagation effects.

Equations

$$P_{\delta n_{\mathrm{e}}}=C_n^2{q}^{-\beta },q_1\le q\le q_2$$

$${\sigma }_{\mathrm{DM}(\nu )}({\nu }_1,{\nu }_2)$$

$${\sigma }_{\mathrm{DM}}(\tau )$$

$${\sigma }_{d\mathrm{DM}/dt}\approx \frac{{\sigma }_{\mathrm{DM}}(\tau )}{\tau },$$

$$R_{d\mathrm{DM}/dt}\equiv \frac{\left|d\mathrm{DM}/dt\right|}{{\sigma }_{d\mathrm{DM}/dt}}.$$

$$C_n^2(r)=1.8\times 10^{10}{\left(\frac{r}{10R_{\circ }}\right)}^{-3.66}{\mathrm{m}}^{-20/3},$$

$${\sigma }_{\mathrm{DM}(\nu ),\mathrm{solar}}=\sqrt{\frac{1}{⟨1/{\sigma }_{\mathrm{DM}(\nu )}^2(s)⟩}},$$

$$t_{\mathrm{DM}}\sim \frac{{\lambda }^2{r}_e{\mathrm{DM}}_l}{2\pi c}\sim \frac{{\lambda }^2{r}_e{n}_{e}L}{2\pi c},$$

$$t_{\mathrm{geo}}\sim \frac{D_l(1-D_l/D_p){\lambda }^4{r}_e^2({\mathrm{DM}}_l^{\prime }{)}^2}{8{\pi }^{2}c}\sim \frac{D_l(1-D_l/D_p){\lambda }^4{r}_e^2{n}_e^2}{8{\pi }^{2}c{\zeta }^2},$$

$$t_{\mathrm{bary}}\sim \frac{(1-D_l/D_p){\lambda }^2{r}_{\oplus }{r}_e{\mathrm{DM}}_l^{\prime }}{2\pi c}\sim \frac{(1-D_l/D_p){\lambda }^2{r}_{\oplus }{r}_e{n}_e}{2\pi c\zeta },$$

$$\delta \mathrm{DM}=-\frac{{ϵ}_{{\nu }_1}-{ϵ}_{{\nu }_2}}{K({\nu }_1^{-2}-{\nu }_2^{-2})},$$

$$D_{\mathrm{DM}}(\Delta t_{\mathrm{d}})=1.47\times 10^{-15}(\mathrm{pc}\mathrm{c}{\mathrm{m}}^{-3}{)}^2{\left(\frac{\nu }{\mathrm{GHz}}\right)}^2.$$

$$n_{\mathrm{ISS}}\approx \left(1+{\eta }_t\frac{T}{\Delta t_{\mathrm{d}}}\right)\left(1+{\eta }_{\nu }\frac{B}{\Delta {\nu }_{\mathrm{d}}}\right).$$

$${\sigma }_{\mathrm{pol}}\sim 1\mu \mathrm{s}\left(\frac{\epsilon }{0.1}\right)\left(\frac{{\pi }_{\mathrm{V}}}{0.1}\right)\left(\frac{W}{100\mu \mathrm{s}}\right),$$

$$D_{\mathrm{DM}}^{(1)}(\tau )$$

$$D_{\mathrm{DM}}^{(2)}(\tau )$$

$${\sigma }_{\delta t_{\infty }}={\left(\frac{{\sigma }_{{ϵ}_{{\nu }_1}}^2+r^4{\sigma }_{{ϵ}_{{\nu }_2}}^2}{(r^2-1{)}^2}\right)}^{1/2},$$

$$t_{\nu }=t_{\infty }+\frac{K\mathrm{DM}}{{\nu }^2}+t_{C,\nu }+{ϵ}_{\nu }.$$

$$\widehat{\mathrm{DM}}=\frac{t_{{\nu }_1}-t_{{\nu }_2}}{K({\nu }_1^{-2}-{\nu }_2^{-2})}.$$

$${\hat{t}}_{\infty }$$

$$=t_{{\nu }_2}-\frac{K\widehat{\mathrm{DM}}}{{\nu }_2^2}$$

$$\delta \mathrm{DM}\equiv \mathrm{DM}-\widehat{\mathrm{DM}}$$

$$\delta t_{\infty }\equiv t_{\infty }-{\hat{t}}_{\infty }$$

$$\delta t_{\infty }=\frac{{ϵ}_{{\nu }_1}-r^2{ϵ}_{{\nu }_2}}{r^2-1},$$

$${\sigma }_{\delta t_{\infty }}^2=⟨\delta t_{\infty }^2⟩=⟨{\left(\frac{{ϵ}_{{\nu }_1}-r^2{ϵ}_{{\nu }_2}}{r^2-1}\right)}^2⟩.$$

$${\sigma }_{\delta t_{\infty }}^2=\frac{{\sigma }_{{ϵ}_{{\nu }_1}}^2+r^4{\sigma }_{{ϵ}_{{\nu }_2}}^2}{(r^2-1{)}^2}.$$

$${\sigma }_{\delta t_{\infty }}=\sqrt{⟨\delta t_{\infty }^2⟩}={\sigma }_{{ϵ}_{\nu }}{\left(\frac{r^4+1}{r^4-2r^2+1}\right)}^{1/2}.$$

Full text

ISM: structure – pulsars: individual: PSR J2219$+$4754

1 Introduction

The propagation of radio pulsar emission through the interstellar medium (ISM) provides a unique probe of ionized material along the path of propagation (Rickett, 1990). Estimates of the dispersion measure (DM), or integrated electron column density, are obtained from precisely measuring the arrival times of pulse emission as a function of radio frequency (Stairs, 2002). Measurements of how these dispersive delays, along with other frequency-dependent delays, change over time allow us to build models of structures in the ISM on a wide range of spatial scales (Armstrong et al., 1995).

While temporal DM variations have been observed for decades (Rankin & Roberts, 1971; Phillips & Wolszczan, 1991, 1992; Backer et al., 1993; Hobbs et al., 2004), frequency-dependent DM due to differing volumes of the ISM probed has only been observed more recently (see e.g., for PSR B1937+21 in Ramachandran et al. 2006, Demorest 2007; see also in Pennucci 2015). A thorough treatment on the theory of frequency-dependent DM was given by Cordes et al. (2016). Ray-path averaging of the radio emission through different volumes, and therefore different electron content, of the turbulent ISM will result in a smoothing of the DM timeseries by a kernel that broadens rapidly at lower frequencies.

Besides dispersive delays, there are a number of other physical mechanisms that cause the propagating radio emission to be delayed. These delays vary as a function of frequency (Foster & Cordes, 1990; Clegg et al., 1998; Shannon & Cordes, 2017). One example is from lensing of the emission around compact over- or under-densities in the ISM (Clegg et al., 1998; Cordes et al., 2017). However, since typical pulsar timing models only account for dispersive delays, there will be errors associated with the estimated DM due to these additional delays (Cordes & Shannon, 2010; Lam et al., 2016b). Therefore, DM estimates taken at two separate frequencies may show differences solely due to biases from these other propagation effects and can be mis-attributed to a frequency dependence in the true dispersive delay.

Donner et al. (2019) have most recently reported three-and-a-half years of timing measurements of PSR J2219$+$4754 (B2217$+$47) with three stations of the International LOFAR (LOw-Frequency ARray) Telescope (ILT). From those data, they have determined frequency-dependent variations in the $\mathrm{DM}$. The temporal variations noted in the DM timeseries for PSR J2219$+$4754 for both the lower- and higher-frequency data qualitatively appear very similar to those predicted by Cordes et al. (2016) for a turbulent medium, especially in that the higher radio-frequency DM data show higher fluctuation/Fourier-frequency structure than the lower radio-frequency timeseries. Donner et al. (2019) concluded, however, that while the data are consistent with arising from the turbulent medium, the DM variations result from small-scale ($\sim$1 AU) structure(s) in the ISM, potentially multiple extreme scattering events (ESEs).

ESEs were originally seen as localized “events” in flux density measurements of compact radio sources (Fiedler et al., 1987; Cognard et al., 1993; Maitia et al., 2003). The first measurements of ESEs observed in DM and scintillation/scattering measurements were performed by Coles et al. (2015), who described ESEs observed towards two pulsars.

This line of sight is particularly unique due to observed “light echoes” in the trailing components of the pulsar’s pulse profile (Michilli et al., 2018) using LOFAR. These additional components were seen to vary with time. Michilli et al. (2018) described several possible mechanisms, with one requiring the lensing of radio emission due to passage near a compact interstellar structure, like those causing ESEs. Such lensing would cause the emission to traverse slightly different paths through the ISM, which could manifest as some combination of dispersive delays and light-travel-time delays. Given the complexity of the observed features, the line of sight toward PSR J2219$+$4754 therefore represents an excellent laboratory to test the predictions of Cordes et al. (2016) to disentangle the properties of the intervening medium.

In this work, using the data on PSR J2219$+$4754 as a test to distinguish between the different physical mechanisms, we carry out analyses to determine whether or not the observations are consistent with lensing events or with the theoretical expectation from a turbulent medium. In §2, we build upon their structure-function analysis to show consistency of the frequency-dependent DM with the theoretical predictions for a Kolmogorov medium in Cordes et al. (2016). In §3, we further examine the reported ESE detection, analyzing the impact of a variety of timing effects on their results. We apply a wide range of analyses discussed in the literature to understanding this specific line of sight, including: (i) determining the probability for the false signature of a gradient to be seen in the DM timeseries; (ii) establishing the importance of the solar wind on the frequency-dependence of DM and testing if a DM gradient is due to the changing Earth-pulsar line of sight; (iii) determining the impact of non-dispersive delays on DM estimation, specifically the amplitude of the delay caused by refraction from a lensing structure; (iv) understanding the systematic bias of pulse profile evolution in both time and frequency on these analyses; (v) accounting for various components of pulse arrival-time uncertainties, demonstrating the importance of high-precision pulsar timing techniques given the sensitivity of the measurements; (vi) developing a method built upon traditional structure function analyses, providing a more robust measure of the expected amplitude of the frequency-dependent DM while further showing inconsistencies with gradients in the DM timeseries; and (vii) implementing this structure function method on the data presented in Donner et al. (2019) to gain additional insights on the line of sight. We found through our analyses that any lensing structure as described by Donner et al. (2019) would affect the various observables in a way that has not been seen, thereby ruling out such a structure. We discuss the impact on precision pulsar-timing experiments in §4 and briefly discuss possible future observations of this type in §5. In the Appendix, we provide a derivation for DM estimation errors in the presence of time-of-arrival (TOA) uncertainties and additional frequency-dependent time delays that is useful for a number of arguments in our work.

2 Structure Function Analysis and Frequency-dependent DM

Structure functions are used in analyses of pulsar DM timeseries, as well as the broader literature on turbulence in general, as a method of constraining the spectral properties of the variations (Cordes & Rickett, 1998). From the structure function, one can derive the amplitude of the electron-density wavenumber spectrum, which directly relates to the size of the fluctuations seen in DM timeseries (Rickett, 1990). Therefore, they can also be used to quantify the magnitude of the differences between DM timeseries measured at two frequencies (Cordes et al., 2016). The variations in the observations presented in Donner et al. (2019) provide a strong data set with which to test for the root-mean-square (rms) fluctuations in the DM differences with both frequency and time.

As part of their analysis, Donner et al. (2019, see their Figure 6) determined the structure function of the $\mathrm{DM}(t)$ timeseries, defined as $D_{\mathrm{DM}}(\tau)=\left<\left[\mathrm{DM}(t+\tau)-\mathrm{DM}(t)\right]^{2}\right>$, where $\tau$ is the time lag separating two observations and the brackets denote the ensemble average. The authors found a power-law structure function with a spectral index consistent with that of Kolmogorov turbulence.The fit was performed over lags $\tau\leq 200$ days. Note that the calculated structure function shows clear evidence of a white noise contribution at lags $\tau\lesssim 30$ days. A white noise contribution to a time series will produce a “plateau” at small lags (e.g., Cordes & Downs, 1985) or for any trends (e.g., linear) in the data (Lam et al., 2016b; Jones et al., 2017) and bias the fit for the spectral index which will itself become a function of $\tau$. For the purposes of this paper, we accept that a Kolmogorov (or near-Kolmogorov) spectrum is an acceptable fit to the DM structure function given the clear overlap with the estimated values shown. We will assume from here that the spectrum is described by a Kolmogorov medium based on the observed consistency with this model in the ISM (e.g., Armstrong et al., 1995).

Following the determination of a Kolmogorov spectrum, the authors then concluded that the cause of the variations may be due to ESEs, which they describe in a context other than localized events. A spectrum steeper than Kolmogorov or discrete structures (or both) are required to produce ESEs and the refractive effects observed in pulsar dynamic spectra (e.g., Roberts & Ables, 1982; Hewish et al., 1985; Cordes & Wolszczan, 1986; Romani et al., 1987). The turbulence in the ISM is often described with an electron-density wavenumber spectrum having a power-law form (see e.g., Armstrong et al., 1995)

[TABLE]

where $q=2\pi/l$ is the wavenumber corresponding to physical scale $l$ and $C_{n}^{2}$ is the spectral coefficient representing the overall amplitude or strength of the turbulence. The low and high wavenumber cutoffs represent outer and inner physical scales, respectively. A Kolmogorov medium corresponds to $\beta=11/3$. Strong refractive effects occur for a density spectral index $\beta>4$ (Cordes et al., 1986), which was not observed for PSR J2219$+$4754. Notably, Donner et al. (2019) referenced Fiedler et al. (1994) and Coles et al. (2015) as examples of ESEs in which larger-scale turbulent structure is identified as the cause of the observed lensing; however, both of those works instead find the causes of their observations to be due to structures consistent with compact canonical ESEs, i.e., smaller-scale lenses as described previously in the literature and above. In the former, the authors did examine the associated large-scale structure of the ISM in the foreground of sources with observed ESEs from Fiedler et al. (1987). In the latter, the outer scale of the ESE was taken to be of order the size scale of the smallest dimension of the lens, which is distinct from the outer scale of the electron-density wavenumber spectrum of the ionized ISM written above, known to be many orders of magnitude larger (Armstrong et al., 1995).

Using the constant value of the amplitude of the wavenumber spectrum Donner et al. (2019) estimated, $C_{n}^{2}=0.9\times 10^{-3}$ m*-20/3*, and a distance estimate of 2.2 kpc from the NE2001 electron density model (Cordes & Lazio, 2002), the corresponding scattering measure is $\mathrm{SM}=\int_{0}^{D_{p}}C_{n}^{2}(z)dz=2.0\times 10^{-3}$ kpc m*-20/3*, an order of magnitude higher than the value predicted by the NE2001 model111We note this measurement (and extrapolated to the estimate of the scintillation timescale discussed later) may then be useful in constraining properties in future electron density models or uncovering interesting turbulence physics along this line of sight.. Following Cordes et al. (2016, see Eq. 12), the rms difference in $\mathrm{DM}$ between two spot frequencies $\nu_{1}$ and $\nu_{2}$, where $\nu_{1}<\nu_{2}$, for a uniform Kolmogorov medium is

[TABLE]

where the frequency ratio $r\equiv\nu_{2}/\nu_{1}$, $F_{11/3}(r)$ is a factor that contains all of the relative frequency dependence (Cordes et al., 2016, see Eq. 11), $G_{11/3}$ is a geometric-only factor where the fiducial value of 0.145 is for a uniform Kolmogorov medium (Cordes et al., 2016, see Table 1), and $D_{p}$ is the pulsar distance. For frequencies at the centers of the two frequency bands from the observations of Donner et al. (2019), we have $\sigma_{\mathrm{DM}(\nu)}(133~{}\mathrm{MHz},169~{}\mathrm{MHz})=0.6\times 10^{-3}$ pc cm*-3*, while for the extremes of the bands, $\sigma_{\mathrm{DM}(\nu)}(118~{}\mathrm{MHz},190~{}\mathrm{MHz})=1.1\times 10^{-3}$ pc cm*-3*. The bottom panel of Figure 5 of Donner et al. (2019) shows rms variations at exactly these levels, suggesting that the observed frequency-dependent DM is consistent with the theoretical predictions of Cordes et al. (2016) for a Kolmogorov medium.

3 Testing the ESE Interpretation

Frequency-dependent DM results from raypath averaging over different volumes of the intervening medium. Refraction by discrete structures along the propagation path can also result in variations in DM timeseries via mis-estimation from non-dispersive delays (Foster & Cordes, 1990; Cordes et al., 2017). These delays can therefore be useful in characterizing such discrete structures (e.g., Lam et al., 2018b). However, care needs to be taken in order to properly account for all effects in timing data.

Following their analysis of the structure function, Donner et al. (2019) proposed an interpretation of the DM variability due to plasma lensing structures. For simplicity, they made an assumption that the variability is due to one to three spherical lenses. Under this assumption, they found the potential parameters of one of the clouds (see their Figure 7). In this section, we will present an alternative explanation from a lens as described by Donner et al. (2019) via characterization of the pulse arrival times in conjunction with the DM timeseries.

Qualitatively, the appearance of the DM time series is similar to that of other pulsars; for example, in the DM time series presented by Jones et al. (2017), one can identify similar features for PSR J0613$-$0200, PSR B1937+21, and potentially other pulsars in that dataset alone.

Evidence of these features being caused by something other than ESEs can be found by considering the expected DM amplitude. Donner et al. (2019) identified the variation in DM between MJD 56950 and 57100 (their Figure 1) as an ESE by virtue of its amplitude $|\delta\mathrm{DM}|\approx 3\times 10^{-3}$ pc cm*-3*. By considering their calculated structure function and relating that to the expected rms variations in the DM (Lam et al., 2016b, Eq. 30), we find

[TABLE]

For $\tau=150$ days, then $\sigma_{\mathrm{DM}}(150~{}\mathrm{days})=0.8\times 10^{-3}$ pc cm*-3*, and $|\delta\mathrm{DM}|$ is much larger than the expected rms variation.

If the DM variation is from a true steep gradient, then we can use an estimate of the rms DM gradient (rather than the rms DM, Lam et al., 2016b),

[TABLE]

and compute the “signal-to-noise” ratio for such a gradient,

[TABLE]

Evaluating for $|d\mathrm{DM}/dt|=2\times 10^{-5}$ pc cm*-3* day*-1*, we find that $R_{d\mathrm{DM}/dt}=3.8$ over the 150-day duration of the purported ESE, suggesting moderate significance. However, while the DM obtained from analyzing the entire frequency range shows a rising linear trend prior to this time period, in the split-frequency DM timeseries, the $<149$ MHz values are consistent with a constant DM within $\sim 1\sigma$, and therefore there is no apparent ingress into the assumed spherically-symmetric lensing material for those frequencies. Therefore, a canonical ESE interpretation seems unlikely, though a potential DM gradient may still be plausible. Note that the egress time period occurs around MJD 57100, right when the DM difference between the timeseries is zero by construction. It is thus unclear if this time period truly marks the end of such an event or not.

We investigated whether this potential DM gradient could be due to the line of sight cutting through different parts of the solar wind. By MJD 56950 moving forward, the suggested time of a peak, the solar elongation is shrinking (for reference, the pulsar’s ecliptic latitude is 52.5∘). Thus, as the pulsar as seen on the sky approaches the Sun we would expect the DM to increase (only slightly due to the high ecliptic latitude) rather than decrease (You et al., 2007b; Jones et al., 2017; Madison et al., 2018). Any solar flare or coronal events that may have occurred during this timespan would also cause an increase rather than a decrease in the DM (Lam et al., 2016b). In addition, we looked at the $C_{n}^{2}$ contribution due to turbulence in the solar wind using the form provided in Spangler et al. (2002),

[TABLE]

where $r$ is the radial distance from the Sun. We integrated the solar wind $C_{n}^{2}$ over the line of sight to determine the SM, shown in Figure 1. With each integration element acting as a thin screen, we determined the frequency-dependent DM error and then combined these to find a total error of

[TABLE]

i.e., the total error is the square root of the reciprocal of the mean of the “weights” ($1/\sigma^{2}$), and we are integrating along the line-of-sight position $s$ for each volume element. We find that while the SM is increasing over the time of the suggested ESE and is of a similar value to the rest of the ISM (fiducial value of 10*-3.5*; Rickett 1990, Cordes & Rickett 1998)—because the material is very close to the Earth in comparison with the pulsar’s distance—the rms frequency-dependent DM is small, well below the measurement uncertainties and even below that of other high-precision timing experiments (see e.g., Jones et al., 2017).

For variations in DM due to the changing ionosphere, in this case because of the yearly modulation from observing the pulsar transitioning between day and night over that time period, the amplitude of the change in DM is at most $\sim$10*-4* pc cm*-3* and therefore is not a significant contributor to the variation here (Lam et al., 2016b).

3.1 Impact of Non-Dispersive Delays

As has been described in the literature, there are multiple frequency-dependent delays that affect pulsar TOAs (Foster & Cordes, 1990; Lam et al., 2016b, 2018b; Shannon & Cordes, 2017). Beyond the traditional dispersive delay due to the integrated electron density ($\propto\nu^{-2}$), there is a geometric delay ($\propto\nu^{-4}$) due to total path length changes and a barycentric-correction delay ($\propto\nu^{-2}$) due to the angle of arrival of the pulsar shifting, i.e., the pulsar’s image appearing from a different direction on the sky.

Following the notation in Lam et al. (2018b), the dispersive delay for a lens of size $L$, electron density $n_{e}$, and dispersion measure $\mathrm{DM}_{l}\sim n_{e}L$, is

[TABLE]

where $\lambda$ is the electromagnetic wavelength, $r_{e}$ is the classical electron radius, and $c$ is the speed of light. The geometric delay is

[TABLE]

where $D_{l}$ and $D_{p}$ are the distance to the lens and pulsar (from the observer), respectively, and $\mathrm{DM}_{l}^{\prime}\sim n_{e}L/(\zeta L)\sim n_{e}/\zeta$ is the DM spatial gradient with the depth-to-length aspect ratio of the lens as $\zeta$. Finally, the barycentric delay is

[TABLE]

where $r_{\oplus}$ is the Earth-Sun distance of 1 AU.

Following Donner et al. (2019), we assume that the cloud is spherical and therefore we set $\zeta=1$. The time for the cloud to pass through the line of sight is $\sim 300$ days. With a proper motion of 22.2 mas/yr (Michilli et al., 2018), the angle subtended on the sky is $\theta=18.2$ mas. The angle sets the physical size of the cloud as $L=\theta D_{l}$. Therefore, the distance of the lens is linearly proportional to the size of the lens (shown by Donner et al. 2019 in the bottom panel of their Figure 7). Given that the DM change from the cloud is simply $n_{e}L$, we have that $n_{e}$ is inversely proportional to the cloud size and inversely proportional to the lens distance (shown by Donner et al. 2019 in the top panel of their Figure 7).

The distance estimate for a possible structure in the ISM causing light “echoes” seen in the pulse profiles for PSR J2219$+$4754  is $D_{l}\approx 1.1$ kpc (Michilli et al., 2018, hereafter Michilli et al. (2018)), a companion paper to the work of Donner et al. (2019). Using the DM change of $\left|\delta\mathrm{DM}\right|\approx 3\times 10^{-3}$ pc cm*-3*, Figure 2 shows the dispersive, geometric, and barycentric delays for such a lensing structure at 1.1 kpc. We see that the dispersive delay dominates over the geometric and barycentric delay. However, when fitting the total frequency-dependent delays observed, the estimated infinite-frequency arrival time is shifted from the true arrival time by $93~{}\mu$s; the estimated DM would then be larger than $\left|\delta\mathrm{DM}\right|$ above by $7.1\times 10^{-3}$ pc cm*-3*, which is unseen in the timeseries unless the baseline DM value is far lower than suggested by Donner et al. (2019).

In Figure 3, we show the perturbation in $t_{\infty}$ and DM (the difference between the “true” DM from the dispersive delay and that from the mis-estimation) after fitting the total delay curve generated by placing a lens at a distance $D_{l}$; again, Figure 2 shows the three frequency-dependent delays along with the total delays when $D_{l}=1.1$ kpc. These perturbations were calculated after fitting each delay curve with the functional form $t_{\nu}=t_{\infty}+K\mathrm{DM}/\nu^{2}$. For reference, the form of these perturbations for two spot frequencies are shown in the Appendix.

If we split the total band into two halves to measure the frequency-dependence of DM, the increases in $\delta\mathrm{DM}$ are then 7.3 and $6.8\times 10^{-3}$ pc cm*-3* for the low and high bands, respectively, or a difference of $5\times 10^{-4}$ pc cm*-3* between the two. This value is approximately half of the difference between the two bands shown by Donner et al. (2019) in their Figure 5. Recall, however, that in this calculation we started with a single DM, and only when estimating the DM in both halves of the band did we recover frequency dependence, solely from the unaccounted for $\propto\nu^{-4}$ geometric delay rather than from a true DM difference. This demonstrates why a complete analysis of the arrival times, such as in Lam et al. (2018b), is crucial to understanding any potential lens/ESEs in the data rather than analyzing only the DM timeseries.

We examined the long-term timing residuals in Michilli et al. (2018, see their Figure 1) for PSR J2219$+$4754 from observations covering 1970 to 2016. Figure 3 shows that the TOA perturbations due to a lens should be of order microseconds or greater depending on the distance, and for a lens at distance 1.1 kpc, we expect a perturbation of $\approx 0.1$ ms if $\propto\nu^{-4}$ delays are not fit for. We replot the residuals as directly shown in Michilli et al. (2018) in Figure 4, starting with the earliest observation epoch in Donner et al. (2019). For reference, decimal year 2015.5 corresponds to the end of the time period in which Donner et al. (2019) stated no baseline DM variations were observed. The blue shaded region denotes the 300-day time period over which the proposed lens occurs. In the middle panel, we show the residuals after a quadratic subtraction so that the scatter is more clearly visible. In the bottom panel, we show the same residuals as in the middle but with an increase in the TOA uncertainties due to jitter ($\approx 10~{}\mu$s) and polarization miscalibration ($\approx 80~{}\mu$s) as described later in §3.3. In red we have injected the systematic effect on the timing given the lens described above, where we simplify input of the offsets as a triangle function starting at 0, decreasing to $-93~{}\mu$s in the middle, and increasing back to 0. Given the scatter in the residuals and the unknown contributions to spin noise, we cannot prove or disprove the presence of a lens at 1.1 kpc. As per Figure 3, a lens closer than 1.1 kpc will show increasingly dramatic dips in the residual timeseries, which would become readily visible.

Michilli et al. (2018) also showed the evolution of the pulsar’s spindown rate. While this curve is also smoothly varying, we see a drop right at the start of 2015 (around MJD 57000), during the time of the proposed ESE. Donner et al. (2019) indicated that a single spin-period derivative (Michilli et al. (2018) showed the equivalent spin-frequency derivative) was fit in their own timing model. The spindown rate shown in Michilli et al. (2018) was not constant over this time period, suggesting that if the spin-frequency derivative was smaller than the average value in the fit by Donner et al. (2019), then the pulsar’s spin was braking more rapidly (i.e. the spin frequency was dropping more rapidly) than expected and the infinite-frequency arrival times should be delayed compared to their model (a linear change would have been absorbed by the shorter-duration fit but we do not see that either). Any fit for DM would therefore be biased by this effect. Visualized in terms of the delay curves in Figure 2, the true arrival time would be delayed, i.e., shifted to the right.

Rather than investigating the impact on the timing residuals, we directly investigated the impact on the two DM timeseries. Given the arguments above, if there was an intervening lens with a true column density equal to $\mathrm{DM}\approx 3\times 10^{-3}$ pc cm*-3* that passed the line of sight, the differences in DM between the two bands could potentially be explained, but not the amplitude change of $\approx 7\times 10^{-3}$ pc cm*-3* in the DM timeseries from the zero value, thereby ruling out a lens with such a column density. It is possible that if a lens with a lower column density were to pass by the line of sight, it would show up with an apparent DM having the appropriate amplitude as estimated by the total time delays. One can search over the phase space of the true column density/DM, size $L$, and distance $D_{l}$ to find the best-fit parameters of a possible lens. As stated previously, for a spherically symmetric lens, we would expect the DM timeseries from both bands to behave symmetrically.

In addition to the three delays described above, general scattering (alternatively, pulse broadening) from a Kolmogorov medium will produce arrival-time delays $\propto\nu^{-4.4}$. Donner et al. (2019) stated that if the frequency-dependent DM is variable in time, then non-$\nu^{-2}$ dispersive delays will vary equally in time. From Cordes et al. (2016), we see that the time- and frequency-dependence of DM arises naturally from ray-path averaging over different volumes of the ISM if there is an effective velocity between the Earth, pulsar, and bulk ISM motions (regarding the effective velocity, see Cordes & Rickett, 1998). Such variations can arise even for a medium with constant $C_{n}^{2}$, which for example can yield a statistically constant scattering timescale even while the DM varies. Changes in pulsar scattering timescales have been observed (e.g., Coles et al., 2015; Levin et al., 2016), which require a change in $C_{n}^{2}$ (or the inner scale of the turbulence, again see Cordes & Rickett, 1998) and therefore the statistics of the DM variations, though the timeseries of each need not be one-to-one correlated. As with the geometric and barycentric delay, a $\propto\nu^{-4.4}$ scattering delay will also cause DM to be incorrectly estimated if not properly modeled. Donner et al. (2019) noted that unmodeled scattering does affect their results but is small enough in amplitude that it does not change their conclusions substantially.

3.2 Impact of Frequency- and Time-Dependent Profile Evolution

The variability of the pulse profile of PSR J2219$+$4754 with time was discussed in Michilli et al. (2018) in the context of light echoes due to propagation effects. Bilous et al. (2016) showed in an observation taken in early 2014 (around MJD 56700) that there was significant frequency dependence in the pulse profile. Using profiles of PSR J2219$+$4754 from Bilous et al. (2016)222Obtained via the European Pulsar Network, http://www.jb.man.ac.uk/research/pulsar/Resources/epn/, we examined the TOA and DM perturbations due to the frequency dependence of the profiles. These four profiles are shown in Figure 5; the trailing components are clearer at lower frequencies than at higher frequencies, suggesting that the shape variations as a function of time as seen in Michilli et al. (2018) are also a function of frequency, at least for some epochs.

Bilous et al. (2016) used their observations to determine a new spin period and DM for each of the pulsars in their census. The initial adjustment of these two parameters maximized the pulse signal-to-noise ratio. Afterward, an average template was generated, which was used to calculate more precise TOAs and perform a subsequent timing analysis to improve their period and DM estimate. They did not account for profile evolution, either intrinsic or from interstellar scattering, in their work but noted it as a potential bias; the profiles used in this analysis are therefore phase-aligned according to their method.

We generated smooth template shapes for each of the four profiles. While pulsar components are often fit with von Mises functions (e.g., Osłowski et al., 2011; Hassall et al., 2012), the circular analogue of a Gaussian function, the automated routine in the PyPulse package (Lam, 2017) to fit multiple such components (see a more thorough description of the procedure conditions in Lam et al. 2019) did not converge for the lower-frequency profiles due to the shapes of the trailing components. We then used a simpler Savitzky-Golay filter (Savitzky & Golay, 1964), implemented in scipy with a moving 11-point cubic polynomial, which produced residuals between the template and the data profile at the level of the rms noise of the off pulse region; slight variations in the parameters did not affect the overall fit significantly. With our four smooth pulse shapes, we used PyPulse to fit the highest-frequency (178 MHz) template to the other three templates in order to compute the arrival-time offsets due to frequency-dependent profile evolution; these were 229.0, 75.7, and 11.2 $\mu$s for the 120, 139, and 159 MHz profiles, respectively. We chose the highest-frequency template as a standard of comparison to minimize the impact of the trailing component shapes as an approximation for the intrinsic pulse shape compared to a pulse shape showing a light echo.

Using these timing perturbations, we then calculated the estimated DM perturbations; see the Appendix for more on DM estimation in the presence of an additional chromatic perturbation, in this case delays due to the pulse shape changes as a function of frequency. Since the method of calculating DM between Bilous et al. (2016) and this analysis are different, we do not expect the absolute DM value to be the same (Lam et al., 2016b). However, using the two lowest or highest frequencies to calculate the DM alone yielded a difference between the two of $1.7\times 10^{-3}$ pc cm*-3*, again of the same amplitude of the variations seen by Donner et al. (2019). This value is an order of magnitude larger than the DM difference they quote between scattered and unscattered profiles. In performing their DM analysis, Donner et al. (2019) accounted for frequency-dependent profile evolution with a frequency-resolved template, which should then mitigate any impact of this profile evolution on the measured DMs. As mentioned previously, since their template was generated from an observation on MJD 57161, we expect the DM at both bands to be equal by construction at that epoch, which is also noted by the authors, though the true absolute DM between the two at that time and over the course of the entire observing span may differ. Nevertheless, the changing profile shapes as a function of time as seen in Donner et al. (2019) and in Michilli et al. (2018) could cause biases in the TOAs and the subsequent DM determinations, especially when coupled with the changes in frequency. Note again that Donner et al. (2019) did state the amplitude of scattering in their data and the biases in their DM; we find that this amplitude does not affect their conclusion that frequency-dependent DM is observed.

3.3 Impact of Arrival-Time Uncertainties

We examined the role of additional components to the TOA uncertainties due to pulse jitter on the DM estimates. Traditionally until the last several years, many analyses considered the TOA uncertainty as arising only from template fitting (e.g., Manchester et al., 2013; Demorest et al., 2013; Desvignes et al., 2016), the process of fitting a smoothed template shape to the data profile. The assumption of matched filtering which underlies this fitting is that the data is a shifted and scaled copy of the template shape. It has been long known that single pulses from pulsars vary stochastically (Craft, 1970), implying that the data shape cannot be an exact copy of a template since the average of a finite number of single pulses will always be slightly different. In general, shape changes due to jitter include contributions both from phase and amplitude variations.

We used the jitter parameter $f_{\rm J}$ defined in Cordes & Downs (1985) for our analysis, defined as the ratio between the single-pulse rms jitter and the equivalent rms (i.e., for a Gaussian pulse, the standard deviation as compared to the full width at half maximum) of the template (Shannon et al., 2014; Lam et al., 2016a). Cordes & Shannon (2010) summarized a number of analyses in the literature, primarily for canonical pulsars, and suggest that $f_{\rm J}\approx 1/3-1/2$. Shannon et al. (2014) found similar values for millisecond pulsars. Lam et al. (2016a, 2019) used a separate jitter parameter that is independent of the pulse shape but found comparable statistics for millisecond pulsars. In many cases, especially for bright canonical pulsars in which single pulses can be detected, jitter is the dominant component of the TOA uncertainty (Lam et al., 2016a).

Again using profiles of PSR J2219$+$4754 from Bilous et al. (2016), calculating the pulse width, and assuming a fiducial value of $f_{\rm J}=1/3$, we found the rms single-pulse jitter for PSR J2219$+$4754 to be 1.2 ms, $\approx 0.2\%$ of pulse phase ($P=0.5385$ s), consistent with some pulsars in Lam et al. (2016a, 2019) though below the average. (Note that millisecond pulsar studies as described above were performed at significantly higher frequencies, and the statistics for millisecond pulsars may be different from that of canonical pulsars.) Since the rms jitter scales as the number of pulses $N_{\rm p}^{-1/2}$, we assumed the pulsar was observed for the median observing time each epoch, 115 minutes, using the LOFAR station most used in Donner et al. (2019), DE605. Given the number of pulses in that time, we find that the TOA uncertainty due to jitter is $\approx 10~{}\mu$s. Using a larger value of the jitter parameter(s) will yield a larger uncertainty. In addition, changes in the integration time will vary the TOA uncertainty; the integration times for DE605 ranged from 2 to 146 minutes, which yield an equivalent rms jitter of 79 $\mu$s to 9 $\mu$s, respectively.

We estimated the value of the TOA uncertainties as follows. Donner et al. (2019) stated a median DM uncertainty of $3.7\times 10^{-5}$ pc cm*-3*. Following the formalism of Lam et al. (2015) and Cordes et al. (2016) of assuming that the (frequency-independent) DM is measured at two spot frequencies and then the infinite-frequency TOA is estimated by removing the dispersive delay, one can calculate the DM difference between the true DM and the estimated DM as

[TABLE]

where $K\approx 4.149\times 10^{9}~{}\mathrm{\mu s~{}MHz^{2}~{}pc^{-1}~{}cm^{3}}$ is the dispersion constant (Lorimer & Kramer, 2012) and $\sigma_{\epsilon_{\nu}}$ is the rms timing uncertainty for frequency $\nu$. (See the Appendix for more on this derivation and its effect on TOA perturbation.) The rms DM uncertainty, $\sigma_{\delta\mathrm{DM}}$, is then $\left<(\delta\mathrm{DM})^{2}\right>^{1/2}$. Assuming that $\sigma_{\epsilon_{\nu}}$ is the same for the two halves of the LOFAR band and that $\sigma_{\delta\mathrm{DM}}=3.7\times 10^{-5}$ pc cm*-3*, we find that $\sigma_{\epsilon_{\nu}}\approx 2.3~{}\mu$s, a factor of four smaller than the jitter uncertainty described above.

If the total TOA uncertainty is then the square root of the quadrature sum of the previous TOA uncertainty of 2.3 $\mu$s and the rms jitter of 10 $\mu$s, we can solve for the corrected rms DM uncertainty and find that $\sigma_{\delta\mathrm{DM}}=1.7\times 10^{-4}$ pc cm*-3*. This is the rms DM uncertainty determined over the whole frequency band. Next we calculate new DM uncertainties assuming that DM is derived from each half of the band separately (making sure to correct the TOA uncertainties for the change in the signal-to-noise ratio by a factor of $\sqrt{2}$ but no change in the jitter as it is roughly frequency-independent over a small frequency range; Lam et al. 2019), as is done by Donner et al. (2019) to find the DM using measurements taken above and below 149 MHz. We split the full band into four and then used the centers of the bottom two frequency channels as our spot frequencies to determine the DM for data taken below 149 MHz. Correspondingly the centers of the top two frequency channels were used to determine the DM for data taken above 149 MHz. We found that $\sigma_{\delta\mathrm{DM}}=3-5\times 10^{-4}$ pc cm*-3*. This range of values is of similar order to those of $\sigma_{\mathrm{DM}(\nu)}$ discussed in §2. Our analysis describing the underestimation of the TOA uncertainties further strengthens the argument that differences in DM between the two frequency bands discussed in Donner et al. (2019) are less significant than were shown. We do note, however, that while the amplitude of these uncertainties adds significantly to the two DM timeseries, this is a white-noise contribution in time and, therefore, cannot explain the systematic offsets between the two that are seen by Donner et al. (2019), i.e., frequency-dependent DM causes non-white-noise shifts in the TOAs and so increasing the error bars does not remove the effect.

We also estimated the contribution to the TOA uncertainty from scintillation noise, also known as the “finite-scintle effect” (Cordes et al., 1990), one of three commonly-analyzed white-noise contributions to the TOA uncertainty on short timescales (Lam et al., 2016a, 2018a). The scintillation timescale $\Delta t_{\mathrm{d}}$ can be found when the structure function of the electromagnetic phase perturbation is equal to 1, or alternatively in terms of the DM structure function (Lam et al., 2016b),

[TABLE]

Given that they observe $D_{\mathrm{DM}}(\tau)=3.1\times 10^{-10}~{}(\mathrm{pc~{}cm}^{-3})^{2}~{}(\tau/\mathrm{day})^{5/3}$, we found that $\Delta t_{\mathrm{d}}=5.7$ s. Due to the time-variability of the trailing components in the profile, we estimated the scattering timescale simply from NE2001 as $\tau_{\mathrm{d}}=0.4$ ms at 150 MHz, or an equivalent scintillation bandwidth of $\Delta\nu_{\mathrm{d}}=0.46$ kHz. The TOA uncertainty from scintillation noise is $\approx\tau_{\mathrm{d}}/\sqrt{n_{\mathrm{ISS}}}$, where $n_{\mathrm{ISS}}$ is the number of scintles (“patches” of increased intensity in the time-frequency plane), given by

[TABLE]

The filling factors $\eta_{t},\eta_{\nu}$ are in the range 0.1-0.3 depending on the properties of the medium, and we have adopted a value of 0.2 for both (Cordes & Shannon, 2010). Given an observation time $T=115$ minutes and a bandwidth $B=71.5$ MHz, we found that contribution of scintillation noise to the TOA uncertainty is 0.4 $\mu$s, much smaller than either the template-fitting component or the jitter component, and therefore it should not factor into the error budget substantially.

Lastly, we looked at the impact of polarization calibration errors on the data (Stinebring et al., 1984). Foster et al. (2015) have shown that polarization leakage can result in significant TOA uncertainties ($\sim$microseconds) for well-timed millisecond pulsars. Gentile et al. (2018) showed that the stability of the Arecibo Observatory system varies quite dramatically between frequencies and epochs and therefore re-calibration must be performed per epoch. Given the errors, we performed an uncertainty analysis expanding upon the calibration procedure of Donner et al. (2019), who followed Noutsos et al. (2015), where PSR J2219$+$4754 itself was used to test the calibration stability of the LOFAR antennas. Noutsos et al. (2015) stated that systematic uncertainties in the polarization leakage are of order 5–10%.

Donner et al. (2019) performed observations at or close to transit, and therefore claim that imperfections in the calibration did not affect the analysis. However, Noutsos et al. (2015) showed that at hour angles far away from transit, the profile differences for PSR J2219$+$4754 can be as large as 30% compared with at transit for the circularly polarized flux. While the observations Donner et al. (2019) performed were closer to transit, we would still expect pulse shape deviations of order a few percent given the observation lengths alone.

Following Cordes et al. (2004), the TOA uncertainty due to pulse shape changes from an incorrect absolute gain calibration is

[TABLE]

where we provide fiducial values as in Lam et al. (2016a, 2018a) for the fractional gain error $\varepsilon$, degree of circular polarization $\pi_{\mathrm{V}}$, and pulse width $W$. Using Eq. 14 as a crude estimator, with a pulse with of $\sim 8$ ms, a circular polarization fraction of 9% (Noutsos et al., 2015), and assuming $\varepsilon\sim 0.1$ (the gain error and polarization leakage are not entirely equivalent quantities but we take the fractional errors above as representative), then the component of the TOA uncertainty is $\sim 72~{}\mu$s, many times larger than the template-fitting errors. If polarization error yields a consistent offset/perturbation in the arrival times, then the net stochastic TOA uncertainty is zero and the frequency-dependent DM analysis should not be affected. However, given analyses such as that of Gentile et al. (2018) regard system stability, we do not expect these calibration errors to be systematic alone. While PSR J2219$+$4754 is not a millisecond pulsar, we see that it is quite likely that polarization calibration errors on the order of microseconds or tens of microseconds are expected (see again Foster et al., 2015), which since comparable to the template-fitting errors, should further be accounted for in dispersive-delay removal.

3.4 Impact of a DM Gradient or Scattering Screen on the Structure Functions

Any additional structure in the DM timeseries beyond that from a turbulent medium will increase the measured structure function (Lam et al., 2016b; Jones et al., 2017). For the following analysis, we will make a distinction between the first-order structure function $D_{\mathrm{DM}}^{(1)}(\tau)$ that we have implicitly discussed previously and the second-order structure function $D_{\mathrm{DM}}^{(2)}(\tau)$, both defined as follows:

[TABLE]

While the first-order structure function will remove any constant term from the timeseries (e.g., the mean), the second-order structure function removes linear terms and can be used to detect discrete changes in any underlying linear trends in the DM timeseries (Lam et al., 2016b). The latter can also be thought of as related to the curvature of the timeseries. For a Kolmogorov medium, it is proportional to the first-order structure function (Lam et al., 2016b, Appendix A). Again, any additional structures beyond a turbulent medium seen in the timeseries will increase these measured quantities.

To understand the impact of a DM gradient or ESE on the structure functions of both orders, we performed simulations as in Lam et al. (2015, 2016b) of red-noise realizations of DM timeseries. We generated 10,000 realizations of 1200-day DM timeseries from a stochastic Kolmogorov medium for which the amplitude was set by the measured structure function $D_{\mathrm{DM}}^{(1)}(\tau)=3.1\times 10^{-10}~{}(\mathrm{pc~{}cm}^{-3})^{2}(\tau/\mathrm{day})^{5/3}$. Next, we added a gradient with slope $2\times 10^{-5}$ pc cm*-3* day*-1* and of length 150 days into each timeseries. To avoid a discontinuity, we added a baseline value of $3\times 10^{-3}$ pc cm*-3* to the higher side of the gradient (i.e., in total we added a slanted step function). Lastly, we added a triangle function representing an ESE into each stochastic realization (separate from the gradient), with slope $2\times 10^{-3}$ pc cm*-3* day*-1* and of length 150 days on either side. Figure 6 (left panel) shows several of these timeseries.

As expected, we see in the first-order structure functions shown in the middle panel that there is an increase in the amplitude as well as a changing slope, though the mean stochastic + gradient structure function is within one standard deviation of the mean stochastic-only structure function. Note that because the length of the actual timeseries is only 1200 days, the range of $D_{\mathrm{DM}}^{(1)}$ at large lag varies significantly from the mean at lags greater than a few hundred days333Note that Donner et al. (2019) calculate the standard deviation of the base-10 logarithm of the structure function, whereas here we calculate the standard deviation of the structure function itself.. Since only a few increments (DM differences) contribute to the averages in the bins at large lags, we do expect a wide variation as we see in Figure 6 and in previous simulations of ours (Lam et al., 2015, 2016b).

If either a density gradient or an ESE is along the line of sight, then the stochastic Kolmogorov component of the measured first-order structure function should be lower in amplitude than previously reported, thereby lowering the measured SM or raising the scintillation timescale. While the slope of the mean first-order structure function for the stochastic + triangle simulations has a steeper slope, the “realization errors” on the structure function for the stochastic-only simulations can lead to both steeper or shallower slopes for a single measured structure function, which must be accounted for when performing these analyses to constrain the consistency of the spectrum with a Kolmogorov medium without bias (Lam et al., 2016b; Jones et al., 2017).

Analyses of the second-order DM structure function have not been performed in the literature. We show the results of the analysis of our simulations in the right panel of Fig 6. We notice two features. First, all three structure functions tend toward the same value at low time lags, which may then allow for a more robust estimate of the scintillation timescale (or alternatively the SM as in Donner et al. 2019) as per Appendix A of Lam et al. (2016b), and therefore the amplitude of the frequency-dependent DM. Second, we see that adding extra components to the DM timeseries on top of the purely stochastic term will produce a more pronounced increase in the value of the second-order structure function at a time lag of $\sim 150$ days, the timescale of the injected structures. There is still some slight overlap in the mean stochastic + gradient versus the mean stochastic-only second-order structure function.

To test our simulations, we used the values of the two DM timeseries estimated for different frequencies as provided in Donner et al. (2019) to calculate both the first- and second-order structure functions for each, shown in Figure 7. The gray regions are the same as in Figure 6 for the stochastic-only simulations, again scaled to the amplitude they estimate of $D_{\mathrm{DM}}^{(1)}(\tau)=3.1\times 10^{-10}~{}(\mathrm{pc~{}cm}^{-3})^{2}(\tau/\mathrm{day})^{5/3}$. Recall that no white noise has been added to our simulations, and therefore there is no flattening of the simulated structure functions as compared with the values estimated from the data. We see significant spread in the values of both structure functions at low lags, due partially to the white noise and partially from the irregular sampling at small lags.

For the first-order structure function, we see no increase in the values for $\tau\gtrsim 100$ days indicative of a discrete structure or gradient as analyzed in our simulations (see again Figure 6). Additionally, we see that the gray regions for the second-order structure function do not well-represent the data for most lags, and we also do not see the equivalent increase at large lags if a discrete structure was present.

By increasing the overall amplitude of the simulated structure functions, we see greater agreement with the data, which can be used to extract information about the medium. As fitting of the structure functions, or even the increments, can be complex as noted previously (Lam et al., 2016b; Jones et al., 2017), we decided instead to visually increase the values only for demonstrative purposes about the utility of these estimates. We increased the simulated structure functions by a factor of 2.5, which are seen in the blue shaded regions in Figure 7. The greater consistency with the data points is easily visible. Since $\sigma_{\mathrm{DM}(\nu)}\propto\mathrm{SM}\propto C_{n}^{2}\propto D_{\mathrm{DM}}^{(1)}$ (see Eq. 2, and also Cordes et al. 2016, Lam et al. 2016b), then a factor of 2.5 increase in the structure function translates to the same increase in $\sigma_{\mathrm{DM}(\nu)}$.

However, the observed rms difference we calculated between the two DM timeseries was $0.7\times 10^{-3}$ pc cm*-3* (see also Figure 6 of Donner et al. 2019), which was in agreement with the prediction for a uniform Kolmogorov medium discussed in §2. Since the rms difference444Note that while typically the structure function is proportional to the variance of a time series at a given lag, in this case we are looking at the proportionality with the rms frequency-dependent DM rather than the rms of the DM timeseries. in the frequency-dependent DM is increasing by a factor of 2.5 due to an increase in the estimated $D_{\mathrm{DM}}^{(1)}$, or alternatively the SM, then we find disagreement with the initial assumption of the line of sight being a uniform Kolmogorov medium.

One common model for describing the ISM along the line of sight is not as uniform but with a thin scattering screen model since the mathematics is simpler and the resulting quantities in a calculation depend only on geometric terms (Cordes et al., 1986; Cordes & Rickett, 1998); the geometry of the screen is independent from the spectral index (i.e., Kolmogorov or otherwise). Such a model can also physically describe regions of larger-scale overdense structure in the ISM. Rather than assuming that the medium along the line of sight is entirely uniform, we can quantify what effect a factor of 2.5 increase in the structure function implies if we assume the thin-screen model. In Eq. 2, we provided the geometric factor $G_{11/3}=0.145$ for a uniform medium via Cordes et al. (2016). If instead we assume a screen at distance $D_{s}$ from the pulsar, at distance $D$, then the geometric factor becomes $G_{11/3}=\left[x(1-x)\right]^{5/6}$, where $x=D_{s}/D$ is the screen’s fractional distance between the pulsar and us. Note that the function is symmetric about the halfway ($x=0.5$) point. Assuming a screen describes the line of sight allows us to simultaneously increase the SM while decreasing the value of $G_{11/3}$ to keep $\sigma_{\mathrm{DM}(\nu)}$ fixed at its observed value. Taking the observed DM difference rms and using Eq. 2 to instead determine $G_{11/3}$ and therefore $x$, we found that $x=0.045$ or $x=0.955$, i.e., the screen is at a distance 0.1 kpc from either the pulsar or the Earth. In the case of the latter, this nominal distance lies a factor of a few beyond that of the Local Bubble (Frisch et al., 2011). An $x\approx 0.98$ would place the screen in the rough appropriate location ($\sim 40$ pc), which in turn decreases $G_{11/3}$ or increases the structure functions by an additional factor of 2, which starts to come in tension with the values shown in Figure 7. However, since the pulsar’s distance is uncertain due to estimation with the NE2001 model, then the Local Bubble could still be a plausible location of a screen.

Since $G_{11/3}$ and therefore $\sigma_{\mathrm{DM}(\nu)}$ increases as $x$ tends toward the midpoint of the line of sight, we can rule out any screen structure much closer than the endpoints of the line of sight since the estimated structure functions would be in greater tension with those observed. Our method demonstrates that we can simultaneously use the variations in the DM timeseries taken at different frequencies along with the variations of the differences between the two DM timeseries to constrain the location of an assumed screen quite strongly. Future observations of this and other pulsars taken at widely separated frequencies with high cadence may be able to jointly model even higher-order structure functions for even further constraining power.

Due to irregular sampling in DM timeseries, common for many pulsar observations, it may be preferred to analyze the individual second-order increments ($\mathrm{DM}(t+\tau)-2\mathrm{DM}(t)+\mathrm{DM}(t-\tau)$, which are squared and averaged over to obtain $D_{\mathrm{DM}}^{(2)}$) as was performed for PSR B1534+12 in Lam et al. (2016b) than the second-order structure function as performed here. However, the high-cadence of observations in Donner et al. (2019) combined with the duration has allowed us to directly calculate the second-order structure function rather than simply the increments. Continued observations, especially with a dense high-cadence program, could help constrain the values of both structure functions at smaller lags.

Note that the large transverse velocity means that the motion of the pulsar across the sky is fairly straight (e.g., the parallax motion is small, see trajectory plots in Jones et al. 2017). The quasi-periodic variations in the DM timeseries that deviate from a purely power-law spectrum, seen in Madison et al. (2018) due to the line-of-sight crossing correlated spatial DM fluctuations, should therefore be small and not impact the power spectrum/structure function significantly. Put another way, while the assumption of the line of sight crossing independent DM fluctuations is broken and we expect short-term rapid variations in the DM that might impact the measured structure function, we expect this change to be negligible for this pulsar.

In conclusion, the consistency of the structure functions with a Kolmogorov model over the range of time and spatial scales observed suggests that the medium is turbulent down to 10s of days or $\sim$1 AU (for an assumed lens at distance 1.1 kpc). As discussed, below this scale the structure function appears dominated by white-noise fluctuations given the cadence of observations, which causes the observed structure function/spectrum to become shallower. The $L=20$ AU lens size (equivalent to the 300 day timescale) should be seen as a steepening in the structure function as shown in Figure 6. In addition, any truly stochastic process with a spectrum steeper than that of a Kolmogorov medium, which could then produce strong refractive effects (Cordes et al., 1986), is also unseen in the structure functions in Figure 7.

4 Impact on Precision Timing

Temporal and frequency-dependent DM variations will have significant impacts on high-precision pulsar timing experiments, including the efforts to detect low-frequency gravitational waves. While these topics have been discussed in the literature (see e.g., You et al. 2007a; Keith et al. 2013; Lee et al. 2014; Lam et al. 2016b; Cordes et al. 2016), we will take the analysis of Donner et al. (2019) further and discuss the impact in the context of our accounting of the various noise processes in the real data.

Donner et al. (2019) estimated the limiting TOA precision by considering observations at 1.4 GHz with a relatively limited 250 MHz bandwidth (20% bandwidth). However, current pulsar backend systems can process up to 800 MHz of bandwidth (see e.g., the Appendix of Dolch et al., 2014) and techniques have been developed to process even larger bandwidths (e.g., Dunning et al., 2015) while simultaneously compensating for the DM (Pennucci et al., 2014). More generally, a common practice is to conduct simultaneous or near-simultaneous observations at two frequencies (e.g., the North American Nanohertz Observatory for Gravitational Waves, NANOGrav, uses combinations of 0.8 GHz and 1.4 GHz for some pulsars and 1.4 GHz to 2.3 GHz for some others; Arzoumanian et al. 2018), from which even higher precision DM estimates can be obtained.

When only purely dispersive delays factor into TOAs measured at two spot frequencies, the rms TOA is given by (see the Appendix but also Lee et al. 2014, Eqn. 12; Lam et al. 2018a, Table 1),

[TABLE]

where again $\sigma_{\epsilon_{\nu}}$ is the rms timing uncertainty for frequency $\nu$. This timing uncertainty includes all sources of white noise that affect the TOA estimation, such as from radiometer noise or jitter (see e.g., Lam et al., 2018a). A critial aspect of Eq. 16 is that the rms TOA is dependent only on $r$ and not a specific frequency. It does not matter whether the timing measurements are acquired at low or high frequencies so long as the individual frequency-channel TOA uncertainties are the same for a given value of $r$.

As an illustration, taking $r\sim 2$ and setting $\sigma_{\epsilon_{\nu}}=\sigma_{\epsilon_{\nu_{1}}}=\sigma_{\epsilon_{\nu_{2}}}$, then $\sigma_{\delta t_{\infty}}=1.4\sigma_{\epsilon_{\nu}}$. Therefore, in cases where the median TOA uncertainty from finite pulse signal-to-noise ratio is small (i.e., well below 1 $\mu$s), such as reported per-pulsar and per-band in the NANOGrav 11-Year Data Set (Arzoumanian et al., 2018), then the requirement for sub-microsecond precision timing is met. Several other sources of error, such as from pulse phase jitter, are known to be much smaller than this limit (Lam et al., 2019). Therefore, using low-frequency timing data to increase $r$ will quantitatively improve the timing of many pulsars as long as the TOA uncertainties for pulses at those frequencies are low enough and unmitigated ISM effects are small, e.g., typically for pulsars with lower DM values (Lam et al., 2018a).

As an extreme case, we considered $\sigma_{\epsilon_{\nu_{2}}}\to 0$. In this case, the overall timing precision will be dominated by the precision at the lower frequency, $\sigma_{\delta t_{\infty}}\approx\sigma_{\epsilon_{\nu_{1}}}/(r^{2}-1)$. For a sufficiently large value of $r$, the timing precision due to DM uncertainty or variations could be made negligible. As a specific example, with relevance to the data presented by Donner et al. (2019), we considered $r\approx 10$, equivalent to $\nu_{1}\approx 140$ MHz and $\nu_{2}\approx 1400$ MHz. In this hypothetical example, a timing precision $\sigma_{\delta t_{\infty}}\approx 10$ ns could be obtained, if $\sigma_{\epsilon_{\nu_{1}}}\approx 1\,\mu\mathrm{s}$. We chose this illustration because a timing precision of order 10 ns is comparable to the expected precision required for the study of gravitational waves.

This timing improvement of course neglects the uncertainties due to frequency-dependent DM as discussed by Cordes et al. (2016). The amplitude of the differences in DM between frequencies also depends on $r$ (as well as the specific frequency choices) but the net effect on the overall timing precision can be quantified and then built into noise models, e.g., via covariance matrices such as constructed in Lam et al. (2018a). Cordes et al. (2016, see Figures 7 and 8) show that for pulsars with $\mathrm{DM}\lesssim 30$ pc cm*-3*, combining data from the 100 MHz and 2 GHz regimes will yield TOA errors that still meet the requirement of sub-microsecond precision, especially as current and future telescopes come online, drastically reducing the template-fitting errors from finite pulse signal-to-noise ratio. Many pulsars used in precision-timing experiments have DM values in this range (Verbiest et al., 2016). Cordes et al. (2016) also discussed fitting a wide range of frequencies versus only two spot values; as expected, the increase in frequency coverage can in many cases improve the timing precision over the case where only two spot frequencies are used, which is no longer the case for many modern precision-pulsar-timing experiments (Verbiest et al., 2016).

Unmitigated chromatic delays will also add to the TOA uncertainties (Lam et al., 2018a, Table 1) and cause frequency-dependent excess noise in the timing residuals (Lentati et al., 2016; Lam et al., 2017). However, even for the pulsar with the highest DM millisecond pulsar used in precision timing, PSR J1903+0327 ($\mathrm{DM}\approx 300$ pc cm*-3*), observed between roughly 1.1 and 2.5 GHz, the long-term rms residual is $4~{}\mu$s, though a significant portion of that rms is again from frequency-dependent excess noise (Lam et al., 2017), which future timing methodologies might be able to mitigate partially (Shannon & Cordes, 2017).

Lastly, we expand on the discussion in Donner et al. (2019) that at higher frequencies, longer-term DM variations are of particular importance to take into account for timing data. Given that the frequency-dependent DM comes from differences in ray-path averaging whereas the trends in DM come from the relative Earth-pulsar motion along the line of sight, we expect the longer-term DM variations to track each other between frequencies, regardless of frequency; Donner et al. (2019) also agree with this given their analysis. However, in terms of the overall timing, it is the lower frequencies that are impacted much more heavily since the dispersive delays are weighted by $\nu^{-2}$. At frequencies much higher than typically used in precision timing, the dispersive delay becomes small and thus any changes in the dispersive delay are also small; as an extreme, X-ray pulsar data do not require DM corrections of any kind.

As an additional consideration, it is important to remember that short-term DM variations that are improperly corrected for can contribute heavily to the overall rms timing, and may not contribute as simple white noise (e.g., Lam et al., 2015). Many pulsars in precision-timing experiments show very rapid timescales for DM to vary (Jones et al., 2017), including from the solar wind (You et al. 2007b; Madison et al. 2018; see also Howard et al. 2016 for the study of a coronal mass ejection with a slow-period pulsar) or structures in the ISM (e.g., Fonseca et al., 2014; Coles et al., 2015; Lam et al., 2018b).

5 Discussion

In this work, we have described the requirements necessary for observational tests separating the effects of frequency-dependent DM from refractive lensing. The observations shown in Donner et al. (2019) provide an excellent test of frequency-dependent DM as laid out by Cordes et al. (2016). While the theoretical treatment considers the ability to perfectly measure DM at a given spot frequency, going forward, analyses of observations of the kind reported in Donner et al. (2019) must account for the wide range of different TOA uncertainty components and ISM propagation effects presented in this work. We showed here that even for a canonical pulsar versus a millisecond pulsar, the recent work done in precision-timing experiments has become relevant given the timing quality of the pulsars and the instrumentation used to access to new types of observations, e.g., low-frequency observations via LOFAR as described here. And while canonical pulsars may not provide constraints on the same tests of fundamental physics as millisecond pulsars, their use in studying variations in the ionized ISM along many different lines of sight (e.g., Petroff et al., 2013) will be unparalleled given the greater population of them over millisecond pulsars, especially if observations covering a large frequency ratio can be leveraged.

While we have focused our analyses on the plasma lens as described by Donner et al. (2019), we have not discussed the potential causes of the light echoes as seen by Michilli et al. (2018). They described the similarities between their pulse profiles and those seen in PSR B0531+21 (the Crab Pulsar), though those profile variations are attributed to structure in the local environment in the nebula surrounding the pulsar (Backer et al., 2000; Lyne et al., 2001). While an interstellar lens would produce a negligible effect if very close to the pulsar, as shown in Figure 3, the model proposed by Backer et al. (2000) involves the pulsar traversing near material (a prism geometry) close enough that emission at different frequencies passes through different electron content at different times (a true dispersive delay, whereas they argue that the refractive geometric delay will be significantly smaller) while that of Lyne et al. (2001) involves reflections of the images. Note that the fact that the Crab Pulsar has both a main pulse and interpulse means that propagation variations will affect both in the same way, making it easier to disentangle from intrinsic profile shape variations. Both of the proposed mechanisms could be examined in more detail with respect to the pulse profiles shown in Michilli et al. (2018). If in a combined analysis the variability is also tied to the DM variations shown in Donner et al. (2019), then the total data set will provide an excellent probe of the material local to PSR J2219$+$4754, though such variations will then need to be disentangled from the variability expected from a turbulent medium with which we have shown consistency.

Identifying ESEs or other “ISM events” seen in pulsar timing data in near-real time will allow for more intensive follow-up observations, including a higher cadence of observations over many frequencies and using different observatories worldwide, especially if dynamic spectra with resolved scintles can be obtained (Hewish, 1980; Stinebring et al., 2001). The characteristic timescale and bandwidth can be used to constrain the location of a lensing structure (Cordes & Rickett, 1998) while the drift rate (“rotation”) of the scintles provides the component of the refractive angle in the direction of the pulsar’s motion (Hewish et al., 1985; Cordes et al., 1986), providing partial information on the geometric time delay. Interferometric observations can help constrain the changing position and sizes of the pulsar image (Blandford & Narayan, 1985), or possibly multiple images (Cordes & Wolszczan, 1986; Cordes et al., 2017), providing additional constraints on the lensing geometry over the line of sight. Any such additional observations will allow us to resolve small-scale structure in the ISM and probe the Galactic population of these lenses.

The NANOGrav Project receives support from NSF Physics Frontiers Center award number 1430284. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Part of this research has made use of the database of published pulse profiles maintained by the European Pulsar Network, available at: http://www.jb.man.ac.uk/research/pulsar/Resources/epn/.

Appendix A DM Estimation with Additional Chromatic Errors

As a useful reference, here we will describe DM estimation from observations taken at two spot frequencies with additional chromatic errors following the formalism of Lam et al. (2015) and Cordes et al. (2016). We can write the TOA at a particular frequency $\nu$ as the infinite-frequency arrival time plus the dispersive delay term. For this calculation, we will also include measurement errors $\epsilon_{\nu}$ and a chromatic (frequency-dependent) timing perturbation $t_{C,\nu}$, such that

[TABLE]

Here, $K\approx 4.149\times 10^{9}~{}\mathrm{\mu s~{}MHz^{2}~{}pc^{-1}~{}cm^{3}}$ is the dispersion constant in observationally convenient units (Lorimer & Kramer, 2012). We estimate the DM by taking TOAs at two frequencies $\nu_{1}$ and $\nu_{2}$ and calculating

[TABLE]

As in the main text, we will define the frequency ratio $r\equiv\nu_{2}/\nu_{1}$ with $\nu_{1}<\nu_{2}$. The estimated infinite-frequency arrival time can then be written in one of two ways as

[TABLE]

We will now solve for the DM difference. Substituting the measured TOAs $t_{\nu}$ into the equation for $\widehat{\mathrm{DM}}$ and subtracting from the true DM, we have

[TABLE]

The TOA perturbation will be

[TABLE]

When the chromatic offsets are zero, we arrive at simply

[TABLE]

which agrees with Eq. 21 in Cordes et al. (2016) assuming the frequency-dependent DM term is zero.

Eq. A6 provides the timing offset but one must consider the TOA uncertainty, $\sigma_{\delta t_{\infty}}$, from the variance

[TABLE]

We have assumed here that $\left<\delta t_{\infty}\right>=0$, which will be true if the errors $\epsilon_{\nu}$ are Gaussian distributed. If they are, and with variance $\sigma_{\epsilon_{\nu}}^{2}$, then the sum of the two terms in the numerator of Eq. A7 will be Gaussian distributed, which when squared will then be chi-squared distributed. Taking the expected value of the resultant quantity yields

[TABLE]

If the TOA uncertainties are the same at both frequencies such that $\sigma_{\epsilon_{\nu_{1}}}=\sigma_{\epsilon_{\nu_{2}}}=\sigma_{\epsilon_{\nu}}$, then we have

[TABLE]