{'fraction_non_alphanumeric': 0.03356529632010892, 'fraction_numerical': 0.023070231836919936, 'mean_word_length': 4.465943152454781, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Frequency Ranges for Sub-Microsecond Precision Pulsar Timing

Precision pulsar timing requires optimization against measurement errors and astrophysical variance from the neutron stars themselves and the interstellar medium. We investigate optimization of arrival time precision as a function of radio frequency and bandwidth. We find that increases in bandwidth that reduce the contribution from receiver noise are countered by the strong chromatic dependence of interstellar effects and intrinsic pulse-profile evolution. The resulting optimal frequency range is therefore telescope and pulsar dependent. We demonstrate the results for five pulsars included in current pulsar timing arrays and determine that they are not optimally observed at current center frequencies. For those objects, we find that better choices of total bandwidth as well as center frequency can improve the arrival-time precision. Wideband receivers centered at somewhat higher frequencies with respect to the currently adopted receivers can reduce required overall integration times and provide significant improvements in arrival time uncertainty by a factor of ~sqrt(2) in most cases, assuming a fixed integration time. We also discuss how timing programs can be extended to pulsars with larger dispersion measures through the use of higher-frequency observations.


INTRODUCTION
Pulsars have been used in some of the most constraining tests of gravity and general relativity (see Will 2014, for an overview), as laboratories for super-dense nuclear equation-of-state experiments (Demorest et al. 2010;Lattimer & Prakash 2016), and as detectors of low-frequency (nanohertz to microhertz) gravitational waves (Zhu et al. 2014;Wang et al. 2015;Lentati et al. 2015;Arzoumanian et al. 2015a;Shannon et al. 2015;Lasky et al. 2016;Arzoumanian et al. 2016;Babak et al. 2016). Observations of pulse times of arrival (TOAs) allow us to develop a timing model that tracks every rotational phase of the pulsar (Verbiest et al. 2009;Cordes & Shannon 2010). Advancing pulsar timing tests will require the highest arrival-time precision possible Lam et al. 2017).
Many frequency-dependent effects distort pulsar signals along the line of sight (Cordes 2013;Stinebring 2013). Intrinsic pulse emission varies as a function of frequency (e.g., Pennucci et al. 2014;Liu et al. 2014). The interstellar medium (ISM) causes both time and frequency-dependent pulse shape and flux changes as well as arrival time delays from numerous optical effects (Palliyaguru et al. 2015;Levin et al. 2016;Lam et al. 2016a). Measurement at the Earth is limited to specific receiver instrumentation and frequency ranges, as well as from observation through specific radio-frequency-interference (RFI) environments (e.g., Verbiest et al. 2016). The TOAs and their uncertainties, and therefore the ultimate achievable timing precision, will be affected by the choice of radio frequencies with which we observe.
This paper considers optimization of timing precision as a function of observing frequency and bandwidth. Our results are generally applicable to any high-timing-precision experiment with pulsars. We apply our results to several pulsars observed by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; McLaughlin 2013) collaboration. NANOGrav observes many pulsars with frequencies ranging from ∼0.3−2.5 GHz but our analysis is focused on a select few.
In §2, we describe the various frequency-dependent effects in pulsar timing and develop the methodology to compute a single TOA uncertainty. Table 1 provides a summary of the different effects we discuss and the general forms of the equations. We discuss implications for specific high-precision millisecond pulsars (MSPs) observed by NANOGrav (hereafter "NANOGrav pulsars" for short) in §3 and then expand our analysis to other MSPs broadly in §4. We conclude in §5.

FREQUENCY-DEPENDENT TOA UNCERTAINTY
In this section, we discuss the various components contributing to pulse TOA uncertainty as a function of frequency. Our goal is to compute the TOAs and uncertainties referenced to the infinite-frequency arrival time, t ∞ . For frequencies ν 1 and ν 2 with ν 1 < ν 2 , we define the bandwidth B = ν 2 − ν 1 and where convenient the frequency ratio r = ν 2 /ν 1 . The observing timespan will be given by T . Throughout the paper, we use the subscript '0' to denote a parameter referenced to a center frequency of ν 0 unless otherwise specified.
Our approach in this paper is to examine the frequency-dependent timing and noise model we use to describe our TOAs and then formulate an uncertainty on t ∞ . Similar to Cordes & Shannon (2010) and , we consider our timing model for frequencydependent TOAs t ν to be where KDM(ν)/ν 2 is the dispersive delay given by a frequency-dependent dispersion measure (DM) with dispersion constant K ≈ 4.149 ms GHz 2 pc −1 cm 3 (Lorimer & Kramer 2012;Cordes et al. 2016), t C,ν represents an additional systematic chromatic delay term, and ν is the TOA uncertainty measured at specific frequency ν. In general, ν will contain components from multiple white-noise sources described by a covariance matrix C νν , where ν and ν denote two separate observing subbands. The three white-noise (i.e., uncorrelated, for this case considered in time) components of our covariance matrix are additive, i.e., where S νν is from template fitting, J νν is from pulse jitter, and D νν is from diffractive interstellar scintillations (DISS). Template-fitting errors are calculated from a matched-filtering procedure using a template shape compared with a data profile containing additive noise (such as from radiometer noise) that varies in both time and frequency (see Lam et al. 2016a, for more details). Pulse jitter results from individual pulses varying stochastically, breaking the matched-filtering assumption of the averaged pulse profile being an exact copy of the template plus additive noise, and thus is uncorrelated between pulses in time though does have some correlation in frequency. For DISS, the stochastically varying ISM along the line of sight similarly results in observed pulses deviating from the template shape which are uncorrelated in time and frequency beyond the characteristic scintillation time and frequency scales (Levin et al. 2016).
All three effects will be discussed in more detail in later subsections ( §2.2.1, §2.2.2, §2.2.3, respectively). Multiple effects can cause TOA uncertainties that are correlated over longer timescales than a single epoch though we do not discuss them here in this paper; errors on any timing parameters will typically be small due to the long-term timing fit.
We will use the combined covariance matrix to determine the white-noise uncertainty averaged over the single observation. Additional error terms will be added in quadrature to this single white-noise uncertainty. We will describe each of the noise components in detail before discussing the DM fitting procedure. For brevity, we will drop the subscripts on all matrices.
Rather than deal with N number of independentlymeasured pulses across frequency channels over our bandpass and therefore N template-fitting uncertainties, we choose to combine the errors into a single number, thus providing a metric to compare TOA uncertainties as a function of center frequency and bandwidth. We compute the single white-noise uncertainty component on t ∞ following the procedure in Arzoumanian et al. (2015a), where the matrix U is the "exploder matrix" defined in Arzoumanian et al. (2014) that in general groups the multifrequency arrival times taken over many observations into specific epochs and for this work specifically groups the frequency-dependent arrival times into a single epoch by setting U = col(1), a column matrix where each element is 1. Again, the covariance matrix C is the sum of the three separate covariance matrices S, J, and D. The white-noise uncertainty σ W is then the square root of the single element of C E . We define the total TOA uncertainty (squared) on t ∞ as where the variances in order are the squares of: the white-noise errors, the DM estimation error, and the telescope error. We wish to reinforce that σ TOA is separate from ν , which is the per-frequency uncertainty on t ν . We will describe the three components in depth in the following subsections.

White-Noise Variations
The three components of the white-noise TOA uncertainties are given by where σ S/N is the template-fitting error, σ J is the jitter error, and σ DISS is the scintillation error. All three are discussed in detail in Lam et al. (2016a). While the three terms have been shown to be uncorrelated between epochs (and thus are white noise), σ J is correlated in frequency and σ DISS has a characteristic time and frequency correlation structure, discussed shortly.

Template-Fitting Error
We estimate TOAs calculated with a template-fitting approach (Taylor 1992). In order to calculate the template-fitting error, we assume that the pulse profile is an exact copy of a template shape with additive noise. The minimum TOA error can be calculated from the effective width of the pulse W eff and the S/N (peak to off-pulse rms, written as S for clarity in equations), where N φ is the number of samples (bins) across pulse phase φ from 0 to 1. Changes in N φ will change the rms noise and thus S; however the product SN 1/2 φ is invariant (Lam et al. 2016a). Use of the equation assumes that W eff is constant over the band. We note that in the low-S/N regime, the true template-fitting error is larger than given by Eq. 6 because the template-matching procedure begins to fail and the error becomes significantly non-Gaussian (Arzoumanian et al. 2015b, Appendix B).
Since the pulse profile changes over the frequency range, the template-matching assumption fails if we are to average the pulse over the band. Pulse profiles tend to broaden at low frequencies intrinsically (though not necessarily), pulse profile components can overlap in pulse phase or appear/disappear altogether, and multipath scattering broadens pulses. Therefore, to compute the template-matching error, we assume that the pulse shape in a sufficiently narrow frequency channel is constant.
The assumption has the added benefit of allowing for a single S for a narrow frequency channel.
The effective width of the pulse at a given frequency can be written as where P is the pulsar spin period, τ d,0 is the scattering timescale at reference frequency ν 0 , and the observed pulse template U obs scaled to a peak of unity is equal to the intrinsic pulse shape convolved with the pulse broadening function (PBF), The differences in the denominator in Eq. 7 are related to the derivative of the template and so pulses with sharp features will have a smaller W eff ; the effective width is not equivalent to the width of a pulse component. Since τ d will change as a function of frequency proportional to ν −22/5 (see §2.2.3), the width of the PBF will also change as a function of frequency. The PBF that conserves pulse fluence is approximately a negative exponential which causes the pulse not only to broaden but the amplitude to decrease (Bhat et al. 2003;Geyer et al. 2017). Therefore, the effective width increases and the pulse S/N decreases, causing drastic changes to σ S/N at lower frequencies. Now we consider the pulsar signal, the numerator of the signal-to-noise ratio. Above frequencies of 100 MHz, the period-averaged flux densities can be described by a simple power law, where I 0 is the intensity constant and α is the spectral index (Lorimer & Kramer 2012). The mean value of α = −1.60 ± 0.03 (Jankowski et al. 2017). For a pulsar signal with a power-law flux density in frequency, the period-averaged signal-to-noise ratio is where the numerator represents the pulsar signal while the denominator and factor in the square root represents the rms off-pulse noise from the radiometer equation (the factor of N pol = 2 represents the number of polarizations). In order to obtain the peak-to-off-pulse signalto-noise ratio S, we must multiply Eq. 10 by a factor of whereŪ obs (ν) is simply the mean of U obs (ν), calculated . The factor of 1.0 in the numerator comes from the template being scaled to peak unity as previously discussed. We write the ratio as S following the convention in Lam et al. (2016a) for brevity in equations. Again, B is the observing bandwidth and T is the observation duration. The factor of A e /2k = (A e /2760 m 2 ) K/Jy in the denominator is the gain that converts the system temperature (T sys ) equivalent into the total measured noise power in flux density units.
In addition to altering the effective width W eff of the pulse profiles, multipath scattering will cause pulse profiles to have lower S/N since the total flux will be broadened across pulse phase. Since the scattering timescale τ d varies as a function of frequency (proportional to ν −22/5 for a Kolmogorov medium), the pulsar signal term in the numerator will therefore take a more complicated, frequency-dependent form, I(ν|τ d,0 ). Free-free absorption has a negligible effect for pulsar-frequency combinations used for sub-microsecond timing precision.
The system temperature T sys is a combination of all temperature contributions measured by the antenna, such as from the Cosmic Microwave Background (CMB; T CMB ), the receiver bandpass (T rcvr ), the Galactic background (T Gal ), etc., and can be written as the sum of those contributions: The response of the receiver is also not constant over the band and will therefore also affect S(ν), however we will assume here that it is constant and group it together with the CMB temperature to form T const . For the frequency ranges we consider, we use the form in Cortés-Medellín (2004) for the Galactic background, At frequencies above about 10 GHz, a thermal Galactic component can become dominant over this nonthermal component. However, for the range of frequencies that we consider here, this thermal component will always be sub-dominant for the vast majority of lines of sight and we do not consider it here. The Galactic component will also have some Galactic latitude and longitude dependence that changes both the coefficient and the spectral index of T Gal (Cortés-Medellín 2004; Kogut et al. 2011). In addition, the Sun will also increase T sys as a function of time; both the intrinsic emission and the angular size on the sky will have frequency dependence (Oberoi et al. 2017). These effects primarily dominate at the lowest frequencies we will consider and cause the template-fitting uncertainty to increase. We will show later that the optimal frequency ranges tend to be at much higher frequencies and these spatial-and time-dependent contributions will be negligible there. Therefore, we will ignore them in our analysis for simplicity though they can be considered for future analyses depending on the applications.
Combining Eqs. 6, 7 and 10, the template-fitting error can be written as where we have explicitly included a scattering timescale that modifies both the pulse shape and intensity. The template-fitting covariance matrix is then diagonal with components

Pulse Phase Jitter
While average pulse profiles at a specific frequency are stable over time, the shapes of individual pulses stochastically vary (Cordes & Downs 1985). Known as pulse jitter, a finite sum of pulses will deviate from the template slightly, breaking the template-fitting assumption of the profile being the sum of the template plus additive noise. Therefore an additional TOA uncertainty must be taken into account. Since the average pulse shapes change as a function of frequency, we expect that pulse jitter will change as a function of frequency such that when summed they yield the observed frequencydependent pulse profile evolution; comparisons between the two regimes, however, have not been studied previously. Shannon et al. (2014) found that the jitter noise decorrelates over a bandwidth of approximately 2 GHz for observations of time-averaged pulses covering ∼0.7−3.6 GHz for PSR J0437−4715. Lam et al. (2016a) found statistically different values of jitter noise for many MSPs in the NANOGrav pulsar timing array (PTA) as a function of receiver band observed. The values are expected to be correlated between frequency bands though an analysis has not yet been undertaken.
For narrow bands over which the rms jitter σ J is constant, the covariance matrix J has all components set to σ 2 J . In general however the jitter covariance matrix is given by where ρ J describes the correlation of the variance between the two frequency bands and here we set ρ = 1 for simplicity, roughly consistent with Shannon et al. (2014).

Scintillation Noise
DISS modulates the observed pulse S/N by a gain g (different from the telescope gain) with a probability density function (PDF) given by where S 0 is the mean S/N, n ISS is the number of scintles, Γ is the gamma function, and Θ is the Heaviside step function (Lam et al. 2016a). When the number of scintles n ISS is large, the PDF tends towards a constant (the mean) S/N, f S (S|n ISS ) → δ(S − S 0 ). The number of scintles n ISS scales as where ∆t d and ∆ν d are the scintillation bandwidth and timescale respectively, B is again the total observing bandwidth, and T is the observing duration. The filling factors η t , η ν are typically in the range 0.1 to 0.3 depending on the scattering geometry along the line of sight (Cordes & Shannon 2010;Levin et al. 2016). For a Kolmogorov medium, the scintillation parameters scale with frequency (Cordes & Rickett 1998;Cordes & Lazio 2002) as In addition to modulating the pulse gain, DISS causes the measured PBF to vary stochastically. Like pulse jitter, the epoch-to-epoch deviation from the template shape introduces an additional error that must be accounted for in the TOA uncertainty (the finite scintle effect, see Cordes et al. 1990 for a more thorough discussion). When added in quadrature to the TOA uncertainty, the scintillation noise component takes the form when n ISS is large and τ d when there is only one scintle or a partial scintle across the band ). The pulse-broadening timescale τ d is related to the scintillation bandwidth by where C 1 is a coefficient of order unity. We assumed in our analysis that the ISM is a fully-uniform Kolmogorov medium and adopted C 1 = 1.16 and η t = η ν = 0.2 (Cordes & Rickett 1998;Cordes & Shannon 2010). The DISS covariance matrix is where like ρ J , ρ DISS describes the correlation of the variance between the two frequency bands though it depends on the scintillation time and frequency correlation scales at both frequencies. In the application of our analysis, we will make a simplifying assumption that ρ DISS is split into two regimes, where there are multiple scintles n ISS in the band and therefore uncorrelated and where n ISS ≈ 1 and so ρ DISS → 1. Our assumption is discussed further in §3 when we apply our analysis to specific MSPs.

DM Fluctuations
One of the largest timing errors from ISM effects comes from DM misestimation. Typically, pulsars are observed at multiple frequencies covering a wide band in order to gain infinite-frequency timing precision. We include three contributions to the DM estimation error: from additive white-noise uncertainties (σ DM ), from frequency-dependent DM (σ DM(ν) ), and from systematic changes in the PBF (σ δt C ). Stochastic changes from the finite-scintle effect are included in the covariance matrix whereas we let the chromatic delay term t C,ν represent unknown changes in the scattering delay timescale over many epochs.
In this work, we consider only a fit for the dispersive delay, which follows most PTA procedures up to now. We follow the formalism in Appendix D of Cordes & Shannon (2010) to calculate the DM estimation error over many frequency channels, taking into account that the covariance matrix we consider is not diagonal. The timing model can be written in matrix form as T = Xθ + , where T = col(t k ) is the column vector of arrival times for the k-th frequency, X is the design matrix given by matrix(1 ν 2 k ) from a fit for only t ∞ and dispersive delays, θ is the parameter vector col(t ∞ , KDM), and are the errors described by the covariance matrix C. In the case of the design matrix we use here, we are assuming a timing model as in Eq. 1 but with a fit only for a constant DM over the observation, i.e., we do not fit for the additional chromatic delay since again we are only considering the model of current PTA procedures and we do not yet consider the frequency-dependence of DM ). The standard error on t ∞ from a fit for DM, which we denote σ DM , is given by the first element of the covariance matrix for the parameter errors (X T C −1 X) −1 . In Table 1, we provide the simplified form in the case where only two individual frequencies are used.
We model systematic variations in the PBF as a powerlaw scaling for the chromatic delay, t C,ν = a C ν −X . Foster & Cordes (1990) assume X = 4 because they assume that the geometric increase in the propagation path is the greatest contribution after the dispersive delay (and the barycentric angle-of-arrival correction, see the discussion in Lam et al. 2016b). A scaling X = 4.4 occurs in the moderate scattering regime for a Kolmogorov medium (Cordes & Shannon 2010), though there are known sources with indices that deviate from such a scaling (e.g., Bhat et al. 2004;Löhmer et al. 2004). Here we take X = 4.4. The parameter estimate for t ∞ is given by the first element of θ = (X T C −1 X) −1 X T C −1 T. If we let T = T 0 + T C , where T 0 describes the components of the TOAs modeled by the traditional DM fit discussed previously and T C describes the systematic t C,ν , then the systematic error is simply the difference between the fit when only DM is accounted for and when the chromatic term is accounted for, i.e., We add this systematic error σ δt C in quadrature to the standard error σ DM . Again, we automatically account for stochastic changes in the PBF from DISS since scintilla-tion noise is included in C. Estimation and uncertainties on the shape of the profiles as a function of frequency will also enter into σ δt C but we have assumed perfect knowledge of the shape changes and do not include the effect here. As with σ DM , we provide the simplified form of σ δt C in Table 1 when only two frequency channels are used.
The frequency dependence of DM is described in Cordes et al. (2016) as a result of ray paths at different frequencies covering different volumes through the ISM. They primarily consider observations at two spot frequencies ν 1 and ν 2 ; observations covering greater frequency ratios r = ν 2 /ν 1 will have greater error in DM estimation since the difference in DM between the two frequencies will statistically grow larger. We combine Eqs. 15 and 25 of their paper to obtain a scaling for the rms timing error between the two frequencies for a uniform Kolmogorov medium, where E 11/3 (r) is a dimensionless quantity that contains all of the relative frequency dependence with respect to the type of medium. The covariance matrix for the total frequency-dependent DM error σ DM(ν) can be constructed simply with components σ DM(ν);pair (ν, ν ). Again, we add this error in quadrature to the previously discussed DM estimation errors, and thus the total DM estimation error (squared) is All terms depend on the center frequency ν 0 and bandwidth B but we drop them for clarity.

Telescope Error
Besides radiometer noise discussed previously, additional noise contributions from the telescope or its backends or receivers, σ tel , such as from RFI, will contribute to the overall TOA uncertainty. Narrow-band RFI such as from artificial satellites will contaminate particular frequency channels in pulsar data that can be excised, making the effective bandwidth observed smaller without introducing additional TOA uncertainties, which we do not account for in our analysis. However, any remaining RFI will increase the TOA uncertainty.
Pulse shape changes from incorrect absolute gain calibration and summation of the polarization profiles into intensity profiles lead to additional TOA uncertainties from typical template fitting (Lam et al. 2016a). For circularly polarized channels, the TOA uncertainty from gain variation is ; also see Stinebring et al. 1984 for an in-depth overview of alterations to the polarization properties of pulse profiles; here we consider Gaussian profiles propagated through the Mueller matrices in Appendix A of that paper) where ε = δg/g is the fractional gain error, π V is the degree of circular polarization, and W is the pulse width.
The fiducial values listed for the polarization parameters match closely to the NANOGrav AO data in the 1.4 GHz frequency band (ε ∼ 0.08, π V ∼ 0.1; P. A. Gentile et al. submitted). The uncertainty σ pol can be a systematic or random error depending on how ε varies in time. All three quantities are frequency dependent and therefore σ pol is also a function of frequency. Instrumental selfpolarization can also induce TOA uncertainties. Cross coupling for a circularly-polarized feed will result in a measured false circular polarization π V 2η 1/2 π L assuming the receiver polarizations are orthogonal, where η is the voltage cross-coupling coefficient and π L is the degree of linear polarization for the pulsar (see again Cordes et al. 2004). The false circular polarization yields a TOA error as in Eq. 26. Cross-coupling for our current set of receivers has not been well-measured and so we exclude the contribution to the TOA uncertainty in this work.
As with the various white-noise components to the timing uncertainty, the total uncertainty from polarization errors can be determined from the covariance matrix and using Eq. 3 to find the single value (squared) for the observation. Since telescope environments and instrumentation can vary wildly, we do not account for other telescope errors in our specific pulsar analysis though in general it will add into the total TOA uncertainty (Eq. 4).

Integration Time Considerations
Increased integration time decreases the white noise components of the TOA uncertainty; most NANOGrav pulsars are dominated by white noise and see large improvements in GW sensitivity by increased time per source ). Increased observing time improves the average S/N by T 1/2 as per Eq 10. The jitter error also improves by the same factor because the number of pulses averaged increases (Lam et al. 2016a). The scintillation error scales according to Eq 18 and 21, also as ∼T 1/2 when the number of scintles is large.
Increased integration time will also affect the number of scintles measured in a single observation. While the average observed flux across epochs remains the same (see Eq. 17), the distribution will change; in the lowscintle limit the median S/N is much lower than the mean and the majority of measured TOA uncertainties will be below the mean TOA uncertainty. Since many MSPs have white-noise TOA uncertainties dominated by template-fitting errors, ensuring that multiple scintles are observed every epoch is imperative for robust timing-model fitting. If the scintle size is of order the integration/bandwidth of the observation, the intensity will often be suppressed as per the PDF in Eq. 17. In the very low-S/N limit, determination of a TOA is not possible at all (again see Appendix B of Arzoumanian et al. 2015b).
In the case where multiple frequency bands are observed non-simultaneously within an epoch, an additional uncertainty due to time-variable DM must be accounted for in the DM estimation process (Lam et al. 2015). While typically a minor effect on the TOA uncertainty, the error grows when the time between observation increases, when ∆t d is small, or when the DM over The circle indicates the approximate ν 0 and B used in current observations at the GBT, with σ TOA (ν 0 = 1.3 GHz, B = 1.2 GHz) = 99 ns. However, because NANOGrav currently uses two observing bands separately to cover the whole frequency range, it takes 60 minutes of on-sky time to obtain 30 minutes of integration across the band; therefore, the overall σ TOA can be improved by a factor of ∼ √ 2 using a single wideband receiver with the same integration time.
the line-of-sight rapidly changes such as when a pulsar is observed close to the Sun (You et al. 2007;Howard et al. 2016;Jones et al. 2017). For optimal DM correction and reduced TOA uncertainty, observing with a wide band simultaneously is preferred.

IMPLICATIONS FOR NANOGRAV
In this section, we describe our analysis of five specific MSPs from the NANOGrav PTA. We considered pulsars observed with the two telescopes primarily used by NANOGrav: the Green Bank Telescope (GBT) at the Green Bank Observatory and the Arecibo Observatory (AO). Each telescope has multiple receivers to cover portions of the total frequency band, with observations at each band being nearly contiguous at AO but with multiple days of spacing at the GBT for observational efficiency. Most pulsars are observed over 2−3 frequency bands chosen to optimize the long-term timing rms and with a ∼30-day cadence though several are observed weekly (Arzoumanian et al. 2018).
Our goal was to calculate σ TOA (ν 0 , B) for a given set of pulsar-and telescope-specific parameters. Since we covered a large range in center frequencies and bandwidths, we chose to calculate over the logarithm of each quantity, i.e., σ TOA (log ν 0 , log B), between 0.1 and 10.0 GHz. We excluded bandwidths greater than 1.9 times the center frequency from the analysis and set the maximum σ TOA to be the pulsar spin periods. We assumed a constant gain of 2 K/Jy at GBT and 10 K/Jy at AO. To simplify our calculations, we chose sub-band bandwidths such that there would be multiple scintles in each sub-band; this is not a common practice for pulsar timing observations. Such a choice also ensures that the pulse S/N would be approximately S 0 and the scintillation bandwidth would be approximately constant over the sub-band. At the highest frequencies however, typically of order a few GHz, scattering transitions from the strong to the weak regime (Rickett 1990;Lorimer & Kramer 2012). In the weak scattering regime where DISS is attenuated by averaging effects, the PDF of the scintillation gain g will be log-normally distributed rather than be in the χ 2 -form given in Eq. 17 (J. M. Cordes et al. in prep). Therefore, even at high frequencies when we would previously consider n ISS ∼ 1, the pulse S/N will tend towards the mean S 0 . However, to simplify our analysis, we ignore larger epoch-to-epoch variations in the S/N, noting that TOA uncertainties will be better on certain epochs and worse for others; a fuller analysis will want to account for such a variation. We note that for some lines of sight such as those through the Galactic plane, scattering will remain in the strong regime beyond our frequency limit of 10 GHz, with the largest effects in directions close to the Galactic center.
Since scintillation noise also depends on the number of scintles in the band but changes between the strong and weak scattering regime, we approximate D as follows. When there are multiple scintles (n ISS ≥ 2) in the frequency sub-band, we assume that the sub-bands are uncorrelated with each other. In the other case, we assume that the sub-bands are completely correlated (ρ DISS = 1) and therefore, the scintillation covariance Gain (K/Jy) 2/10 g 2 10 2 2 Ae = 27600 for AO, 5520 for GBT a Values provided are referenced to 1 GHz. The Galactic background temperature is set to scale proportional to ν −2.75 . We set Tconst = 30 K, ε = 0.08, and π V = 0.1. We used an integration time of T = 30 minutes. b Scaled from Keith et al. (2013). c Estimated from Cordes & Lazio (2002). d Scaled from Levin et al. (2016) assuming C 1 = 1.16. e Scaled from Champion et al. (2008). f For an approximation of t C , we used scattering timescale variations taken as half the maximum range in Levin et al. (2016), otherwise we assumed δτ d,0 = 0.5τ d,0 . See text for more details.
g For PSR J1713+0747, we list GBT and then AO parameters when two numbers are provided.
matrix will take the block-diagonal form (28) For our analysis, we chose 100 sub-bands for each calculation of σ TOA (ν 0 , B).
Without prior knowledge of the amplitude of the extra chromatic delay t C (ν), we used half the total variation (maximum minus minimum) in τ d,0 given in Levin et al. (2016) as an estimate for the amplitude of the chromatic delay away from the nominal scattering timescale. We set the frequency dependence proportional to a fiducial value of ν −4.4 , allowing us to include this additional systematic uncertainty in our analysis. In cases where the variation was not measured, we used δτ d,0 = 0.5τ d,0 , which is roughly the amplitude for pulsars with measured variations. Table 2 lists the parameters we included and the values we adopted for our selection of MSPs. For simplicity, we assumed that the telescope parameters were constant over all frequencies. In addition, we chose W eff , W 50 , and σ J to be constant over all frequencies, using the 1.4 GHz values in Lam et al. (2016a). In reality, the values do vary across frequency but we do not have information at the lowest and highest frequencies in our analysis. Interpolating values over the frequency range we do have information for does not change our results substantially. In the next subsections, we discuss the results for the pulsars individually. The code that we used in our analysis is freely available at https://github.com/mtlam/FrequencyOptimizer 5 .

PSR J1713+0747
The precision timing for PSR J1713+0747 has led to its usage in a number of test of gravity and it typically dominates the sensitivity of PTAs towards gravitational waves (e.g., Zhu et al. 2015;Arzoumanian et al. 2015a). Current observations of the pulsar by NANOGrav are conducted at both the GBT and AO with a weekly cadence because of its high precision timing. We therefore perform our calculations for this pulsar at both telescopes but separately. Figure 1 shows the results of our optimal frequency and bandwidth analysis for PSR J1713+0747 for 30-minute observations at the GBT. The solid contours represent TOA uncertainties of 2, 1, 0.5, 0.2, and 0.1 µs, in order of increasing darkness and thickness. The two dashed blue contours show the 10% and 50% increase above the minimum σ TOA . The hatched region represents where B > 1.9ν 0 . We remind the reader that in our calculations for systematic delays, we only correct for a DM term and no scattering estimation is performed, thus mimicking current observing schemes. The increase in TOA uncertainty at the highest bandwidths comes partially from the lack of scattering correction (though a full scattering correction still yields a scattering misestimation term, for more see Figure 9 of Cordes & Shannon 2010) and partially from an increased contribution from frequencydependent DM. If epoch-to-epoch DM variations were not removed, the TOA uncertainties will drastically increase towards larger bandwidths as well. The black star in Figure 1 shows the minimum TOA uncertainty is σ TOA (ν 0 = 9.1 GHz, B = 16.6 GHz) = 85 ns. The black circle shows NANOGrav's current frequency coverage assuming a continuous frequency range from ∼0.7 to ∼1.9 GHz, i.e., a center frequency of ∼1.3 GHz and a bandwidth of ∼1.2 GHz. In reality, there is a slight gap in frequency coverage from ∼0.9 to 1.1 GHz and high cadence (∼ weekly instead of monthly) observations only cover the higher frequencies above the gap (Arzoumanian et al. 2018), thus complicating the analysis. Ignoring gaps, the TOA uncertainty is σ TOA (ν 0 = 1.3 GHz, B = 1.2 GHz) = 99 ns.
The above analysis assumes the standard NANOGrav receiver setup, which requires 60 minutes of on-sky time in order to obtain 30 minutes of effective integration time uniformly across the entire band. Were a broadband receiver available, capable of covering the entire 0.7 GHz to 1.9 GHz band, the overall σ TOA could be improved by a factor of ∼ √ 2 using a single wideband receiver with the same effective integration time. We see that for this pulsar and observing setup, we lie within the phase-space valley very close to the minimum for the average S/N epoch though some slight improvements can be made by going to higher frequencies. Figure 2 shows the same analysis but completed for AO. We see that because of the sensitivity of AO, the same TOA uncertainty covers a larger frequencybandwidth region. An increase in the gain improves the pulse S/N overall, which allows for reduced TOA uncertainties at higher frequencies where the pulsar powerlaw spectral flux is diminished but also at lower fre- quencies where scattering begins to reduce the flux as well; our analysis demonstrates why considering every contribution to the TOA uncertainty is important. For PSR J1713+0747 observed with AO, the TOA uncertainty at the minimum is σ TOA (ν 0 = 10.0 GHz, B = 17.8 GHz) = 63 ns. At AO, NANOGrav observes PSR J1713+0747 between ∼1.2 and ∼2.4 GHz, i.e., ν 0 = 1.8 GHz and B = 1.2 GHz. The estimate of the TOA uncertainty for these parameters is σ TOA (ν 0 = 1.8 GHz, B = 1.2 GHz) = 66 ns. We note that since NANOGrav observes PSR J1713+0747 with both telescopes, determining the combined optimal frequency ranges becomes more complex since the telescope gains and T const values will be different, though the TOA uncertainty will be weighted more towards AO's improved sensitivity.
3.2. PSR J1909−3744 PSR J1909−3744 is one of the best-timed pulsars. It has the second lowest amount of white noise and the lowest measured excess noise beyond the white noise . Therefore it is one of the most sensitive components in PTAs to low-frequency gravitational wave detection (e.g., Arzoumanian et al. 2014;Shannon et al. 2015). Figure 3 shows the analysis for observations of PSR J1909−3744 using the GBT. Again assuming the GBT parameters of a center frequency of ∼1.3 GHz and a bandwidth of ∼1.2 GHz, we see that we obtain similar results to PSR J1713+0747 when observed with the GBT, although note the narrowing of the optimal region in the phase space. As reported for PSR J1713+0747, the minimum TOA uncertainty for PSR J1909−3744 is σ TOA (ν 0 = 8.1 GHz, B = 14.8 GHz) = 50 ns and the estimate given the current frequency coverage is σ TOA (ν 0 = 1.3 GHz, B = 1.2 GHz) = 59 ns. PSR J1903+0327 is the pulsar with the highest DM (297.52 pc cm −3 ) currently observed by NANOGrav. It shows considerable long-term red noise, both frequencydependent and independent, suggesting unmodeled ISM and possibly intrinsic effects . Figure 4 shows the analysis for observations of PSR J1903+0327 using AO. Like PSR J1713+0747 observed with AO, NANOGrav observes this pulsar with ν 0 = 1.8 GHz and B = 1.2 GHz, denoted by the black dot.
The minimum TOA uncertainty is σ TOA (ν 0 = 9.8 GHz, B = 13.2 GHz) = 1.0 µs and the estimate given the current frequency coverage is σ TOA (ν 0 = 1.8 GHz, B = 1.2 GHz) = 44.0 µs. Note that since the Galactic latitude of the pulsar is ≈ −1 • , the pulsar emission remains in the strong-scattering regime as discussed at the beginning of the section. Therefore, the S/N at high frequencies will be reduced for a majority of the epochs and the observed pulses will be scintillated high a small fraction of the time. For a typical epoch, the reduction in S/N implies that the center frequency will be pushed to lower frequencies where the S/N is higher. More work is needed to determine the effect of epoch-toepoch variation in the scintillation gain on the long-term timing rms.
We see that the current projected TOA uncertainty is quite poor for this pulsar and an order of magnitude larger than measured in previous work Arzoumanian et al. 2018). The dominant factor for this discrepancy is the unknown scattering timescale variation parameter δτ d,0 . Since the variations (and the value itself) were unmeasured in Levin et al. (2016) and we have adopted a fiducial value of 50% of τ d,0 when we cannot estimate them, we arrive at such a large number for σ TOA . If we reduce the variation to 10%, the contours shift to the left some, such that we match with justify using such a low value. However, while we see that there is a large discrepancy, the overall trend in the phase space remains quite clear. For this pulsar, observing at higher frequencies and bandwidths will improve its timing quality. Scattering corrections may be able to improve the overall timing of this pulsar even further, as well as others with high DM (Shannon & Cordes 2017).
3.4. PSR J1744−1134 On the other extreme, PSR J1744−1134 is the pulsar with the lowest DM (3.09 pc cm −3 ) currently observed by NANOGrav. We show it as an example of a "typical" NANOGrav MSP with TOAs spanning the full length of our observing campaign and with scintillation parameters and rms timing residuals (TOAs minus timing model) comparable to many other NANOGrav pulsars. It does show frequency-independent excess white noise beyond the white noise in its timing residuals though no detectable red noise over the first eleven years of observation Arzoumanian et al. 2018). Figure 5 shows the analysis for observations of PSR J1744−1134 using the GBT. The frequency coverage is the same as for PSR J1909−3744. Note that at the highest bandwidths, we see the best TOA uncertainties for a given center frequency because the various delay misestimation terms do not dominate for a pulsar with such a low DM unlike in the previous examples. Therefore, observing with the highest possible bandwidths becomes important for pulsars similar to PSR J1744−1134, many of which are currently being observed in PTA efforts. The minimum TOA uncertainty is σ TOA (ν 0 = 8.3 GHz, B = 15.5 GHz) = 120 ns and the estimate given the current frequency coverage is σ TOA (ν 0 = 1.3 GHz, B = 1.2 GHz) = 140 ns.
3.5. PSR J1643−1224 PSR J1643−1224 is a high-DM pulsar (62.41 pc cm −3 ) with evidence of red noise in the timing residuals possi-bly from uncorrected ISM effects Lam et al. 2017;Arzoumanian et al. 2018); the pulsar lies beyond the HII region Sharpless 2-27 (Gum 73) and thus plasma along the line of sight may be affecting the timing noise along with the determination of pulsar parameters (Blitz et al. 1982;van Leeuwen 2007;Matthews et al. 2016;Reardon et al. 2016). Corrections of interstellar scattering delays can improve the timing of the pulsar (Lentati et al. 2017). The pulsar has also shown variations in the pulse profile interpreted as a change in the magnetosphere and such variations can also cause deviations from the timing model . Figure 6 shows the analysis for observations of PSR J1643−1224 using the GBT. The frequency coverage is the same as the other GBT pulsars. As with PSR J1903+0327, at the current center frequency, NANOGrav's timing noise is sub-optimal and would benefit from an increase in the center frequency. We stress that observation of pulse profiles at large bandwidths allows for estimation of other parameters, such as the scintillation timescale and bandwidth, and therefore it is important to collect as much data as possible and then choose what data will go into a timing analysis. As stated previously though, increases in the center frequency and the bandwidth jointly will lead to an improvement in σ TOA ; such increases become important for pulsars like PSR J1643−1224 with moderate DMs and chromatic noise. For PSR J1643−1224, the minimum TOA uncertainty is σ TOA (ν 0 = 10.0 GHz, B = 16.6 GHz) = 1.3 µs (recall that we strict our analysis to ν 0 ≤ 10.0 GHz) and the estimate given the current frequency coverage is σ TOA (ν 0 = 1.3 GHz, B = 1.2 GHz) = 20.8 µs. The unknown scattering timescale variation δτ d,0 for this high DM pulsar is again the biggest contributor to the difference between the estimated current σ TOA and the observed values in the literature.

BROAD IMPLICATIONS FOR PRECISION PULSAR TIMING
We have seen in the previous section our frequency optimization analysis for five example pulsars. Now we extend our analysis broadly to more MSPs. Since the pulsar flux I 0 and DM are the most dominant parameters, we chose to vary these parameters while holding all other parameters fixed for a fiducial pulsar.
We assumed a spectral index α = −1.6, a constant effective pulse width W eff = 500 µs,Ū obs = 20, and rms jitter σ J,30 min = 100 ns. We also assumed that the observations are performed with AO. We saw slightly lower optimal center frequencies when using the GBT gain of 2 K/Jy but do not show the plots since they are qualitatively similar. Note this trend is different from our results for PSR J1713+0747, though we discuss the possible reasons shortly. For the scattering timescale, we used the scaling from Bhat et al. (2004), log τ d,ms = −6.46 + 0.154 log DM + 1.07(log DM) 2 −(3.86 ± 0.16) log ν GHz , to vary τ d and assumed δτ d,0 = 0.5τ d,0 , allowing us to input a chromatic term t C into the analysis. The maximum-to-minimum variation in τ d found in Levin et al. (2016) was typically found to be of this order. Assuming a uniform Kolmogorov medium, the scintillation Optimal ν 0 (GHz) Figure 7. Optimal center frequency ν 0 as a function of flux density at 1 GHz I 0 and DM for pulsar parameters discussed in §4. Here we allow for any possible bandwidth. The black region denotes where ν 0 ≥ 10 GHz and is thus excluded from our analysis. The bands in the middle of the plot are a result of the minimum σ TOA moving near our bandwidth cutoff B = 1.9ν 0 . The minimum optimal ν 0 of this phase space is around 1.9 GHz. The "ripple" feature around the highest frequencies results from the minimum σ TOA wandering around the B = 1.9ν 0 cutoff. The feature in the top left represents slight numerical wandering of the minimum σ TOA in a large valley in the ν 0 , B phase space and is over-represented in area in this figure due to the logarithmic axes.
timescale is (Cordes & Rickett 1998) where V p,eff⊥ is the pulsar's effective velocity perpendicular to the line of sight, which includes components from the pulsar and Earth velocities relative to their respective Galactic rotating frame and the random ISM velocity. Assuming the pulsar velocity is the dominant component, we can convert the equation to use the total proper motion, µ T . We assumed that D = 1 kpc and µ T = 10 mas yr −1 as fiducial values (Matthews et al. 2016). Figure 7 shows the optimal center frequency for our fiducial pulsar as a function of I 0 and DM allowing for any possible bandwidth. Pulsars with low DM, such as PSR J1744−1134, obtain better σ TOA when observed at lower center frequencies. Conversely, pulsars with high DM, such as PSR J1903+0327, are dominated by ISM-related timing effects, which are partially mitigated by the higher center frequencies. However, all ν 0 were 1.9 GHz, indicating the preference to observe with higher center frequencies than typical current observing configurations.
Especially at larger optimal frequencies, we found typically that the optimal bandwidths were near or on the maximum of our allowed range (B ≤ 1.9ν 0 ). Figure 8 shows the optimal center frequency where we considered more realistic engineering constraints for a single receiver, with B = 0.5ν 0 up to a maximum bandwidth Optimal ν 0 (GHz) Figure 8. Optimal center frequency ν 0 as a function of flux density at 1 GHz I 0 and DM for pulsar parameters discussed in §4.
Here we set the maximum allowed bandwidth to be half the center frequency (B = 0.5ν 0 ) up to a bandwidth of 4 GHz. The black region denotes where ν 0 ≥ 6 GHz for improved dynamic range in the plot. The minimum optimal ν 0 of this phase space is around 0.48 GHz. As in Figure 7, the "ripple" now in the top left represents the minimum σ TOA wandering near the bandwidth cutoffs.
of 4 GHz. We see that because of the constraints, lower optimal frequencies are preferred and are limited by the bandwidth restrictions. We stress that the bandwidth restrictions are not hard limitations but serve as approximate guidelines that may be improved upon in the future. As is current procedure, multiple telescopes and/or receivers can be used to observe a source within a single epoch and therefore the optimal bandwidths can be larger than what a single receiver will allow. Increases in telescope collecting area will improve the pulse S/N and therefore fainter pulsars can be observed with similar σ TOA as brighter pulsars. Since the S/N of pulses will improve with collecting area (i.e., the telescope gain), larger telescopes will have optimal center frequencies that may be preferentially higher than their smaller counterparts for objects with higher DMs but may in fact be lower for local objects with very low scattering. For any pulsar, observations that combine data from a variety of different telescopes can be tuned to cover the optimal frequency range for a given timing analysis and therefore provide the highest TOA precision possible.

DISCUSSION
For the first time, we have demonstrated a comprehensive methodology for calculating TOA uncertainties as a function of center frequency and bandwidth and therefore have determined the optimal frequencies for a pulsar to be observed given a set of pulsar-, Galactic-, and telescope-dependent parameters. So far we have only discussed the uncertainty on a TOA for a single epoch; an analysis of the long-term rms of the timing residuals must include additional noise terms that have not yet been fully accounted for, including systematic timing changes from temporal variations in parameters such as t C . Our optimization analysis is very generic in its applications and can be used in any study that requires high-precision pulsar timing, including: tests of gravity, experiments in nuclear physics, determination of astrometry and stellar dynamics, detection of gravitational waves, probes of the ISM, and others.
Looking at the NANOGrav pulsars examined here, we conclude that many MSPs timed in PTA observing programs (e.g, the European Pulsar Timing Array and Parkes Pulsar Timing Array collaborations, see Desvignes et al. 2016 andReardon et al. 2016, respectively) are done so with non-optimal frequency-bandwidth combinations. The range of optimality is much larger when observations are radiometer-noise-limited however. For many MSPs, especially those at higher DMs, PTA observations would highly benefit from wideband systems typically centered at somewhat higher observing frequencies than the receiver combinations used so far. Very large increases in bandwidth can also help with ancillary data products that can provide additional information into the overall timing and noise model (e.g., in estimating perepoch scintillation parameters); we do not suggest that wide-bandwidth systems be limited in size but rather that frequency ranges for should be carefully chosen for optimal TOA estimation specifically. Wide bandwidth systems can also help improve how optimal center frequencies are chosen and reduce the necessary integration time per pulsar used by current multiple-band systems, providing an instant reduction in the TOA uncertainty by ∼ √ 2, assuming a fixed integration time. While we have discussed how to optimize the frequency ranges for a given pulsar, optimization of an entire PTA requires consideration of additional factors. In general, the specific scientific goals will dictate the optimization of the array. For example, continuous gravitational wave sources benefit from observations of a few well-timed pulsars (Arzoumanian et al. 2014). A stochastic background of gravitational waves from sources such as supermassive black hole binaries or the inflationary epoch require many pulsars distributed across the sky (Arzoumanian et al. 2016). In all cases, however, even small increases in TOA precision over many pulsars and many epochs can make a large difference in the detection and eventual characterization of GW sources.
In practice, implemented receiving systems tend to show a decrease in sensitivity (increase in system temperature) as the bandwidth increases above an octave. This penalty has been only poorly quantified, however, and we have not tried to represent it in this work. However, this penalty would likely shift the optimal frequency slightly down and to the left in the σ TOA (ν 0 , B) phase space if a larger than octave bandwidth is implemented. As new receivers and new telescopes become available, our formalism should be considered for the development of radio instrumentation used for high-precision pulsar timing.
was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.