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

Noise characterization for pulsar-timing applications accounts for interstellar dispersion by assuming a known frequency dependence of the delay it introduces in the times of arrival (TOAs). However, calculations of this delay suffer from misestimations due to other chromatic effects in the observations. The precision in modeling dispersion is dependent on the observed bandwidth. In this work, we calculate the offsets in infinite-frequency TOAs due to misestimations in the modeling of dispersion when using varying bandwidths at the Green Bank Telescope. We use a set of broadband observations of PSR J1643−1224, a pulsar with unusual chromatic timing behavior. We artificially restricted these observations to a narrowband frequency range, then used both the broad- and narrowband data sets to calculate residuals with a timing model that does not account for time variations in the dispersion. By fitting the resulting residuals to a dispersion model and comparing the fits, we quantify the error introduced in the timing parameters due to using a reduced frequency range. Moreover, by calculating the autocovariance function of the parameters, we obtained a characteristic timescale over which the dispersion misestimates are correlated. For PSR J1643−1224, which has one of the highest dispersion measures (DM) in the NANOGrav pulsar timing array, we find that the infinite-frequency TOAs suffer from a systematic offset of ∼22 μs due to incomplete frequency sampling, with correlations over about one month. For lower-DM pulsars, the offset is ∼7 μs. This error quantification can be used to provide more robust noise modeling in the NANOGrav data, thereby increasing the sensitivity and improving the parameter estimation in gravitational wave searches.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
In 2010, the National Radio Astronomy Observatory (NRAO) launched the Green Bank Ultimate Pulsar Processing Instrument (GUPPI; DuPlain et al. 2008), a digital signal processor designed for pulsar observations with the 100 m Robert C. Byrd Green Bank Telescope (GBT).Its large bandwidth coherent dedispersion observation modes constituted a significant improvement over previously available backends (see Table 1), such as the Green Bank Astronomical Signal Processor (GASP; Demorest 2007).
The advent of receivers with larger bandwidths allows for better estimates of the dispersion measure (DM).These estimates are obtained by precisely measuring the arrival times of pulsar emission as a function of radio frequency (e.g., Manchester & Taylor 1977;Stairs 2002;Lorimer & Kramer 2004).Therefore, sampling more frequencies (up to a point; see e.g., Lam et al. 2018) allows for a better modeling of interstellar dispersion and mitigates misestimates introduced by other frequency-dependent delays.Conversely, DM misestimates will be more pronounced in observations with narrowband receivers because a smaller frequency space is sampled to model the dispersion of the signal.By comparing simultaneous observations in both frequency regimes, we can quantify the errors introduced by using a narrower frequency band.
Misestimates in the DM can arise from using narrowband frequency sampling (Shannon & Cordes 2017), incorrect temporal correlations in the DM when estimating it through interpolation techniques (Lee et al. 2014), pulse-amplitude modulation from interstellar scintillation that will shift the reference frequency (Lam 2016), the combination of asynchronously observed multifrequency data (Lam et al. 2015), or frequency-dependent DMs due to interstellar scattering (Cordes et al. 2016).These misestimates will cause red noise in the timing residuals (Keith et al. 2013).If it is unaccounted for, this noise can severely hinder the resulting timing precision (Lam et al. 2018).This is of special relevance for all high-precision pulsar-timing experiments involving frequency-dependent effects, such as gravitational wave (GW) searches (Agazie et al. 2023a;Antoniadis et al. 2023;Reardon et al. 2023;Xu et al. 2023), calculating pulse-broadening functions employing CLEAN deconvolution algorithms (Young & Lam 2024), monitoring interstellar scattering delays (Turner et al. 2021), and studying jitter in millisecond pulsars (Lam et al. 2019).Therefore, modeling and accounting for DM misestimates when using narrowband receivers is essential for providing realistic timing errors.
Observations with unaccounted-for errors due to DM misestimations are not suitable for high-precision timing experiments, so they are usually not included in studies like this (e.g., Archibald et al. 2018;Antoniadis et al. 2023).However, this approach reduces the available time baseline of pulsar observations, therefore decreasing our sensitivity to long-period gravitational waves.In this work, we perform narrow-bandwidth DM estimations using the GUPPI data set and compare the offset with the broader-bandwidth values to estimate the errors in the timing parameters.In turn, these error estimates can be used to provide more reliable calculations of the covariance matrix that is used for weighting narrowband times of arrival (TOAs), rendering them suitable for higherprecision timing applications.
This type of analysis will be of special interest as a new generation of wideband receivers approaches.A wider observing bandwidth will enable us to calculate misestimations due to incomplete frequency sampling, which were previously unaccounted for because of bandwidth limitations.Calculating these differences in the timing error due to different observing bandwidths will be of foremost importance for integrating current and future sets of broadband observations.In particular, the introduction of the VErsatile GBT Astronomical Spectrometer (VEGAS; Bussa et al. 2012) can duplicate all the capabilities of the GUPPI backend, but also allows for wider instantaneous bandwidths of up to 3.8 GHz.
In Section 2.1 we describe the different frequency-dependent delays affecting the signal propagation, and how they are incorporated into our timing model.In Section 3 we provide information on the data collection and reduction methods used for the NANOGrav observations with the GBT.In Section 4 we describe the main pulsar we analyzed in this work as a case study, PSR J1643−1224.We also outline how we obtained a narrowband data set from a broadband data set to carry out an optimal comparison between the two.In Section 5 we summarize the results and implications of this work.The processed data products presented here are publicly available33 as of the date this work is published (Sosa Fiscella & Lam 2023).

Frequency-dependent Timing Delays
As a radio signal is emitted at a pulsar and propagates through the ISM, it will encounter various frequencydependent 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 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/ν 2 , 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 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 /2πm e c ; 4.149 ms GHz 2 pc −1 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 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 n , where DMX i is the correction corresponding to the observing epoch i.
In addition to the ν −2 offsets due to dispersion,34 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 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 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(ν)/ν 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 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 plasmalens-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.

Simulation of Chromatic Delays
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/ν 2 + 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/ν 2 that only accounts for ν −2 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.

Pulsar Backends and Observations
In the present analysis, we used observations from NANOGrav's 12.5-year data set release (Alam et al. 2021a) obtained using the GBT at the Green Bank Observatory in West Virginia, USA.Two radio receivers at separated frequency bands were used to perform the observations: one band centered near 820 MHz, and the other band near 1400 MHz (see Table 1).The observations were performed with a monthly cadence using both receivers.However, these two separate frequency ranges were not observed simultaneously.Instead, the observations were separated by a few days due to the need for a physical receiver change at that telescope.The typical observation duration was about 25 minutes.
Two generations of pulsar backend processors were used for real-time coherent dedispersion and folding of the signal: 1. GASP was used from the start of the NANOGrav observing program in 2004 until its decommissioning in 2012.It decomposed the signal into contiguous 4 MHz channels over a bandwidth of 64 MHz (Ferdman 2008).2. Starting in 2010, GASP was replaced by GUPPI, a wideband system that can process up to 800 MHz in bandwidth using smaller 1.5625 MHz channels, and this significantly improved the timing precision relative to GASP (Ford et al. 2010).During the transition from GASP to GUPPI, precise measurements of time offsets between the instruments were made and included in the residual calculation (Alam et al. 2021a).
The observations were calibrated and analyzed using standard pulsar processing techniques as implemented in the code PSRCHIVE (Hotan et al. 2004) within the NANOGrav data reduction pipeline (Demorest 2018).In brief, the backend divides the telescope passband into narrow spectral channels, undertakes coherent dedispersion of the signals within each channel, and folds the resulting time series in real time using a pulsar-timing model.The data were thus transformed into folded pulse profiles as a function of time, pulsar phase, radio frequency, and polarization.These profiles have 2048 phase bins across the pulsar spin period, a frequency resolution of 4 MHz (GASP) or 1.5 MHz (GUPPI), and a time resolution (subintegration time) of 1 and 10 s, respectively (Arzoumanian et al. 2018).
Care was taken to remove all artifacts that will result in a frequency-dependent systematic TOA bias.This includes removing image rejection artifacts that could arise from running two interleaved analog-to-digital conversion schemes if the gain of the two converters is not identical.Furthermore, the data set cleaning pipeline also involved systematically removing radio-frequency interference, excluding low-signalto-noise ratio TOAs , where we observe a characteristic dispersion curve as a result of dispersion (dashed gray line), or dispersion and scattering (solid black line).Bottom panel: Expected timing residuals Δt = t obs − t pred when the "incomplete" model t pred (ν) = a + b/ν 2 is fitted to the dispersion-and scatteringaffected TOAs (solid black line in the upper panel) using the GASP bandwidth (green), the GUPPI bandwidth (orange), and the full range of frequencies (blue).
observations affected by calibration or digitization errors, and manual inspection of the data sets.The full details regarding data collection, calibration, pulse arrival-time determination, and noise modeling for the NANOGrav 12.5-year data set are provided in Alam et al. (2021a).

PSR J1643−1224
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 −20 ) 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 −3 ).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 −1 , 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 banddependent 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 epochaveraged 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.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.49MHz.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.
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.

Biases in the Timing Parameters
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/ν 2 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  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 ν −2 , as presented in Figure 4. We can then model the effects of the chromatic delays assuming the residuals can be described 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.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 2 , 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

{ } { } values with their corresponding fitting errors e 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 ) , is the fitting error from fitting r k,i using the BB data, and e r NB k i , 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, Dr k , from each of the differences Δr k,i .
In the left panels of Figure 5, we plot the resulting values of 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 s m =  The frequency ranges corresponding to the two GASP receivers, Revr_800 and Revr1_2, are shaded in green and red, respectively.Because the simplified timing model does not include DMX and additional chromatic corrections, the residuals are not distributed around R = 0, but follow a dispersion curve that can be fitted using Equation (4).The fitted values for this window and their errors are presented in the text boxes.
their corresponding fitting errors e Dr k i , .We observe that for r 2 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.

Autocovariance Functions
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 D -D r r k i 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 + τ) 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 NB between the values that were fitted at epoch i using the BB and the NB data sets.We also report the mean value of all the differences, Dr k , and the standard deviation of each set of residuals, s Drk .Since only differences are relevant, the mean value of the residuals is subtracted in each case.For each parameter, we present in the right panel histograms of the number of residuals divided by their corresponding fitting errors e Drk i , (see Equation ( 6)).
function is then given by The resulting ACFs for each of the fitted parameters are shown in Figure 6 in panels (a), (b), and (c).
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 t = R 0 ( ).Moreover, we can also find the characteristic time τ 0 over which the values of D -D r r k i k , are correlated.For this purpose, we assume that the ACF follows an exponential function given by t = t t -R be 0 ( ) , where b and τ 0 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 τ 0 = 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 τ 0 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.

Other Pulsars
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 −3 ; Jacoby et al. 2007) and PSR J1744−1134, which is one of the lowest-DM pulsars observed by NANOGrav (3.14 pc cm −3 ; Demorest et al. 2013).
The BB-NB offsets in r ∞ , r 2 , 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, s D ¥ r , 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 s m = D ¥ s 14.12 r , which is only ∼ 1.3 times smaller than s D ¥ r 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 s Dr 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 s Dr 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 −3 , 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 s Dr k as a function of the pulsar DM or its timing rms can be established at this point.

Conclusions
In this work, we have shown that when the pulsar frequency spectrum is sampled using NB radio receivers, biases are introduced in the timing parameters describing chromatic delays.These biases arise from leftover chromatic components that are absorbed into r ∞ , a parameter that quantifies the cumulative achromatic effect of the other timing parameters.These misestimates could contribute to the structures in the residuals and the variations in TOA uncertainty observed in Figure 2, thus limiting the precision of pulsar-timing campaigns using NB radio receivers.
This effect is dependent on the DM of a given pulsar, and for a high-DM pulsar such as PSR J1643−1224, we find an offset as large as 22.2 μs.Since timing models depend on these parameters to calculate the timing residuals, this analysis provides an estimate of the error that is consequently introduced in the residuals as a result of DM misestimations.For J1643−1224, we find that the errors in observations at different epochs are correlated within about one month of each other.
For low-DM pulsars such as J1744−1134, the bias in the timing parameters reduces to 8.1 μs.However, in Figure 7 we do not find a clear linear dependence of the error as a function of the DM; thus we infer that the exact offset is highly dependent on the pulsar properties.Nevertheless, for a typical DM value (40 pc cm −3 ), the systematic offset will be ∼5 μs.
The resulting estimates of the TOA uncertainties due to NB frequency sampling and their temporal correlations (both of which are contained in the ACF calculated in Section 4.3) can be used to provide more reliable calculations of the covariance matrix for the TOA data.This matrix is used to account for both uncorrelated and correlated noise across the NB TOAs obtained for a single observing epoch by weighting the TOAs accordingly (Agazie et al. 2023c).Better estimations of the TOA uncertainties will result in NB and BB TOAs being weighted differently, as shown in our work; therefore, providing these uncertainties in the covariance matrix will allow for a more robust accounting for the errors introduced by TOAs with NB biases.In general, there remains a mismatch in the phenomenological noise model used by most PTA collaborations and the physically motivated noise models they are working to develop (Agazie et al. 2023c).The source of uncertainty considered in this work manifests in an unknown way as a component in the phenomenological noise model, but it can be exactly added in the physical noise-model case and with the appropriate time correlation, as discussed in Section 4.3.
The most immediate application of these results will be quantifying the previously unaccounted-for error in the NANOGrav legacy observations in order to incorporate them into the current data set.Legacy data comprise years of observations that are already available, and which could significantly strengthen current evidence for a detection of the long-period GW background (Agazie et al. 2023a).For PTAs such as are observed by NANOGrav, where the GW background power dominates the intrinsic pulsar noise and is therefore in the strong-signal regime, the signal-to-noise ratio will scale as the square root of the time baseline (Siemens et al. 2013;Pol et al. 2022).However, for PTAs in the weak-signal regime, where the lowest frequency of the stochastic background power spectrum is below the white-noise level, the signal-to-noise ratio scales with time as roughly the fourth power of the time baseline.This is the case of the Indian Pulsar Timing Array (Antoniadis et al. 2023) and the Pulsar Monitoring in Argentina collaboration (Gancio et al. 2020).Therefore, incorporating even a few years of legacy observations will result in a significant increase in their sensitivity.
The new generation of BB radio receivers will substantially improve our estimations of the chromatic parameters resulting from the improved lever-arm effect of a larger bandwidth (see the limitations discussed in Lam et al. 2018).The subsequent effects on achromatic timing parameters discussed in this paper will be minimized with systems like this.In particular, we find that DM misestimates are significantly mitigated by using BB radio receivers, such as the ultrawideband receiver (Bulatek & White 2020) at the GBT.Moreover, the upcoming Deep Synoptic Array-2000(DSA-2000) telescope, which will produce pulsar observations across the entire 0.7-2 GHz band (Hallinan et al. 2019), is expected to provide major improvements in current ISM models and to significantly reduce the residual errors due to biased DM parameters.Repeating the analysis presented in this work using DSA-2000 BB observations could potentially contribute to better constraining DM misestimates.
Even as we move toward ultrawideband systems, this work has the potential to improve already existing NB data sets.In particular, we expect DM-induced biases to still be prevalent in GASP-and GUPPI-based observations (see Figure 2).NB observations are also predominant, for example, in the European Pulsar Timing Array data release 1.0 (Desvignes et al. 2016), which subsequently played a major role in the International Pulsar Timing Array (IPTA) data release 2 (Antoniadis et al. 2022).In particular, the Effelsberg Radio Telescope processes observations with a bandwidth up to 112 MHz (Backer et al. 1997) and the Lovell Radio Telescope up to 128 MHz, the Nançay Radio Telescope uses a coherent dedispersion backend of the same family as GBT's ASP-GASP with a bandwidth of either 64 or 128 MHz, and the Westerbork Synthesis Radio Telescope uses bandwidths of 10, 80, or 160 MHz.Therefore, quantifying the DM misestimates introduced by these NB systems is of the utmost relevance for the search of an isotropic stochastic GW background in the IPTA data set.
These results are also of interest for other radio-astronomical facilities that are currently using NB receivers.For example, the Argentine Institute of Radio Astronomy has two single-dish telescopes capable of performing daily pulsar monitoring, which could contribute to improving the IPTA sensitivity to single sources of GWs (Lam & Hazboun 2021).However, their observations use instantaneous bandwidths of 112 MHz and 56 MHz (Zubieta et al. 2023).Therefore, using this type of analysis to account for DM misestimates might prove of foremost importance in achieving the timing precision required for them to be contributing to future IPTA data sets.
(see details in NANOGrav Collaboration et al. 2015), removing outliers identified by a Bayesian analysis of residuals (see details in Arzoumanian et al. 2018), removing

Figure 1 .
Figure 1.Upper panel: Artificial set of TOAs created as t obs(ν) = K × DM/ ν 2 + t C /ν 4.4 , where we observe a characteristic dispersion curve as a result of dispersion (dashed gray line), or dispersion and scattering (solid black line).Bottom panel: Expected timing residuals Δt = t obs − t pred when the "incomplete" model t pred (ν) = a + b/ν 2 is fitted to the dispersion-and scatteringaffected TOAs (solid black line in the upper panel) using the GASP bandwidth (green), the GUPPI bandwidth (orange), and the full range of frequencies (blue).

Figure 2 .
Figure 2. Timing residuals for PSR J1643−1224 using the NANOGrav 12.5 yr data set (Alam et al. 2021a).The predominant data acquisition backend instrument over any given time period is indicated at the top of each figure, and the vertical dashed lines indicate the times at which instruments changed.The colored points indicate the receiver: Rcvr_800 MHz (blue for GASP, green for GUPPI) and Rcvr1_2 (orange for GASP, pink for GUPPI).Top panel: Residual arrival times for all TOAs.The points are semitransparent, and opaque regions arise from the overlap of many points.Middle panel: Average residual arrival times shown in full scale.Each observation is composed of many simultaneously obtained NB TOAs at different frequencies.Bottom panel: Close-up of residuals around zero.

Figure 3 .
Figure 3. Observing frequency as a function of the MJD for each observation in our data set.The typical frequency bands corresponding to each of the GASP receivers, Revr_800 and Revr1_2, are highlighted in purple and red.The full set of observations taken with GUPPI, covering a wide frequency range, is presented in purple.The subset of the observations that were used to emulate a NB data set is presented in red.
fitted using the BB and the NB observations in each epoch.The error associated with each Δr k,i is given by

.
For each fitted parameter, we plot in the right panels of Figure5histograms of the residuals D -

Figure 4 .
Figure4.Example of fitting the r ∞ , r 2 , and r α parameters using the residuals that are obtained by subtracting the predicted TOAs from the simplified timing model and either the BB (panel (a)) or NB (panel (b)) set of observed TOAs of J1643−1224 that fall within the window from MJD = 55305.17896to 55307.17351.The frequency ranges corresponding to the two GASP receivers, Revr_800 and Revr1_2, are shaded in green and red, respectively.Because the simplified timing model does not include DMX and additional chromatic corrections, the residuals are not distributed around R = 0, but follow a dispersion curve that can be fitted using Equation (4).The fitted values for this window and their errors are presented in the text boxes.

Figure 5 .
Figure 5. Fitted parameters for PSR J1643−1224 in panel (a) and PSR J1744−1134 in panel (b).For each fitted parameter r k in Equation (3), each panel represents at the left the parameter residual D = r r r

Figure 6 .
Figure 6.Autocovariance functions for PSR J1643−1224.Panels (a), (b), and (c) show the ACFs for the values of D -D ¥ ¥ r r , D -D r r 2

Figure 7 .
Figure 7. Top panel: Standard deviation of Δr ∞ , Δr 2 , and Δr α for the different pulsars we have included in this study as a function of the pulsar DM.Each pulsar is represented with a different marker, and each parameter is represented with a different color.The DM values were obtained from Demorest et al. (2013).Bottom panel: Standard deviation of Δr ∞ , Δr 2 , and Δr α as a function of the weighted rms of post-fit timing residuals.The rms values were obtained from Alam et al. (2021b).