The NANOGrav 12.5-Year Data Set: Dispersion Measure Mis-Estimation with Varying Bandwidths

As a radio signal is emitted at a pulsar and as it propagates through the ISM, it will encounter various frequency-dependent delays. As a result, pulses with frequency $\nu$ will arrive at Earth’s position at a time $t_{\nu}$ that is delayed with respect to the expected time $t_{\infty}$ for a signal of infinite frequency (i.e., assuming no chromatic effects). In this section we will describe the various frequency-dependent effects on the TOAs responsible for this delay, and how they are accounted for in our timing models.

For each pulsar, TOAs are calculated for all frequency channels recorded with a given receiver using a single standard template profile. However, pulse shapes vary with frequency (Kramer et al., 1998; Pennucci et al., 2014) even with no intervening ISM, so when compared against a single-frequency template this introduces small systematic frequency-dependent perturbations in the TOAs. These changes are modeled as polynomials in log-frequency, described by the $\mathrm{FD}_{k}$ parameters as (NANOGrav Collaboration et al., 2015)

While propagating through the ionized plasma, the pulsar signal encounters ionized plasma and electron density variations along the way, resulting in a frequency-dependent index of refraction. As a result, the radiation will suffer a first-order chromatic delay – the interstellar dispersion – which is the largest frequency-dependent effect due to the ISM. For a cold, unmagnetized plasma, a pulse observed at frequency $\nu$ is delayed compared to one at infinite frequency by an amount $t_{\mathrm{DM}}=K\times\mathrm{DM}/\nu^{2}$, where the dispersion measure ($\mathrm{DM}$) is the line-of-sight (LOS) integral of the free electron along the line of sight to a pulsar. The DM can be quantified as $\mathrm{DM}=\int_{0}^{d}n_{e}(l)dl$, where $n_{e}$ is the free electron density along the line of sight $l$, and $d$ is the pulsar distance. The dispersion constant is given by $K=e^{2}/2\pi m_{e}c\simeq 4.149~{}\mathrm{ms}~{}\mathrm{GHz}^{2}\mathrm{pc}^{-1}~{}\mathrm{cm}^{3}$.

Turbulent and bulk motions within the ISM, solar wind, differences in the relative velocity of the pulsar and the Earth, and stochastic variations from pulsar motion can cause the line of sight to sample electron-density fluctuations on a variety of scales (Lam et al., 2016; Cordes & Rickett, 1998; Phillips & Wolszczan, 1991). The result is a DM that varies with time, changing on timescales of hours to years. To model these short-scale variations, we use a stepwise model for variation in DM, in which DM is allowed to have independently varying values in time intervals. The time intervals range in length from $0.5$ to $15$ days, depending on the telescope and instrumentation. The offset from the globally fixed fiducial $\mathrm{DM}$ value is given by the epoch-dependent $\mathrm{DMX}$ parameters (e.g., Jones et al., 2017; Shapiro-Albert et al., 2021) The ensuing timing correction is given by $t_{\mathrm{DMX}_{i}}=K\times\mathrm{DMX}_{i}/\nu^{2}$, where $\mathrm{DMX}_{i}$ is the correction corresponding to the observing epoch $i$.

In addition to the $\nu^{-2}$ offsets due to dispersion²²2Angle-of-arrival variations from refraction will also yield $\nu^{-2}$ delays (Foster & Cordes, 1990) but since they are entirely covariant with the dispersive correction, they are absorbed in the fit and we will ignore them further., there are also a variety of frequency-dependent effects for which the perturbations scale with radio frequency obeying a $\nu^{-\alpha}$ power law (e.g., Lam et al., 2016). This includes geometric perturbations such as delays due to incorrectly referencing the arrival time of the pulse at the solar system barycenter ($\alpha=2$; Foster & Cordes, 1990). It also comprises pulse scattering, associated with the variable path length due to refraction which results in the signal reaching the observer along different geometrical paths (e.g., Bansal et al., 2019; Lewandowski et al., 2013). As a result, the pulse will arrive at the observer over a finite interval and it will be enveloped by a pulse broadening function. The thin-screen scattering model (e.g., Shannon & Cordes, 2017) considers an isotropic homogeneous turbulent medium with a Gaussian ($\alpha=4$; Lang, 1971) or a Kolmogorov ($\alpha=4.4$; Romani et al., 1986) distribution of inhomogeneities. Finally, interstellar scintillation will introduce a random component in the TOA delay whose variance is strongly frequency dependent. A simple model for the single-epoch TOA delay introduced by these chromatic effects is:

By incorporating all these effects into our timing model, we quantify the time of arrival (TOA) as a function of frequency $\nu$ as

In addition to these delays, there are non-power-law frequency-dependent timing effects that can modify the pulse arrival times. As previously stated, the DM is defined as the line-of-sight integral of the electron density. However, the line of sight can change as a function of frequency as a result of ray paths at different frequencies covering different volumes through the ISM (Cordes et al., 2016). Therefore, DM itself is frequency-dependent, i.e., $t_{\rm DM}=K\times\mathrm{DM}(\nu)/\nu^{2}$, and this dependence will introduce timing errors that cannot be mitigated solely by increasing the observing bandwidth (Lam et al., 2018). Instead, in this analysis we fit only for a constant DM over the observation, and aim to quantify the mis-estimation in the DM that are introduced by other chromatic effects in the observation altering the DM fit. Other non-power-law frequency-dependent timing effects may be present at small amplitudes along specific lines of sight as well; for example, refraction through plasma-lens-like structures in the ISM is manifested in distinctive DM and geometric path variations that result in such delays (see e.g., Eq. 17 in Cordes et al., 2017). These effects result in higher-order corrections in the measured TOAs that we can neglect for the precision required in this work.

In Fig. 1 we present a qualitative representation of the timing residuals (differences between the observed times of arrival and the predictions from the timing model) that would be obtained by applying a simple timing model to three sets of TOAs that cover either GASP’s bandwidth, GUPPI’s bandwidth, or the full range of frequencies from $0.7~{}\mathrm{GHz}$ to $1.9~{}\mathrm{GHz}$. For each backend, an artificial set of frequency-dispersed TOAs was fabricated by simplifying Eq. 3 as $t_{\mathrm{obs}}(\nu)=\mathrm{K\times DM}/\nu^{2}+t_{\mathrm{C}}/\nu^{4}$, using frequencies $\nu$ in the backend’s bandwidth. We choose $\mathrm{K\times DM}=1~{}\mathrm{GHz}$ for the dispersion coefficient and $t_{\mathrm{C}}=0.01~{}\mathrm{GHz}$ for the chromatic coefficient. The “observed” TOAs were used to fit a simple timing model $t_{\mathrm{pred}}(\nu)=a+b/\nu^{2}$ that only accounts for $\nu^{-2}$-delays in order to simulate the effects of having other chromatic effects being absorbed into the DM fit. This function is then evaluated at the same frequencies to obtain “predicted” TOAs. By subtracting both sets of TOAs we obtain the residuals $\Delta t=t_{\mathrm{obs}}-t_{\mathrm{pred}}$ presented in the figure. We observe that the residuals vary significantly depending on the bandwidth of the backend that was used to take the observations. Effectively, a larger bandwidth allows for a larger sampling of the frequency space over which to model the dispersion. The GUPPI backend combines non-simultaneous observations around two frequency bands, one centered near 820 MHz and another near 1400 MHz, to cover most of the bandwidth from 0.7 GHz to 1.9 GHz. Even then, the resulting timing residuals differ from the estimation we would expect if the receiver covered the full band. As a result, we expect a bias in our measurements even when more advanced backends are used.

We expect that the effects of DM mis-estimations on timing residuals should be more clearly discernible in highly dispersed pulsars. Therefore, the main focus of our work are the observations of PSR J1643$-$1224, a 4.62 ms-period, high spin-down ($dP/dt=1.85\times 10^{-20}$) pulsar in a 147-day binary orbit with a white dwarf companion. This pulsar is of particular interest because lies behind the Hii region Sh 2$-$27, which has an inferred diameter of $0.034$ kpc assuming spherical symmetry (Harvey-Smith et al., 2011). As a result, its pulses suffer high interstellar dispersion ($\mathrm{DM}=62.3~{}\mathrm{pc~{}cm^{-3}}$). Furthermore, PSR J1643$-$1224 has been shown to have significant scattering and profile shape variations (Shannon et al., 2016a; Lentati et al., 2017).

PSR J1643$-$1224 has been observed by NANOGrav for over $12.7$ yr using the GBT with a nearly monthly cadence (Alam et al., 2021a). During this time, its timing residuals have been reported to exhibit significant red noise (red noise spectral amplitude at $f=1~{}\mathrm{yr}^{-1}$, $A_{\mathrm{red}}=1.619~{}\mathrm{\mu s}~{}\mathrm{yr}^{1/2}$) , which may include contributions from unmodeled interstellar-medium propagation effects (Arzoumanian et al., 2018). A dip in the timing residuals between 2015 February 21 and March 7 was found by Shannon et al. (2016b), which is associated with a sudden change of pulse profile. When not modeled, this dip affects upper limits on the stochastic gravitational-wave background (GWB), so it has been included in subsequent timing models.

In Fig. 2 we present the timing residuals for PSR J1643$-$1224 using NANOGrav’s 12.5-yr data release. We observe in the middle panel that during the time period when GASP was the predominant data acquisition backend, the epoch-averaged residuals are highly correlated and track each other. Most importantly, this trend is not present in the GUPPI data set. A possible explanation for this correlation are unaccounted offsets caused by DM mis-estimation in GASP’s narrower bandwidth. In addition, we notice that the residuals in GUPPI’s data set, especially those in the $800$ MHz band, generally exhibit larger uncertainties and a bigger spread than those taken with GASP in the same frequency band. These features are in agreement with the behavior predicted by Fig. 1 for a timing model that is affected by DM mis-estimations due to other chromatic effects: when we calculate residuals sampling a wider frequency range we expect some frequencies to be more heavily affected by such mis-estimations and, therefore, to exhibit larger residuals and uncertainties (orange curve in Fig. 1), but those same frequencies might not be present when we sample narrow frequency ranges (green curve). These results demonstrate the need for an in-depth analysis of the effects of interstellar dispersion and the mis-estimations in its modeling on timing residuals.

For the purposes of this work, we started with a set of $11592$ observations taken with the GUPPI backend as early as 2010 until late 2017, covering a time baseline of $\sim 7.5$ yr. We have separated these observations into two data sets, the original and a modified version:

• •

  The full set of GUPPI broadband observations, which covers a bandwidth of $1157.49$ MHz. For practical purposes, we consider that the timing solution obtained from this data set provides the “best estimation” DM parameters that later on will be compared against those resulting from the narrowband approximation.

• •

  In addition, we created an artificial set of narrowband observations by using Astropy (Astropy Collaboration et al., 2022) to filter out all the GUPPI broadband observations outside the frequency ranges of GASP’s two receivers: Rcvr_800 ($792$-$884$ MHz) and Rcvr1_2 ($1340$-$1432$ MHz). In doing so, we emulate the data set that would have resulted from continuing to use GASP’s bandwidth during the same time period.

These two sets of observations, alongside with the corresponding frequency ranges, are presented in Fig. 3.

Using the time windows that are specified by the DMX parameters (see Sec. 2.1), we have grouped the observations into different time windows, each of them with observations up to 6 days apart. In order to accurately measure the pulsar’s dispersion properties on monthly timescales, and to account for any evolution in these frequency-dependent properties over time, we only considered windows that contain observations in both frequency bands and discarded all observations in windows that do not cover both bands. As a result, the set of broadband observations reduces to 9604 and the set of narrowband observations to 1511. Table 2 summarizes the data sets used for this work.

To isolate the biases introduced in the timing parameters by fitting them using narrowband observations, we created a simplified timing model that includes corrections for long-term interstellar dispersion and frequency-dependent profile evolution (the $K\times\mathrm{DM}/\nu^{2}$ and $t_{\mathrm{PE,\nu}}$ terms in Eq. 3, respectively) but ignores the short-scale variations in the DM that are normally corrected by the $\mathrm{DMX}$ parameters, as well as additional chromatic variations (the $K\times\mathrm{DMX}/\nu^{2}$ and $t_{\mathrm{C,\nu}}$ terms). This simplified timing model is loaded into PINT (Luo et al., 2021), a Python package used for pulsar timing and related activities. We then used this model to generate predicted TOAs at each of the observing dates and frequencies, which are then subtracted from the observed TOAs to produce timing residuals. As a result of ignoring the short-scale DM variations and additional chromatic delays, the residuals within each of the time will display exhibit a curve close to $\nu^{-2}$ as presented in Fig. 4. We can then model the effects of the chromatic delays assuming the residuals can be described by

where $r_{\mathrm{\infty},i}$ is the residual expected at epoch $i$ if no chromatic effects were present, $r_{\mathrm{2},i}=K\times\mathrm{DMX}_{i}$ is the correction for short-scale DM variations, $r_{\mathrm{\alpha},i}$ is the correction due to a scattering-based delayed, corresponding to $\alpha=4.4$. We used LMFIT (Newville et al., 2014) to fit Eq. 4 to the residuals within each time set, thereby obtaining best-fit values and uncertainties for $r_{\mathrm{\infty}}$, $r_{\mathrm{2}}$, and $r_{\mathrm{\alpha}}$. Such a fit is data set-dependent and any unaccounted biases in the observations, such as DM mis-estimations, will result in biased fitted parameters.

This process was repeated using the narrowband and the broadband data sets separately (see Sec. 3). As a result, for each time window at epoch $i$, we obtained two sets of $\{r_{\infty,i},r_{\mathrm{2},i},r_{\mathrm{\alpha},i}\}=\{r_{k,i}\}_{k=\infty,2,\alpha}$ values with their corresponding fitting errors $\varepsilon_{r_{k,i}}$: one set resulting from using the broadband (BB) data set and another from using the narrowband (NB) one. For each $r_{k,i}$ we computed the parameter residuals

In the left panels of Fig. 5(a) we plot the resulting values of $\Delta r_{k,i}-\overline{\Delta r_{k}}$ with their errors as a function of the MJD at the middle of the window. If the frequency bandwidth played no role in modeling interstellar dispersion, we would expect to obtain $\Delta r_{k,i}=0$ for all the parameters $k$ and all epochs $i$. However, the fact that we find non-zero deviations between the broadband and narrowband sets of fitted parameters is indicative that using narrower frequency bandwidth leads to mis-estimations in modeling these parameters for PSR J1643$-$1224. In order to quantify such mis-estimations, we calculate the standard deviation for each set of parameter differences, which we will hereby use as a measurement of the error introduced in the residuals due to incomplete modeling of interstellar dispersion. The largest offset, corresponding to $r_{\infty}$ is $\sigma_{\Delta r_{\infty}}=22.2~{}\mu$s.

For each fitted parameter, in the right panels of Fig. 5(a) we plot histograms of the residuals $\Delta r_{k,i}-\overline{\Delta r_{k}}$ divided by their corresponding fitting errors $\varepsilon_{\Delta r_{k,i}}$. We observe that for $r_{2}$ and $r_{\alpha}$ the residuals are consistent with zero, which is indicative of a small offset introduced in the fit of these parameters when using narrowband data. However, the histogram for $r_{\infty}$ reveals a skewed distribution (significantly more noticeable for J$1643-1224$) which suggests a systematic offset in the estimations of the infinite frequency arrival times.

Next, in order to study variations in the behaviour of the ISM, we analyze whether the chromatic delay exhibits correlations in time. In that case, we expect such a pattern to be revealed in the autocovariance function (ACF) of this quantity. Therefore, we compute the ACF of $\Delta r_{k,i}-\overline{\Delta r_{k}}$ (see Eq.5) for each parameter $r_{k}$ between consecutive time windows. This process involves correlating a signal $f(t)$ with a delayed copy $f(t+\tau)$ of itself as a function of delay the $\tau$. However, since we have a discrete and irregularly sampled signal, we calculated a binned ACF. For a given delay $\tau$, we averaged the correlation $f(t_{m})f(t_{n})$ between all point pairs that are spaced apart by a time difference $t_{m}-t_{n}$ within the range $\tau\pm\Delta\tau/2$, where $\Delta\tau=30$ days is the bin width and the normalization constant is given by the number of point pairs that satisfy this condition, $N_{\tau}$. More precisely:

As a byproduct of the ACFs, we also found the characteristic time $\tau_{0}$ over which the values of $\Delta r_{k,i}-\overline{\Delta r_{k}}$ are correlated. For this purpose, we assume that the ACF follows an exponential function given by $R(\tau)=be^{-\tau/\tau_{0}}$, where $b$ and $\tau_{0}$ are free parameters that we fit to the ACF. In doing so, we find a characteristic time of $\tau_{0}=22.4\pm 6.6$ days, which corresponds to roughly the observing cadence. The fitted function and the derived parameters are presented in Fig. 6 panel (d).

For completeness, we repeated this analysis for a sample of other pulsars observed by NANOGrav using GBT. This includes high-DM pulsars such as PSR J1600$-$3053 ($52.33~{}\mathrm{pc}~{}\mathrm{cm}^{-3}$, Jacoby et al. 2007) and PSR J1744$-$1134, which is one of the lowest-DM pulsars observed by NANOGrav ($3.14~{}\mathrm{pc}~{}\mathrm{cm}^{-3}$, Demorest et al. 2013).

The broadband-narrowband offsets in $r_{\infty}$, $r_{\mathrm{2}}$, and $r_{\mathrm{\alpha}}$ for PSR J1744$-$1134 are presented in Fig. 5(b). We find that this pulsar yields smaller mis-estimations than those obtained for PSR J1643$-$1224. In particular, the mis-estimation in the infinite-frequency time of arrival, $\sigma_{\Delta r_{\infty}}$, is $4.4~{}\mu$s, which is $\sim 4.5$ times smaller than the value obtained for PSR J1643$-$1224. On the other hand, for J1600$-$3053 we find $\sigma_{\Delta r_{\infty}}=14.12~{}\mu$s, which is only $~{}\sim 1.3$ times smaller than $\sigma_{\Delta r_{\infty}}$ for PSR J1643$-$1224..

In order to study whether these differences are a result of J1744$-$1134 being less affected by interstellar dispersion, we also extend this analysis to other pulsars covering a broad range of DM values. The surveyed pulsars and the resulting values are summarized in Fig. 7. For each pulsar, we present the mis-estimation $\sigma_{\Delta r_{k}}$ as a function of the pulsar $\mathrm{DM}$ (top panel) and as a function of the RMS of its timing residuals (bottom panel), given that the former is another observational property of the pulsar that could have an effect on the mis-estimations. The results show no obvious dependence of the mis-estimation $\sigma_{\Delta r_{k}}$ with the pulsar $\mathrm{DM}$ or its residual RMS. If we set an arbitrary cut-off value close to the middle of our DM range, at $\mathrm{DM}=30$ pc cm^-3, we find that the average mis-estimation for lower-DM pulsars is $6.92~{}\mu$s, and for higher-DM pulsars it is $12.01~{}\mu$s. Broadly speaking, it can then be expected that high-DM pulsars will be generally more affected by dispersion mis-estimations in narrowband observations. However, the exact behavior is dependent on the specifics of the ISM in the pulsar line-of-sight, and no precise dependence of $\sigma_{\Delta r_{k}}$ as a function of the pulsar $\mathrm{DM}$ or its timing RMS can be established at this point.