The Green Bank North Celestial Cap Survey. IX. Timing Follow-up for 128 Pulsars

The Green Bank North Celestial Cap survey is one of the largest and most sensitive searches for pulsars and transient radio objects. Observations for the survey have finished; priorities have shifted toward long-term monitoring of its discoveries. In this study, we have developed a pipeline to handle large data sets of archival observations and connect them to recent, high-cadence observations taken using the Canadian Hydrogen Intensity Mapping Experiment telescope. This pipeline handles data for 128 pulsars and has produced measurements of spin, positional, and orbital parameters that connect data over observation gaps as large as 2000 days. We have also measured glitches in the timing residuals for five of the pulsars included and proper motion for 19 sources (13 new). We include updates to orbital parameters for 19 pulsars, including nine previously unpublished binaries. For two of these binaries, we provide updated measurements of post-Keplerian binary parameters, which result in much more precise estimates of the total masses of both systems. For PSR J0509+3801, the much improved measurement of the Einstein delay yields much improved mass measurements for the pulsar and its companion, 1.399(6) M ⊙ and 1.412(6) M ⊙, respectively. For this system, we have also obtained a measurement of the orbital decay due to the emission of gravitational waves, ṖB=−1.37(7)×10−12 , which is in agreement with the rate predicted by general relativity for these masses.


Introduction
Since the discovery of pulsars in the late 1960s (Hewish et al. 1968), they have been the subjects of intense study.As these stars rotate, their magnetic fields beam radio waves from their poles like a lighthouse.On Earth, this rotation is detectable as a series of pulses.Precise measurements of these pulses at Earth track the passage of time in the pulsar's frame, and the timing models compare this clock to those on Earth.Deviations from the predicted arrival times of pulses encode information about the pulsar and material along the line of sight to it in a process called "pulsar timing" (e.g., Lorimer & Kramer 2012).
The astrophysical applicability of pulsars is wide (e.g., Backer 1975Backer , 1984;;Hewish 1975;Detweiler 1979;Taylor & Weisberg 1982;Taylor et al. 1992).For instance, pulsar astronomy has provided some of the most stringent constraints on the behavior of ultra-dense matter (Özel & Freire 2016;Bogdanov 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.et al. 2019;Hu et al. 2020;Lattimer 2021;Pang et al. 2021).Their immense moments of inertia make them some of the most stable clocks in the Universe (Manchester 2017;Yin et al. 2017).The most stable millisecond pulsars (MSPs) can be used as an ensemble to search for low-frequency gravitational waves (Detweiler 1979;Desvignes et al. 2016;Reardon et al. 2016;Arzoumanian et al. 2020;Antoniadis et al. 2022;Agazie et al. 2023aAgazie et al. , 2023bAgazie et al. , 2023cAgazie et al. , 2023d, 2023e;, 2023e;Afzal et al. 2023).
All these topics benefit from ongoing searches for pulsars, which continue to increase the known population.These searches implement ever-improving technology and searching algorithms that have discovered over 3000 pulsars to date (e.g., Manchester et al. 1978Manchester et al. , 2001Manchester et al. , 2005;;Deneva et al. 2013;Stovall et al. 2014;Sanidas et al. 2019;Cruces et al. 2021;Sengar et al. 2023, etc.).Many of these pulsars have been timed to high precision, and have been used to characterize the pulsar population (Faucher-Giguère & Kaspi 2006;Lorimer et al. 2006;Levin et al. 2013;Lorimer et al. 2015).
This timing process is hampered by many components of pulsar evolution.Young pulsars are born following the supernova of massive progenitors; these highly magnetized neutron stars rotate the least stably and slow the most quickly (Antonelli et al. 2023).This slowing is due at least in part to magnetic dipole radiation, but typically does not follow so simple a model.Many pulsars, especially young pulsars, exhibit some form of timing noise (Parthasarathy et al. 2019;Lower et al. 2020;Singha et al. 2022).This can manifest as a long-term, random drift in the observed spin of a pulsar that deviates from "pure" spin-down due to magnetic dipole radiation (a power law in spin with an index of 3).Measurements have shown that this model is very rarely sufficient, suggesting a more complex relationship between these parameters (Lower et al. 2020).More generally, the adolescent pulsar's spin will follow a random walk in ν and  n.From the timing side, this noise is often removed by fitting and subtracting polynomials from the residuals-a process known as "polynomial whitening" (Hobbs et al. 2004).
Along with this drift, some pulsars show abrupt changes in spin period called "glitches" (Espinoza et al. 2011;Fuentes et al. 2017;Zhou et al. 2022aZhou et al. , 2022b;;Basu et al. 2022).Theories to explain these events largely utilize either of the following: 1.A starquake model where the crust of the neutron star abruptly changes shape and the change in its moment of inertia drives the observed momentum-conserving change in rotation (Alpar et al. 1994;Lai et al. 2018;Lu et al. 2023).2. Differential rotation of the solid crust and a superfluid interior that is coupled by the pinning of vortices in this superfluid (Alpar 1977;Link et al. 1992;Haskell et al. 2020;Layek et al. 2023;Melatos & Millhouse 2023).
While the process that causes pulsars to glitch has still not yet been definitively identified, their impacts on timing have been addressed in many studies (e.g., Lower et al. 2021;Dunn et al. 2022;Zubieta et al. 2023).In a simple case, the glitch appears in timing residuals as an instantaneous change in the slope.This change is also accompanied by a change in the spin-down rate of the pulsar, which introduces an additional accumulation of phase offset in the residuals.Many observed glitches are then followed by an exponential recovery of the affected parameters back to their pre-glitch values.All of these components manifest as changes in pulsar spin parameters, and complicate measurements of spin/spin-down.
On the other hand, the rotation of old pulsars is typically much more stable, particularly those that reside in binary systems.In many cases, these pulsars will harvest rotational momentum from the orbit via the Roche-lobe overflow of the companion.This process decreases the rotation period of the pulsar, potentially down to the millisecond level.It also greatly dampens their magnetic fields, reducing the long-term spindown rate from dipole radiation (Bisnovatyi-Kogan & Komberg 1974;Shibazaki et al. 1989).These pulsars, called millisecond pulsars (or MSPs), have been utilized for the majority of highprecision pulsar science (Alam et al. 2021;Arzoumanian et al. 2023;Falxa et al. 2023;Miles et al. 2023).Aside from their stability, the timing procedure is sensitive to Doppler shifts due to binary motion.These shifts can be used to determine three of the Keplerian parameters of the orbit, including the orbital period, the pulsar's semimajor axis, and the initial phase of the binary (relative to the time of measurement).These can be used to place limits on the masses in the system via the Keplerian mass function.Some binary orbits are compact enough that the pulsar reaches relativistic speeds during the orbit, introducing additional complexities in the timing procedure that can be used to directly measure the mass of the pulsar and its companion (Demorest et al. 2010;Özel & Freire 2016;Cromartie et al. 2020).Given that this measurement is nearly impossible to make otherwise, binary pulsars are unique and powerful laboratories for gravitation and ultra-dense matter; and in those few sources where the masses can be measured via detection of multiple independent post-Keplerian parameters (as for the double pulsar system PSR J0737−3039 and the triple system PSR J0337+1715), stringent constraints are placed on the underlying theory (Archibald et al. 2018;Voisin et al. 2020;Kramer et al. 2021).
Precise measurements of these parameters is hampered by infrequent sampling, as models that describe pulse arrival times well over a single observation may fail to do so for subsequent observations when there are long wait times between scans.The development of high-cadence radio instruments like the Canadian Hydrogen Intensity Mapping Experiment (CHIME; CHIME/Pulsar Collaboration et al. 2021), which observe large swathes of the sky on a daily basis, helps to avoid this loss of phase connection.In this way, CHIME has already proven itself a very capable pulsar-timing instrument, especially with sources that are difficult to time with episodic timing campaigns (Fonseca et al. 2021;Good et al. 2021;Dong et al. 2023).
Here, we discuss the continued efforts of a pulsar survey that has reached the end of its observing program: the Green Bank North Celestial Cap pulsar survey (GBNCC).We discuss the survey's conclusion and current objectives in Section 1.1.In Section 2, we describe the sample of pulsars we have examined in this study and how we produce data products, including the addition of high-cadence CHIME observations.Section 3 outlines the procedure used to generate timing data products and use them to refine models.Section 4 covers the results from our study and highlights some particularly interesting measurements.Finally, Section 5 summarizes our results.We also include an Appendix with timing residuals for all sources.The work below includes comparisons with the ATNF Pulsar Catalog (v1.70;Manchester et al. 2005).28

GBNCC Survey Overview and Completion
Over the last decade, the GBNCC survey has covered the sky north of δ = −40°with 124,852 pointings of 2 minutes duration with the 100 m Green Bank Telescope (GBT) at 350 MHz.Full specifications of survey observations and processing are detailed in previous GBNCC publications (Stovall et al. 2014;Lynch et al. 2018b).This choice of central frequency provides additional sensitivity to steep spectrum sources, in particular older pulsars outside of the Galactic plane (McEwen et al. 2020).The survey is one of the largest pulsar surveys in terms of sky coverage, covering 36,430 deg 2 .
To date, the GBNCC survey has discovered 195 sources, including 33 MSPs, 24 binaries, 24 rotating radio transients (McLaughlin et al. 2006), and a fast radio burst (FRB; Lorimer 2018), many of which have been published in a series of survey papers (Stovall et al. 2014;Lynch et al. 2018b;Kawash et al. 2018;Aloisi et al. 2019;Parent et al. 2020;Agazie et al. 2021;Fiore et al. 2023;Swiggum et al. 2023).These discoveries have included some of pulsar astronomy's most exotic sources and important measurements, including the high mass of PSR J0740 +6620 (Cromartie et al. 2020;Fonseca et al. 2021), eclipsing black widow pulsars that are ablating their companions (Swiggum et al. 2023), nulling pulsars (Anumarlapudi et al. 2023), and double neutron star binaries (Lynch et al. 2018b;Aloisi et al. 2019;Swiggum et al. 2023).These studies and others have amassed data on GBNCC discoveries, much of which is currently available from the public-facing GBNCC GitHub repository. 29uring the summer of 2022, the final GBNCC survey positions were observed (aside from 500 points which will be reobserved for reasons related to radio frequency interference (RFI) and observation complications).Efforts to identify new pulsars are continuing, both with our existing pipeline and with new pipelines that improve sensitivity to certain regions of parameter space (Sengar et al. 2023).For this study, we instead focus on the continued timing of prior discoveries by connecting archival data with new observations.Note that there are also 12 sources included in this paper that were discovered in the GBT350 survey of the northern Galactic plane (Hessels et al. 2008).They were followed up at 820 MHz using the GBT during an earlier timing campaign (project code AGBT13B_290; PI: R. Rosen), and are included among the CHIME/Pulsar sources (for details about the CHIME/Pulsar system, see CHIME/Pulsar Collaboration et al. 2021).

Observations
Following identification, candidates were observed in test scans using either the GBT (earlier discoveries) or CHIME/ Pulsar (starting in early 2020).In the former case, two to three 10 minute scans were taken over the course of a few weeks during other survey observations to confirm discoveries.New pulsars were then included in follow-up proposals with the GBT to establish phase-connected timing solutions.The observing specifications of these proposals vary slightly to account for complex timing models, but typically included ;10 follow-up observations (10-20 minutes on source per scan) with the GBT using either the 350, 820, or 1400 MHz receiver.Depending on the success of the timing campaign and potential scientific benefits, some pulsars were included in subsequent proposals using other instruments like the Low Frequency Array (LOFAR; van Haarlem et al. 2013) and Arecibo.
Following an agreement between the CHIME/Pulsar and GBNCC collaborations in late 2019, candidates (and pulsars) identified in GBNCC data could be followed up using the CHIME/Pulsar instrument with minimal latency due to its continuous and commensal manner of operation.During the spring semester of 2020, many previously discovered (and in some cases, previously published) GBNCC pulsars were observed using the 600 MHz receiver on CHIME with cadences of 1-7 days.The scan length depends on the declination of the source, as pulsars close to the southern horizon are visible for a shorter time compared to those directly overhead.At a minimum (i.e., δ  0°), scans were 10 minutes; the longest scans lasted approximately 1 hr.Despite the slight decrease in sensitivity between the GBT versus CHIME, the high cadence of these observations provides unprecedented sensitivity to the pulsar behavior, and models can be updated on a regular basis.The reduction in beam size further constrains the error in position for the sources, particularly when gridding observations were taken.Normal 350 MHz survey observations continued during the intervening years, and new candidates were confirmed/rejected using both GBT test scans and CHIME/Pulsar daily observations.As the CHIME/Pulsar schedule filled, the priority (and cadence) of some sources was reduced to accommodate those with incomplete/exotic models.Eight pulsars discovered with declinations below CHIME's horizon (−10°) were also included in our timing analysis.These sources were observed solely with the GBT in dedicated timing campaigns.

GBNCC/CHIME Pipeline
As part of the observing activities automatically undertaken by CHIME/Pulsar, a simplified data-reduction and pulse time-ofarrival (TOA) extraction procedure is applied to each observation in order to assess data quality and generate initial timing data products.The procedure utilizes a Python-based workflow that interacts with the data via the standard PSRCHIVE (van Straten et al. 2012) Python interface.The automated workflow, running on the Cedar supercomputing cluster,30 will identify any new observation data products that were transferred from the observatory to the storage system.Subsequently, standardized RFI mitigation is performed on the data using clfd (Morello et al. 2019), including the excision of known corrupted frequency channels for that observation obtained via CHIME/Pulsar utilities.Two reduced copies of the cleaned archive are then retained, one with 32 subbands and 1 minute subintegrations, and another one with 1024 frequency channels and fully averaged in time.
Using daily CHIME/Pulsar observations greatly increased the data volume for many GBNCC sources.To best utilize these data, a pipeline was established between the data taking and data processing via GitHub.As new data are taken, singleepoch and frequency-averaged TOAs are generated using preliminary timing models and a detection pulse template.These TOAs are then uploaded in batches to a shared private GitHub repository, where they can be accessed by GBNCC.With the new data, timing models are updated using fitting tools in PINT (Luo et al. 2021).When a solution has been updated, it can then be passed back to CHIME/Pulsar for future observations as necessary.In some cases, new timing models change the pulse profile enough to justify updating the TOA template, and therefore requires re-extraction of TOAs from all previous observations.Phase-connected solutions offer positional precision well within the observing beam size (∼0.5°at 400 MHz, ∼0.25°at 800 MHz), particularly after a year of observing.

TOA Excision
The data set contains 78,585 TOAs from four telescopes (GBT, Arecibo, LOFAR, and CHIME).As the timing models develop, measured TOA uncertainties will reduce as the profile is improved.In some cases, position improvements from phase connection dramatically improve the signal-to-noise ratio (S/N) of detections.In others, individual CHIME observations resulted in low-S/N detections, even with a timing-derived position.In these cases, multiple observations were summed together before producing TOAs, resulting in fewer (but more precise) TOAs.The wealth of data coming from CHIME allows us to be more selective with TOAs, and we implement TOA zapping as a part of the timing pipeline.Many of the zapped TOAs are due to faint detections and excess RFI.There are also a few epochs of known TOA issues which have been removed more broadly from our data set.As discussed in Andersen et al. (2023), improper packaging of CHIME/Pulsar data prior to MJD 58550 resulted in corrupted data.These TOAs were zapped from all of the sources.

Timing Procedure
The timing baselines for pulsars included in this data set are highly source dependent.A large number of sources (95) include archival data (pre-2019) from the GBT and other telescopes, and nearly all (120) sources have data taken with CHIME (post-2019).There are eight pulsars not visible with CHIME, and so solutions were determined from GBT observations only.For this project, data are only included up to 2023 May, and a minimum baseline of 1 yr was imposed.
After excising corrupted data, predicted pulse TOAs were compared with the time stamps of observed pulses.Correlated deviations indicate errors in the timing model that are fit using a least-squares algorithm.Determining the appropriate parameters to include in this fit is difficult, as covariance is high for any parameters that influence the residuals on timescales larger than the baseline.The covariance of sinusoidal (position, velocity, and orbital motion) and quadratic (rotation rate and deceleration) components limits the prediction power of timing models when the baseline is short.In many cases, the gap between GBT timing campaigns and CHIME/Pulsar daily scans is >2000 days.Fitting models over this gap can be very constraining, but requires a long enough baseline on one side of the gap to reduce phase uncertainty at the other side to within a single rotation.
When starting a timing solution from scratch, a preliminary solution with discovery parameters is used to fit TOAs over a short time period (typically a single observation).This fit only constrains the rotation rate of the pulsar as it is the only parameter with significant influence on short timescales.The improved frequency is then used to model TOAs at the MJD of the nearest observation, and a fit to both days provides a more precise measurement.The regular cadence of CHIME/Pulsar observations can introduce a subtle aliasing error into the model, as the pulsar rotation can become degenerate with the observation spacing determined by sidereal rate (see in CHIME/FRB Collaboration et al 2020).To mitigate this, single CHIME observations can be split into multiple TOAs to find an appropriate starting model to fit to subsequent days.This procedure is repeated until the model no longer appropriately fits the data and additional parameters must be included.When multiband data exist for a pulsar, additional constraints are set on its dispersion measure (DM).This constraint benefits from CHIME's wide observing band.Errors in DM lead to a smeared pulse and reduced S/N, so implementing new DM measurements further improves TOA precision.
Included among the pulsars in this project are 19 binary pulsars.These sources require special attention in the timing process, especially when the orbital period is 1 yr.For short orbital periods (where the observation duration t obs spans 1 orbital period P B ), orbital parameters can be constrained almost immediately afterwards.When the orbital period is long (P B 1 yr), covariance with Earth's motion and the pulsar's intrinsic spin-down become significant.For these cases, fitting timing residuals with polynomial functions of order 2 captures the apparent spin derivatives due to the orbital motion, and can be used to approximate the orbital parameters.Ultimately, these derivatives are replaced by orbital parameters when possible.

Results
We have assembled timing data for 128 sources; timing residuals for all sources are given in the Appendix and a periodperiod derivative plot is shown in Figure 1.Many of the sources included have had detailed timing analyses published in previous/ upcoming GBNCC publications.For these sources, we focus on maintaining phase coherence over long baselines and measurements of new parameters.For others, we provide new timingderived solutions for spin and position parameters.

New Pulsars
The majority of pulsars timed in this study have been published in preceding GBNCC survey papers (Stovall et al. 2014;Lynch et al. 2018a;Kawash et al. 2018;Aloisi et al. 2019;Agazie et al. 2021;Swiggum et al. 2023), but there are 30 sources that have been discovered within the past 2-3 yr for which no timing solution has been previously available.These sources, listed in Tables 1 and 2, include a wide range of periods and DMs.Because they were discovered following the CHIME/GBNCC data sharing agreement, many of these sources only have TOAs from CHIME observations.Along with the GBNCC discoveries, we provide timing solutions for the sources that were discovered in the GBT350 Northern Galactic Plane (NGP) survey (Hessels et al. 2008).These 12 sources are marked in Tables 1 and 2.

Glitches and Timing Noise
In this data set, sources with spin derivatives of order >2 and/or apparent discontinuities in timing residuals were examined as possible glitching pulsars.To model such a glitch, up to four parameters are used: the epoch of the glitch (given in MJD), the magnitude of the change in spin frequency (Δf, given in hertz), and the magnitude of the change in first and second frequency derivatives (  n D , in Hz s −1 and n D in Hz s −2 ).In some cases (especially those where the glitch occurs in a low-observingcadence epoch), only one of these parameters can be reliably measured.
With a high observing cadence (and a high-significance glitch with magnitude typically of order 10 −10 ν), glitches are relatively simple to resolve: A precise enough pre-glitch model can make the sudden period change quite stark.Sources that experienced glitches during the CHIME/Pulsar data are in this regime, as near-daily cadences greatly reduce our uncertainty on the glitch epoch.However, many of the sources included in this project have two or more temporally disparate observing epochs, and observations on either side of gaps are generally taken using different instruments.Distinct observing systems can introduce additional pulse phase uncertainty due to both instrumental and profile evolution effects; to mitigate this, we typically include a free "jump" parameter between different systems/frequencies.This greatly improves fits between multiple telescopes, but it does reduce sensitivity to glitches in the gaps, as the jump is highly covariant with the glitch epoch.We also acknowledge that the GBNCC survey has been shown to be more sensitive to older pulsars (McEwen et al. 2020), which are less likely to glitch than their younger counterparts (Espinoza et al. 2011;Millhouse et al. 2022).Despite this, the long baselines (and large sample of pulsars) increase the odds of observing these phenomena even in unlikely sources.
To that end, we find some pulsars in our data set seem to contain evidence for glitches and report on their characteristics here.We compare the measured glitch parameters to all published glitches in the ATNF glitch database in Figure 2 and find that all measured glitch magnitudes appear well in line with previous measurements, both in Δν and  n D . 31Fit parameters for all included glitches are given in Table 3, and residuals before and after adding glitches are shown in Figures 3-5  PSR J1923+4243 was included in two timing campaigns in 2013 (Lynch et al. 2018b) prior to observations at CHIME.This campaign resulted in a well-constrained timing model, making phase connection to the CHIME/Pulsar observations that began in 2020 trivial.However, following MJD ;59500, the solution began to drift, and the rms of the residuals began to grow.Initial attempts included adding higher-order frequency derivatives; while this did improve the fit, the parameters were not stable, and changed as more data were added.We substituted the frequency derivative terms for glitch parameters (glitch epoch and glitch magnitude), and the fit dramatically improved (reduced χ 2 = 64, with 333 degrees of freedom, DOF, initially and 2.3 for 331 DOF).We used an F-test to compare the fits; comparing the model after adding glitch parameters to the model without a glitch, the F-statistic is ≈10 −241 , indicating that the fit improvement is unlikely due to chance.
While this fit is much better than pre-glitch with spin frequency derivatives, some additional structure remained in the residuals.In an attempt to reduce this, we also tried a model that included terms for the second frequency derivative and glitch recovery (first and second frequency derivatives).Given that glitches typically correlate with significant timing noise (Lower et al. 2021), this choice is reasonable.The resulting fit no longer contained any apparent structure, the reduced χ 2 dropped to 1.2 (DOF now 329), and an F-test comparing it to the pre-glitch fit returned ≈10 −281 .Comparing instead to the simpler glitch model (containing only the epoch and magnitude of the glitch), the F-statistic is ≈10 −44 in preference for the more complex model.We plot the timing residuals before and after including the glitch model in Figure 6, and the parameters for the glitch are given in Table 3.

J0212+5222
The timing residuals for PSR J0212+5222 developed similarly to those of PSR J1923+4243: solutions after ;59200 required an increasing number of spin derivatives and still showed some signature in the TOAs.We replaced the derivatives (third order and above, maintaining the first and second) with a glitch model as described for J1923+4243, and the model did improve.However, the residuals did not become decisively flat as they had for the prior source.In particular, with the improved fit, another glitch-like event became evident in the archival GBT observations.We added a second glitch component and refit only archival data; we measured a second glitch with reasonable magnitude close to MJD 56551 (F-test comparison before/after this glitch in archival data returned 10 −20 ).We did not measure any significant recovery terms for this glitch.Including the magnitude and epoch parameters in the model when we fit all TOAs flattened the residuals.The residuals before and after adding the glitches are included in Figure 3, and the parameters for the glitches are provided in Table 3.

J2029+5459
As with the other glitching sources, fits to PSR J2029 +5459ʼs timing residuals using standard timing parameters were insufficient in removing structure from the TOAs.Unlike many other sources, the phase JUMP between archival GBT data and CHIME/Pulsar observations was not enough to solve the source, and so we focused on the higher-cadence CHIME/ Pulsar observations.A glitch was inserted near MJD 59536, and the CHIME/Pulsar TOAs flattened.However, connection to the archival data was still not possible, and it appeared as though there was a discontinuity in spin frequency between the two observing epochs that was too large to be due to quadratic spin-down.A second glitch was inserted in the gap (MJD 58000) and phase connection was made possible.Due to the lack of data, the epoch of this glitch is entirely covariant with the GBT-CHIME jump, and so we are unable to precisely measure it.The F-test comparing the pre-and post-glitch fits returned 10 −230 .Residuals for the data are shown in Figure 4, and the measured parameters are included in Table 3.

J2202+5040
Identifying the glitch progressed similarly to other sources; a model without a glitch included required several high-order polynomials to flatten TOAs, and was exchanged for a model with a glitch near MJD 59650.Similar to J2029+5459, this sufficiently flattened residuals in the CHIME data, but connection to archival GBT data required an additional glitch in the data gap.A marked difference in this source's solution is the sign of the glitches, which are both negative in Hz s −1 (indicating a spin-down glitch).Residuals are shown in Figure 5 and glitch parameters are included in Table 3.

J2351+6500
Timing residuals for this source hinted at the presence of a glitch shortly after MJD 59500, when the otherwise simple model no longer flattened the residuals.Fits with a glitch immediately resolved the problem.The measured glitch parameters (included in Table 3) are similar to those measured in PSRs J0212+5222 and J1923+4243, where the change in ν is positive (i.e., spin-up glitch) and the change in  n is negative.This is apparent in the timing residuals (Figure 7), where the epoch of the glitch precedes a downward trend that recovers quadratically.

Timing Noise and Unconfirmed Glitch Candidates
The vast majority of pulsars in this data set have been solved, meaning that their timing models correctly predict the arrival times of pulses extending into the future.There is a small subset of pulsars for which this is not true; in many of these cases, timing residuals contain additional structure that is not well fit by standard rotational/positional/binary parameters.In fact, data for PSR J1954+3852 can be fit with as many as four significant spin derivatives.Binary motion can manifest as apparent spin derivatives in timing data when the orbital period is much larger than the data span (for one such example, see Bassa et al. 2016), but no such model improves the fit.We assume that the residual structure is due to timing noise; the daily cadence at CHIME gives us unprecedented sensitivity to variations in period on short timescales.This can be quantified in a simple case where we assume the uncertainty in polynomial fits scales as N TOA s .The TOA uncertainty will depend on the telescope used and the pulsar's intrinsic brightness; given that the CHIME system's sensitivity is approximately 70% that of the 820 MHz system used for the GBT and that pulsar flux density scales with observing frequency (average spectral index of −1.4; Bates et al. 2014), we can directly compare the (analytic) uncertainties on timing noise measurement.We have assumed the pulse width does not change significantly between the observing frequencies, which is not always true.However, between 600 and 820 MHz, this change is likely negligible.GBT timing campaigns for canonical pulsars typically include roughly one observation per month with duration 30 minutes.Comparing this to weekly 10 minutes CHIME scans, the GBT measurement uncertainty is approximately twice that of CHIME; to match CHIME uncertainty, 100 minutes scans would be needed at GBT. Usage of a lower frequency at GBT can help this, but only for sources with steep spectral indices -and CHIME cadences/durations are generally greater than those used above.
One metric used to quantify timing noise is the stability parameter, defined as (Arzoumanian et al. 1994) where ν is the pulsar spin frequency, n is its second derivative, and T 8 = 10 8 s is the baseline over which the former parameters are measured.This study has produced measurements of n for 15 sources; for these and all other sources, we check for baselines in CHIME-only data that span at least 0.8 T 8 and fit for n over this period.These measurements are shown in Figure 8 with our bestfit line (with 3σ uncertainty) and the line published in Arzoumanian et al. (1994) for comparison.We note that our measurements deviate from those of Arzoumanian et al. (1994) slightly, particularly for sources with larger spin-down (i.e., younger pulsars).Models that did not previously include n (plotted with green triangles in Figure 8) agree more closely with those results.In general, though, the spread is fairly large.

Proper Motions
Measurements of pulsar proper motions over short timescales can be difficult to reliably disentangle from position, as errors in the two components will manifest as annual sinusoids.Breaking this degeneracy is possible with high cadence and long observation baselines, which we have achieved for many sources included in this study.Table 4 provides all of the sources for which we have measured proper motion.Many of these sources have been measured in previous studies; we  utilize this archival data to greatly supplement our measurements.Sources for which no previous measurements exist are marked and explained in the table footnotes.
For each measurement of proper motion, we also estimate the pulsar's transverse velocity using the DM distance of the pulsar, which we calculate using both NE2001 (Cordes & Lazio 2002) and YMW+16 (Yao et al. 2017).With a distance and a velocity measurement, we can further refine the pulsar's rate of spin-down by accounting for two effects: Galactic acceleration and transverse motion (Shklovskii 1970).These terms can be calculated and removed from the measured period derivative of a pulsar to calculate a more accurate magnetic field and spin-down luminosity.For a full description of the procedure used to estimate these terms, see Lynch et al. (2018b).Swiggum et al. (2023) also incorporates some new measurements and tools to help with such calculations.Figure 4. Glitch in PSR J2029+5459 timing residuals.The epoch of the glitch is marked with a vertical dashed line, and the horizontal axis is broken for clarity.For this source, we also fit a glitch in the gap between the GBT data and CHIME data.We are not able to constrain the epoch of this glitch, so it is omitted from the plot.
Of the 19 sources for which we measure proper motion, 13 are new measurements.The remaining six sources (J0645+5158, J1125+7819, J1641+8049, J1710+4923, J1816+4510, and J1955+6708) have published proper motions in the ATNF catalog (v1.70;Manchester et al. 2005).Some of these measurements differ from published values by as much as 79σ, or nearly 40 mas yr −1 .We include in Figure 9 a comparison between these previous measurements and those coming from this study.PSRs J0645+5158, J1125+7819, and J1816+4510 are fairly consistent, and the other three are very disparate.PSR J1641+8049 is a spider pulsar (for an overview, see Fruchter et al. 1988), which is actively ablating its companion.These systems often exhibit variable orbital parameters due to their dynamic environment (see, e.g., Polzin et al. 2019)  .Glitch in PSR J2202+5040 timing residuals.The epoch of the glitch is marked with a vertical dashed line, and the horizontal axis is broken for clarity.For this source, we also fit a glitch in the gap between the GBT data and CHIME data.We are not able to constrain the epoch of this glitch, so it is omitted from the plot.twice that used in Lynch et al. (2018b).The timing baselines for PSRs J1955+6708 and J1710+4923 have doubled, as well, which may have resulted in the significant change.

Binary Pulsars
The GBNCC survey has been an excellent tool for finding binary pulsars, with a total of 24 such discoveries.Here, we have included long-term timing solutions for 18 of its binaries, plus an additional binary pulsar from the GBT350 survey (Hessels et al. 2008).These sources span a wide range of orbital periods (0.1 P B 528 days) and eccentricities (10 −6 e 0.6), encompassing many evolutionary pathways and timing phenomena.We plot each binary system's orbital period and minimum companion mass in Figure 10.Tables 5 and 6 include the measured/derived binary parameters for all sources included in this data set.The first includes sources which use the ELL1 binary model (Lange et al. 2001), a parameterization that is appropriate for sources in nearly Figure 8. Timing stability parameter Δ 8 for GBNCC sources.For each source, we limited the baseline to ≈10 8 s in CHIME-only data.Red X markers indicate sources where n had been incorporated in the timing model previously; green triangles indicate sources where we added n to the model to determine a limit.The black dasheddotted line is the fit line from Arzoumanian et al. (1994), and the blue dashed line is a best-fit to the red points.We also include the 3σ region for the best-fit line, which shows apparent discrepancy with the prior fit; we attribute this to the number of points observed and implicit scatter in the Δ 8 parameter.circular orbits.The second includes sources that use either the more generalized binary model (BT; Blandford & Teukolsky 1976) or the relativistic parameterization (DD; Damour & Deruelle 1985, 1986) for sources with significant pulsar acceleration.Many of these binary systems have been published elsewhere, though few data sets have the long baselines/high cadences reported here.We include parameters for nine new binary solutions, and focus on several sources of particular interest below.

Post-Keplerian Parameter Measurements
Interestingly, despite the dramatic increase in timing baselines for the pulsars in this study, only one post-Keplerian term has been added to our binary models: the orbital decay term  P B for PSR J0509+3801, which we measure to be −1.37(7)× 10 −12 s s −1 .For all binaries, we calculated estimates of these parameters using the equations 8.48-8.52 in Lorimer & Kramer (2012) and determined how these the inclusion of these estimates would impact our timing solutions.Aside from the few sources with post-Keplerian measurements from previous works (see below), none of these parameters are predicted to be measured significantly.We also added them to the fit parameters as a test, and no measurements were significant above 2σ.PSRs J0509+3801 (Lynch et al. 2018b) and J1759+5036 (Agazie et al. 2021) are both in relativistic orbits with prior measurements of the advance of periastron  orbit w -the former also has a previous measurement of the Einstein delay parameter γ.We Notes.Pulsar proper motions are given in the ecliptic frame.For each source, we calculate quoted parameters using distance predictions from both NE2001 (Cordes & Lazio 2002) and YMW+16 (Yao et al. 2017) and estimate the Galactic potential using results from Guo et al. (2021).Additional parameters include the total velocity v tot , apparent spin derivatives from Galactic motion  P Gal and Shklovskii  P S corrections, the remaining intrinsic spin-down  P int , magnetic field B, characteristic age τ c , and spin-down luminosity  E. Values in parentheses are the 1σ uncertainty in the last digit as reported by PINT.Uncertainties on distances are generally ;30%.a New measurements of proper motion.b Proper-motion measurements for PSR J1710+4923 imply a negative intrinsic spin-down.This is most likely due to an error in the distance estimate and leads to an erroneous magnetic field and characteristic age.For this reason, we omit terms for this source from the final three columns.Note that this was also mentioned in Lynch et al. (2018b).
improve previous constraints on these parameters, and these improvements are included in Table 7.The improvement in precision for PSR J0509+3801ʼs binary parameters also further constrain the masses in the orbit (previously published in Lynch et al. 2018b), measured to be 1.399( 6) M e and 1.412( 6) M e , respectively, for the pulsar and its companion (shown in Figure 11).All mass measurements are consistent with those made in previous works, including the new  P B measurement (though it is not yet a strong constraint compared to those made from  w and γ).For PSR J1759+5036, our  orbit w measurement of 0.1289(4) deg yr −1 is more constraining than the previously published value (0.127(10) deg yr −1 ; Agazie et al. 2021).This further constrains the total system mass to 2.679(12) M e , which agrees within 2σ of the mass published in Agazie et al. (2021).Considering that the timing baseline for this source has approximately doubled with the inclusion of CHIME observations, such a change is not unexpected.

J2038+3447
PSR J2038+3447 (J2038) was originally discovered in the NGP survey (Hessels et al. 2008) and subsequently observed in 820 MHz observations that spanned ;226 days.The first of these observations included gridding to localize the pulsar to within an  820 MHz beam.Subsequent observations were centered on the improved position at a varying cadence (generally in groups with a few weeks between groups).Timing analysis of these data uncovered a significant period drift between observations much larger than expected from spin-down, suggesting binary motion.However, the data span was insufficient to fully disentangle its rotational, positional, and binary parameters.Further 820 and 350 MHz observations with the GBT were conducted in late 2019/early 2020, but phase connection was still not reached.
J2038 was among the first sources to be added to the CHIME/Pulsar regular observations following the CHIME/ GBNCC data sharing agreement, and has since been observed with near-daily cadence.Given the many observations and lack of binary constraint, we assumed the binary period was ?1 yr.To handle this, we opted to include measurements of frequency derivatives in our model to capture the orbital motion (for a similar case, see Kaplan et al. 2016 andBassa et al. 2016).These apparent period changes are actually due to binaryinduced Doppler shifts in the pulsar signal, and can therefore be mapped back onto the orbital phase/orbital period plane.After ;100 days of CHIME/Pulsar observations, this technique allowed us to predict a binary model that was phase-connected with the archival GBT data.Further observations broke covariance between position and spin-down and fully solved the system.
Determination of these parameters showed that J2038 is a relatively slow rotator (P ; 160 ms) in a large, nearly circular orbit (P B ; 462 days, A 1 ; 142 light seconds, e < 10 −3 ).These measurements place it in a somewhat sparse region of the orbital period/spin period plane, as we show in Figure 12.Nearby binaries are either in highly eccentric orbits with mainsequence companions or circular orbits with degenerate (mostly carbon-oxygen or helium white dwarf) companions.Given J2038ʼs low eccentricity and its minimum/median mass     6) M e and 1.412( 6) M e .We include the 3σ error regions (3σ and 1σ lines are also shown) for each parameter and highlight the maximum likelihood with a red star.The green pentagon shows the result when fitting for the masses using a model described in Taylor & Weisberg (1989), which assumes that general relativity (GR) is the correct theory of gravity (DDGR).The new measurement of ( )  P 1.37 7 10 B 12

= -
´does not strongly constrain the masses, but is consistent with GR predictions for orbital decay due to the emission of gravitational waves.The uncertainty on  P B is larger than predicted contributions from Galactic/Shklovskii terms; regardless, we include these corrections in our estimates using the DM distance of 1.6 kpc.All curves are plotted under the assumption that GR is correct. of 0.36/0.43M e (calculated from the binary mass function), the companion is most likely to be a He-core white dwarf.

Conclusion
In this work, we have analyzed timing data for 128 sources discovered in GBNCC and the GBT NGP surveys.When available, we incorporate both archival data from GBT, Arecibo, and LOFAR with recent high-cadence CHIME/ Pulsar data.We have connected TOA-generation and pulsartiming tools to produce a user-friendly, self-auditing pipeline for managing large and disjoint baselines.Using these data, we fit for spin and position parameters for all sources, including 42 newly published timing solutions.Timing residuals for all sources are provided in Appendix 5.For 19 sources (13 of which are new), we measure proper motion and infer kinematic corrections to the observed spin-down.The data set has also provided measurements of binary parameters for 19 sources and has uncovered seven glitches in four pulsars.We briefly analyzed timing noise in the population, though we leave a more involved analysis for a future work.We provide timing models, TOAs, and configuration files for reproduction of results presented here on the GBNCC GitHub page. 32he use of CHIME follow-up data for the collection of GBNCC-discovered pulsars has provided a new probe into the short-timescale variability of their timing models.Preceding timing campaigns were generally disparate in time, often with weeks to years between subsequent observations (if timing follow-up had been conducted at all).This limits the ability to measure even the "primary" model parameters (i.e., spin, spindown, and position), as phase connection can easily be lost between subsequent observations when the gap is large enough.In this way, high-cadence CHIME observations help to reach a stable solution more quickly than episodic GBT scans.For young pulsars, measurements of these parameters can be contaminated with glitches and timing noise; again, this is mitigated with high-cadence observations where glitches are visible in the timing residuals.
For those sources with long baselines, connection to CHIME data across large gaps (and between different observatories) signals that the timing model in use is robust.This requires very accurate timing positions and proper motions (which we see in Section 4.3) and a complete binary solution where applicable (Section 4.4).Hence, CHIME's timing capabilities improve upon previous timing campaigns by facilitating this phase connection. .

Figure 1 .
Figure 1.Spin period vs. period derivative for all sources timed in this work.Black and red markers indicate sources we have timed, with the latter indicating binary sources.We also plot all pulsars from the ATNF Pulsar Catalog (Manchester et al. 2005).Dashed and dashed-dotted lines indicate constant characteristic ages and magnetic fields, respectively.

Figure 2 .
Figure 2. Glitch magnitude vs. characteristic age.Fractional changes in both spin frequency (top panel) and frequency derivative (bottom panel) are shown in red, while published glitches from the ATNF glitch catalog are shown in blue.

Figure 3 .
Figure 3. Glitches in PSR J0212+5222 timing residuals.Dashed vertical lines highlight the epochs of the two measured glitches; we use a broken horizontal axis for clarity.

Figure 5
Figure5.Glitch in PSR J2202+5040 timing residuals.The epoch of the glitch is marked with a vertical dashed line, and the horizontal axis is broken for clarity.For this source, we also fit a glitch in the gap between the GBT data and CHIME data.We are not able to constrain the epoch of this glitch, so it is omitted from the plot.

Figure 6 .
Figure 6.Glitch in PSR J1923+4243 timing residuals.The dashed vertical line indicates the epoch of the glitch.For clarity, the horizontal axis is broken to omit the data gap between MJD 56700 and 59000.

Figure 7 .
Figure 7. Glitch in PSR J2351+6500 timing residuals.The epoch of the glitch is marked with a vertical dashed line, and the horizontal axis is broken for clarity.

Figure 9 .
Figure 9.Comparison of proper-motion measurements to previously published values.For each of the six sources with measurements of proper motion, we plot a hollow marker for the previous value and a filled marker for our measurement.All points have error bars, but they are small in most cases.The dramatic change in proper motion for PSR J1641+8049 is supported by optical results published by Mata Sánchez et al. (2023).

Figure 10 .
Figure 10.Orbital period vs. companion mass.Bold points indicate pulsars in our sample, and fainter points come from the ATNF catalog.Marker type indicates what binary model is used.The dashed magenta line shows the relationship between white dwarf mass and orbital period given in Tauris & Savonije (1999).

Figure 11 .
Figure 11.Mass-mass diagram for PSR J0509+3801.With three measurements of post-Keplerian parameters (advance of periastron, Einstein delay, and orbital decay), the masses of the pulsar and its companion are measured to be 1.399(6) M e and 1.412(6) M e .We include the 3σ error regions (3σ and 1σ lines are also shown) for each parameter and highlight the maximum likelihood with a red star.The green pentagon shows the result when fitting for the masses using a model described inTaylor & Weisberg (1989), which assumes that general relativity (GR) is the correct theory of gravity (DDGR).The new measurement of ( )  P 1.37 7 10

Figure 12 .
Figure 12.Comparison of PSR J2038+3447ʼs spin and orbital periods to other known binaries.We include points for carbon-oxygen companions (CO, circles), double neutron star companions (DNS, filled X's), helium/helium-CO companions (He and He/CO, triangles and squares), and main-sequence companions (MS, stars).The filled and outlined star indicates J2038.Filled plus sign markers highlight the other binaries published in this paper.For all points, colors indicate the orbital eccentricity.

Figure 13 .
Figure13.Timing residuals for sources in this study.TOAs are colored by their observatory: blue corresponds to CHIME, green to GBT, red to LOFAR, and yellow to Arecibo.

Table 1
Spin Parameters for Newly Solved Pulsars Hessels et al. (2008)heses are the 1σ uncertainty in the last digit as reported by PINT. a Pulsars discovered inHessels et al. (2008).
. As noted in Mata Sánchez et al. (2023), propermotion measurements for this source from Lynch et al. (2018b) are likely overestimated, and are at odds with measurements from optical data.This is likely what is reflected in our new measurement, as the timing baseline in this work is approximately All timing models presented here use the ELL1 binary model, which is appropriate for low-eccentricity orbits.Values in parentheses are the 1σ uncertainty in the last digit as reported by PINT. a Published rotational models, but no previously measured binary parameters.
b Not published prior to this work.

Table 6
Binary Pulsars that Utilize DD/BT Models All timing models presented here use DD/BT binary models.Values in parentheses are the 1σ uncertainty in the last digit as reported by PINT.
a Published rotational models, but no previously measured binary parameters.