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

As a radio signal is emitted at a pulsar and propagates through the ISM, it will encounter various frequency-dependent delays. As a result, pulses with frequency ν will arrive at Earth's position at a time t_ν that is delayed with respect to the expected time t_∞ for a signal of infinite frequency (i.e., assuming no chromatic effects). In this section, we describe the various frequency-dependent effects on the TOAs that cause 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 when no interstellar medium (ISM) intervenes. When compared against a single-frequency template, this variation introduces small systematic frequency-dependent perturbations in the TOAs. Following the modeling presented in Alam et al. (2021a), these perturbations are modeled as polynomials in log-frequency that are described by the FD_k parameters as

$$\begin{eqnarray}&&{t}_{\mathrm{PE}}=\displaystyle \sum _{k=1}^{n}{\mathrm{FD}}_{k}\mathrm{log}{\left(\displaystyle \frac{\nu }{1\,\mathrm{GHz}}\right)}^{k},\end{eqnarray} \tag{ 1 }$$

where t_PE is the profile evolution timing perturbation in units of seconds, ν is the observing frequency, and FD_k are the FD model parameters (Zhu et al. 2015). The number of terms needed varies for any given pulsar, but n = 2 parameters suffice to describe the pulse evolution with frequency for PSR J1643−1223. There is no k = 0 term because this would be a constant phase offset that is removed when the mean is subtracted from the timing residuals.

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 ν is delayed compared to one at infinite frequency by an amount t_DM = K × DM/ν², where the dispersion measure (DM) is the line-of-sight (LOS) integral of the free electron along the LOS 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 LOS l, and d is the pulsar distance. The dispersion constant is given by K = e²/2π m_e c ≃ 4.149 ms GHz²pc⁻¹ cm³.

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 LOS to sample electron-density fluctuations on a variety of scales (Phillips & Wolszczan 1991; Cordes & Rickett 1998; Lam et al. 2016). The result is a DM that varies with time, changing on timescales of hours to years. To model this time-varying dispersion, we use a stepwise model in which the 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 DM value is given by the epoch-dependent 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 DMX_i is the correction corresponding to the observing epoch i.

In addition to the ν⁻² offsets due to dispersion, ³⁴ there are also a variety of frequency-dependent effects for which the perturbations scale with radio frequency obeying a ν^−α 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 (α = 2; Foster & Cordes 1990). It also comprises pulse scattering, associated with the variable path length due to refraction, which means that the signal reaches the observer along different geometrical paths (e.g., Lewandowski et al. 2013; Bansal et al. 2019). 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 (α = 4; Lang 1971) or a Kolmogorov (α = 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 (Cordes et al. 1990). A simple model for the single-epoch TOA delay introduced by these chromatic effects is

$$\begin{eqnarray}&&{t}_{{\rm{C}},\nu }=\displaystyle \sum _{i=1}^{{N}_{c}}{C}_{i}{\nu }^{-{\alpha }_{i}},\end{eqnarray} \tag{ 2 }$$

where the i = 1, N_c additional terms model chromatic TOA variations with unique power-law spectral scalings α_i and power spectrum amplitudes C_i .

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

$$\begin{eqnarray}&&{t}_{\nu ,i}={t}_{\infty ,i}+\displaystyle \frac{K\times (\mathrm{DM}+{\mathrm{DMX}}_{i})}{{\nu }^{2}}+{t}_{\mathrm{PE},\nu }+{t}_{{\rm{C}},\nu ,i},\end{eqnarray} \tag{ 3 }$$

where the subscript i denotes the epoch.

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 LOS integral of the electron density. However, the LOS can change as a function of frequency because ray paths at different frequencies cover different volumes through the ISM (Cordes et al. 2016). Therefore, DM itself is frequency dependent, i.e., t_DM = K × DM(ν)/ν², 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 misestimation in the DM that is introduced by other chromatic effects in the observation that alter the DM fit. Other non-power-law frequency-dependent timing effects may be present at small amplitudes along specific LOSs 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 delays like this (see, e.g., Equation (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 Figure 1 we present a qualitative representation of the timing residuals (differences between the observed TOAs 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 the GASP bandwidth, the GUPPI bandwidth, or the full range of frequencies from 0.7 GHz to 1.9 GHz. For each backend, an artificial set of frequency-dispersed TOAs was created by simplifying Equation (3) as t_obs(ν) = K × DM/ν² + t_C/ν^4.4, using frequencies ν in the backend bandwidth. We choose K × DM = 1 GHz for the dispersion coefficient and t_C = 0.01 GHz for the chromatic coefficient. The "observed" TOAs were used to fit an "incomplete" timing model t_pred(ν) = a + b/ν² that only accounts for ν⁻² delays in order to simulate the effects of having other chromatic effects that are absorbed into the DM fit. The fitted functions are then evaluated at the same frequencies to obtain "predicted" TOAs. By subtracting both sets of TOAs, we obtain the residuals Δt = t_obs − t_pred presented in the figure. If the timing model were complete, we would obtain zero residuals in all cases. Instead, 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 the dispersion is modeled. The GUPPI backend combines nonsimultaneous observations around two frequency bands, one band centered near 820 MHz, and the other band 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.

Zoom In Zoom Out Reset image size

We expect that the effects of DM misestimations 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 × 10⁻²⁰) pulsar in a 147-day binary orbit with a white dwarf companion. This pulsar is of particular interest because it 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 (DM = 62.3 pc cm⁻³). Furthermore, PSR J1643−1224 has been shown to have significant scattering and profile shape variations (Shannon et al. 2016; 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 (with a spectral index γ_red = −1.3 and a spectral amplitude at f = 1 yr⁻¹, A_red = 1.619 μs yr^1/2), which may include contributions from unmodeled ISM propagation effects (Arzoumanian et al. 2018). Moreover, this pulsar has repeatedly been reported to exhibit an unusual chromatic timing behavior, including significant band-dependent noise (Lentati et al. 2016), pronounced structures in the autocovariance plots of its timing residuals (Perrodin et al. 2013), and suspected unaccounted-for chromatic timing biases (Arzoumanian et al. 2015). A dip in the timing residuals between 2015 February 21 and March 7 was found by Shannon et al. (2016), which is associated with a sudden change in the pulse profile. The authors report that when it is not modeled, this dip affects the upper limits on the stochastic GW background (see also Goncharov et al. 2020), so it has been included in subsequent timing models.

In Figure 2 we present the timing residuals for PSR J1643−1224 using the NANOGrav 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 misestimates in the GASP narrower bandwidth. In addition, we notice that epoch-averaged residuals in the GUPPI data set, especially those in the 800 MHz band, generally exhibit larger uncertainties and a broader spread than those taken with GASP in the same frequency band. In particular, the variance in the individual residuals in the 800 MHz band is ∼13.74 with GUPPI but ∼17.82 with GASP with similar error medians (∼2.5), so the larger error bars in the epoch-averaged residuals must result from the broader spread in the individual TOAs. These features are in agreement with the behavior predicted by Figure 1 for a timing model that is affected by DM misestimates 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 misestimations, and therefore, to exhibit higher residuals and uncertainties (orange curve in Figure 1), but these same frequencies might not be present (green curve in Figure 1). These results demonstrate the need for an in-depth analysis of the effects of interstellar dispersion and the misestimates in its modeling on timing residuals.

Zoom In Zoom Out Reset image size

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 ∼7.5 yr. We have separated these observations into two data sets, the original and a modified version:

• 1.  
  The full set of GUPPI broadband (BB) 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 (NB) approximation.
• 2.  
  In addition, we created an artificial set of NB observations by using Astropy (Astropy Collaboration et al. 2022) to filter out all the GUPPI BB observations outside the frequency ranges of the two GASP 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 the GASP bandwidth during the same time period.

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

Zoom In Zoom Out Reset image size

Using the time windows that are specified by the DMX parameters (see Section 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 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 BB observations reduces to 9604 and the set of NB 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 NB observations, first we created a simplified timing model that includes corrections for long-term interstellar dispersion and frequency-dependent profile evolution (the K × DM/ν² and t_PE,ν terms in Equation (3), respectively) but ignores the short-scale variations in the DM that are normally corrected for by the DMX parameters, as well as additional chromatic variations (the K × DMX/ν² and t_C,ν 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 (t_pred) at each of the observing dates and frequencies, which were then subtracted from observed TOAs (t_obs) calculated from each observation using a template-matching procedure. The difference will produce the timing residuals r = t_obs − t_pred.

As a result of ignoring the short-scale DM variations and additional chromatic delays, the residuals within each of the time windows will exhibit a curve close to ν⁻², as presented in Figure 4. We can then model the effects of the chromatic delays assuming the residuals can be described by

$$\begin{eqnarray}&&{r}_{i}(\nu )={r}_{\infty ,i}+{r}_{2,i}{\nu }^{-2}+{r}_{\alpha ,i}{\nu }^{-\alpha },\end{eqnarray} \tag{ 4 }$$

where r_2,i = K × DMX_i is the correction for temporal DM variations, and r_α,i is the correction due to a scattering-based delayed, corresponding to α = 4.4. Analogously to Equation (3), r_∞,i was devised to quantify the cumulative effects of the other achromatic timing parameters; the physical interpretation of r_∞,i is the residual expected at epoch i if no chromatic effects were present.

Zoom In Zoom Out Reset image size

We used LMFIT (Newville et al. 2014) to fit Equation (4) to the residuals within each time set, thereby obtaining best-fit values and uncertainties for r_∞, r₂, and r_α independently per epoch. Any unaccounted biases in the observations, such as DM misestimates, will result in biased fitted parameters. Note that most timing softwares fit all epochs jointly with the rest of the timing model and thus appropriately include long-term timing correlations. This would reduce the uncertainties in our measurements at the expense of significantly increasing the complexity of this analysis; for the exploratory purposes of this work, we therefore consider that this improvement is small and ignore correlations in the fits between epochs. Ultimately, if the frequency bandwidth did not introduce biases in the parameter estimation, we would find no offsets between the two data sets, regardless of any correlations.

This process was repeated using the NB and BB data sets separately (see Section 3). As a result, for each time window at epoch i, we obtained two sets of $\{{r}_{\infty ,i},{r}_{2,i},{r}_{\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 BB data set, and another from using the NB data set. For each r_k,i , we computed the parameter residuals

$$\begin{eqnarray}&&{\rm{\Delta }}{r}_{k,i}={r}_{k,i}^{\mathrm{BB}}-{r}_{k,i}^{\mathrm{NB}},\quad k=\infty ,2,\alpha \end{eqnarray} \tag{ 5 }$$

between the values fitted using the BB and the NB observations in each epoch. The error associated with each Δr_k,i is given by

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{\Delta }}{r}_{k,i}}=\sqrt{{\left({\varepsilon }_{{r}_{k,i}}^{\mathrm{BB}}\right)}^{2}+{\left({\varepsilon }_{{r}_{k,i}}^{\mathrm{NB}}\right)}^{2}},\end{eqnarray} \tag{ 6 }$$

where ${\varepsilon }_{{r}_{k,i}}^{\mathrm{BB}}$ is the fitting error from fitting r_k,i using the BB data, and ${\varepsilon }_{{r}_{k,i}}^{\mathrm{NB}}$ is the fitting error from using the NB data. Moreover, since we are only interested in the relative differences between the two sets of observations, we have subtracted the mean value of all the differences, $\overline{{\rm{\Delta }}{r}_{k}}$, from each of the differences Δr_k,i .

In the left panels of Figure 5, we plot the resulting values of ${\rm{\Delta }}{r}_{k,i}-\overline{{\rm{\Delta }}{r}_{k}}$ with their errors as a function of the MJD in the middle of the window. If the frequency bandwidth played no role in modeling interstellar dispersion, we would expect to obtain Δr_k,i = 0 for all the parameters k and all epochs i. However, the fact that we find nonzero deviations between the BB and NB sets of fitted parameters is indicative that using a narrower frequency bandwidth leads to misestimates in modeling these parameters for PSR J1643−1224. In order to quantify these misestimates, 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_∞, is ${\sigma }_{{\rm{\Delta }}{r}_{\infty }}=22.2\,\mu s$.

Zoom In Zoom Out Reset image size

For each fitted parameter, we plot in the right panels of Figure 5 histograms of the residuals ${\rm{\Delta }}{r}_{k,i}-\overline{{\rm{\Delta }}{r}_{k}}$ divided by their corresponding fitting errors ${\varepsilon }_{{\rm{\Delta }}{r}_{k,i}}$. We observe that for r₂ and r_α , the residuals are consistent with zero, which is indicative of a small offset introduced in the fit of these parameters when using NB data. However, the histogram for r_∞ reveals a skewed distribution (significantly more noticeable for J1643 − 1224), which suggests a systematic offset in the estimations of the infinite-frequency arrival times.

Next, in order to study the biases in the behavior of the ISM over time due to using NB receivers, we analyze whether the chromatic delay exhibits temporal correlations. In this case, we expect the pattern to be revealed in the autocovariance function (ACF) of this quantity. Therefore, we compute the ACF of ${\rm{\Delta }}{r}_{k,i}-\overline{{\rm{\Delta }}{r}_{k}}$ (see Equation (5)) for each parameter r_k between consecutive time windows. This process involves correlating a signal f(t) with a delayed copy f(t + τ) of itself as a function of the delay τ. However, since we have a discrete and irregularly sampled signal, we calculated a binned ACF. For a given delay τ, we averaged the products f(t_m )f(t_n ) between all point pairs of the signal f that are separated by a time difference t_m − t_n within the range τ ± Δτ/2, where Δτ = 30 days is the bin width, and the normalization constant N_τ is given by the number of point pairs that satisfy this condition. If we let R(τ) be the ACF, this function is then given by

$$\begin{eqnarray}&&R(\tau )=\displaystyle \frac{1}{{N}_{\tau }}\displaystyle \mathop{\displaystyle \sum _{{t}_{m}}\displaystyle \sum _{{t}_{n}}}\limits_{\mathop{\unicode{x0FE38}}\limits_{| {t}_{m}-{t}_{n}| \in [\tau \pm {\rm{\Delta }}\tau /2]}}f({t}_{m})f({t}_{n}).\end{eqnarray} \tag{ 7 }$$

The resulting ACFs for each of the fitted parameters are shown in Figure 6 in panels (a), (b), and (c).

Zoom In Zoom Out Reset image size

The information for the parameter misestimates, which we quantified in Section 4.2 as the standard deviation of the time series, is contained in the ACF as $\sqrt{R(\tau =0)}$. Moreover, we can also find the characteristic time τ₀ over which the values of ${\rm{\Delta }}{r}_{k,i}-\overline{{\rm{\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 τ₀ are free parameters that we fit to the ACF. The fitted function and the derived parameters are presented in panel (d) in Figure 6. We find a characteristic time of τ₀ = 22.4 ± 6.6 days, which approximately corresponds to the observing cadence. A better estimate of this timescale could be obtained by calculating the ACF using a finer bin width Δτ, but this would introduce noise in the resulting ACF. In order to gain better insight into the underlying correlation structure, higher observation cadences are needed, such as the cadence that is currently being employed by the CHIME Telescope (CHIME/Pulsar Collaboration et al. 2021).

The derived characteristic timescale implies that the misestimates in the timing parameters due to incomplete frequency sampling cannot be treated as independent in time, but rather they have a weak but non-negligible correlation among different epochs over a timescale as large as 3 weeks. Observing campaigns with a cadence longer than τ₀ will thus mitigate the correlation in the misestimates due to NB sampling, thereby reducing the sources of red noise in the observations at the expense of the number of observing samples. However, one can expect that the misestimations will still be present in the form of uncorrelated noise.

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 pc cm⁻³; Jacoby et al. 2007) and PSR J1744−1134, which is one of the lowest-DM pulsars observed by NANOGrav (3.14 pc cm⁻³; Demorest et al. 2013).

The BB-NB offsets in r_∞, r₂, and r_α for PSR J1744−1134 are presented in Figure 5. We find that this pulsar yields smaller misestimates than those obtained for PSR J1643−1224. In particular, the misestimates in the infinite-frequency time of arrival, ${\sigma }_{{\rm{\Delta }}{r}_{\infty }}$, is 4.4 μs, which is ∼4.5 times smaller than the value obtained for PSR J1643−1224. On the other hand, for J1600−3053, we find ${\sigma }_{{\rm{\Delta }}{r}_{\infty }}=14.12\,\mu s$, which is only ∼ 1.3 times smaller than ${\sigma }_{{\rm{\Delta }}{r}_{\infty }}$ for PSR J1643−1224.

In order to study whether these differences arise because J1744−1134 is less strongly 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 Figure 7. For each pulsar, we present the misestimates ${\sigma }_{{\rm{\Delta }}{r}_{k}}$ as a function of the pulsar 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 misestimations. The results show no obvious dependence of the misestimates ${\sigma }_{{\rm{\Delta }}{r}_{k}}$ on the pulsar DM or its residual rms. When we set an arbitrary cutoff value close to the middle of our DM range, at DM = 30 pc cm⁻³, we find that the average misestimates in the r_∞ parameter for lower-DM pulsars is 6.92 μ s, and for higher-DM pulsars, it is 12.01 μs. Broadly speaking, it can then be expected that high-DM pulsars will generally be more strongly affected by dispersion misestimations in NB observations. However, the exact behavior is dependent on the specifics of the ISM in the pulsar LOS, and no precise dependence of ${\sigma }_{{\rm{\Delta }}{r}_{k}}$ as a function of the pulsar DM or its timing rms can be established at this point.

Zoom In Zoom Out Reset image size