Study of 72 pulsars discovered in the PALFA survey: Timing analysis, glitch activity, emission variability, and a pulsar in an eccentric binary

Survey data were collected with the ALFA receivers, which consisted of seven beams with half power widths of 3.6^′, and data were recorded with the Mock spectrometers⁴⁴4www.naic.edu/~astro/mock.shtml. This backend processed two overlapping 172-MHz bands which, once combined and the two polarizations summed, provided 323 MHz of bandwidth centered at 1375.5 MHz and 960 frequency channels sampled every 64 $\mu$s. PALFA targeted two regions of the Galactic plane ($|b|\,<\,5^{\circ}$): the “inner” Galaxy ($32^{\circ}\,\lesssim\,l\,\lesssim\,77^{\circ}$) and the “outer” Galaxy ($168^{\circ}\,\lesssim\,l\,\lesssim\,214^{\circ}$). Integration times were 268 s and 180 s for inner and outer Galaxy observations, respectively. Additional details on survey observations for PALFA, including the strategy for surveying the inner and outer Galaxy regions, are provided in Section 7.

Survey data have been processed and searched by three independent pipelines. The first is a reduced-resolution “Quicklook” pipeline performed in near real-time on site that enables rapid discovery of bright, nearby pulsars (Stovall, 2013). The second is a full-resolution PRESTO⁵⁵5https://github.com/scottransom/presto-based pipeline (Ransom, 2011). The latter processes data on the Béluga supercluster, a Compute Canada/Calcul Québec facility hosted at the École de technologie supérieure in Montréal. It searches for dispersed periodic signals in the Fourier domain as well as in the time domain with a Fast-Folding Algorithm (Parent et al., 2018). Sporadic pulses such as those emitted by RRATs and FRBs are searched for with a single-pulse pipeline (Patel et al., 2018). Lastly, data are searched for pulsars, particularly those in compact orbits, using the Einstein@Home⁶⁶6https://einsteinathome.org/ pipeline (E@H) described in Allen et al. (2013).

Approximately 40% of the 72 objects we study in this work were first found by either the Quicklook or E@H pipelines. The remaining sources were only found by the full-resolution pipeline. The pipeline search algorithm that identified each source is listed in Table 1 along with the date of discovery. The pulsars’ time-integrated pulse profiles are shown in Figure 1.

Timing observations were largely conducted at AO and, for the 33 strongest sources, with the 76-m Lovell Telescope at Jodrell Bank Observatory (JBO) in Macclesfield, UK. Until 2014, some limited amount of timing data was collected for six of our pulsars with the 100-m Green Bank Telescope (GBT) at the Green Bank Observatory (GBO) in West Virginia, US. In 2020 May, we began follow-up observations on five other sources with the Canadian Hydrogen Intensity Mapping Experiment (CHIME) telescope⁷⁷7https://chime-experiment.ca/en (CHIME/Pulsar Collaboration et al., 2021), located at the Dominion Radio Astrophysical Observatory (DRAO) in British Columbia, Canada. More information on observations for individual sources such as observations sites, cadence, and time span of observations can be inferred from plots introduced in Section 3.1.

At AO, timing data were collected with the L-Wide receiver and recorded with the PUPPI⁸⁸8https://www.naic.edu/puppi-observing/ (Puerto Rican Ultimate Pulsar Processing Instrument) backend at a central frequency of 1380 MHz with a nominal 800-MHz bandwidth split into 4096 channels sampled every 40.96 $\mu$s. Following RFI excision, the usable bandwidth was typically $\sim$600 MHz and integration lengths ranged from 300 to 1500 s, depending upon the pulsar brightness. Initially, data were recorded in incoherent search mode but as a pulsar ephemeris improved, we switched to fold-mode observations where data were coherently dedispersed at the pulsar DM and folded into 10-s subintegrations in real-time modulo the instantaneous pulsar period, producing time-integrated pulse profiles for all frequency channels.

Timing data obtained with the Lovell Telescope were processed using a Digital Filterbank (DFB) which Nyquist samples a 512 MHz band at 8-bit resolution and channelised it into 1024 channels using a polyphase filter. After RFI cleaning, the resultant band spans the range 1350 to 1700 MHz. The data are incoherently dedispersed and folded into 10 s long sub-integrations with 1024 pulse phase bins. Observation durations ranged from 30 to 60 minutes depending on the source flux density. More information can be found in Lyne et al. (2017b).

Observations with GBT were recorded with the Green Bank Ultimate Pulsar Processing Instrument (GUPPI, DuPlain et al. 2008) backend, which processes 800 MHz of bandwidth centered at 1.5 GHz. Coherent-dedispersion mode data were recorded into 512 frequency channels and folded into subintegrations every 10 s.

As for CHIME observations, the stationary instrument operates in the 400-800 MHz frequency range and beamformed data were collected with the CHIME/Pulsar backend. The latter discretizes the observing band into 1024 baseband channels which are then coherently dedispersed. The “Digital Signal Processing for Pulsars” (DSPSR) suite⁹⁹9http://dspsr.sourceforge.net/ is used to fold the data into 10 or 30 s subintegrations, depending on the source’s DM. The duration of any given observation at CHIME is limited by the transit time as set by the declination, $\delta$, of each pulsar. On average, the five objects we followed up with CHIME (all with $16^{\circ}\,\lesssim\,\delta\,\lesssim\,27^{\circ}$) were observed for $\sim$ 15 minutes each per day. More information on the CHIME/Pulsar system can be found in CHIME/Pulsar Collaboration et al. (2021).

In order to calculate the times of arrival (TOAs) of radio pulses in our data sets, we first excised RFI. Where only AO search-mode data were available, the cleaned raw data were folded at the topocentric period and dispersion measure (DM) that yielded the strongest detection with PRESTO’s prepfold tool. We then produced a standard profile template by fitting one or more Gaussian components to the integrated profile that was detected with the highest signal-to-noise ratio (S/N). TOAs were then extracted by cross-correlating the folded data with the standard template using the get_TOAs.py program from PRESTO, which fits for a linear phase gradient in the Fourier domain to determine the shifts between the profiles and the standard template (Taylor, 1992).

When analyzing fold-mode data, we constructed improved, high-S/N standard templates by summing in phase the pulse profiles from multiple observations with the psradd tool from the PSRCHIVE¹⁰¹⁰10http://psrchive.sourceforge.net/ software package. As profiles were being combined, weights were applied based on the signal-to-noise ratio of the pulsar signal in each data set. We then set profile baselines to zero to create “noise-free” templates before smoothing the final profiles with PSRCHIVE’s psrsmooth tool. TOAs were then extracted using the Fourier phase gradient approach with pat, also from the PSRCHIVE package. Given the difference in observing frequencies and the time delays introduced by backend systems, a separate standard profile template was generated to analyse data collected at each site.

The timing analysis was carried out with the TEMPO  software package¹¹¹¹11http://tempo.sourceforge.net/, which implements a $\chi^{2}$ minimization technique to compute the best-fit parameters of a timing model. The JPL DE436 planetary ephemeris¹²¹²12https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/spk/de436s.bsp.lbl and the UTC(NIST) time standard¹³¹³13https://www.nist.gov/ were used. The basic timing model for each source was parameterized by the pulsar period and first period derivative, $P$ and $\dot{P}$, right ascension, declination, and DM. To fit for DM, we extracted TOAs at different frequencies for a number of epochs using between two and eight frequency subbands per epoch, depending on the signal strength. Arbitrary time offsets between TOAs collected at different observatories were also allowed (where applicable) as timing parameters. Orbital parameters were also included in the timing model of the binary pulsar PSR J1954+2529, which we describe later in Section 4.7. In some cases, we additionally fit for higher-order frequency derivatives to reduce scatter due to the timing noise, an effect that arises from spin irregularities that are intrinsic to the pulsar (Manchester & Taylor, 1977). For the five pulsars that displayed glitch activity (Section 4.5), glitch parameters were also computed with TEMPO.

Best-fit timing parameters are provided in Table 3 and corresponding post-fit residuals, with RMS values that range from 86 $\mu$s to 5.4 ms, are shown in Figure 3. Other measured and inferred pulsar properties are listed in Table 4.

Coherent solutions were obtained for all but three pulsars: PSRs J1855$+$0626, J1858$+$0239 and J1928$+$1725. The first pulsar, PSR J1855$+$0626, is an intermittent pulsar with too few detections to enable phase connection. We discuss its intermittency and emission properties in Section 4.3.1. PSR J1858$+$0239 is a pulsar whose unstable average profiles appeared to show evidence of mode-changing behavior, which we discuss in Section 4.1. For that source, timing observations were conducted solely at AO and we monitored the pulsar for $\sim$1.5 years. We were able to phase connect a subset of TOAs from $\sim$ 10 consecutive epochs within a dense-observing timing campaign (spanning approximately one month) during which time the average pulse profile remained fairly stable and could thus be analysed with one common standard profile. Within that TOA set, we detected a significant spin-down rate (see $\dot{P}$ reported in Table 3). However, when we attempted to connect time gaps between other consecutive pairs of closely spaced observations in which the pulsar displayed similar average profiles, we could not produce a consistent solution. Although we clearly observed a change in the observed spin frequency, an accurate $\dot{P}$ measurement requires that the timing coherence extends over at least one year in order to break the covariance between spin and astrometric parameters. Furthermore, we suspect that magnetospheric activity may coexist with torque variability (Kramer et al., 2006), thus the putative $\dot{P}$ reported in Table 3 may not be representative of the pulsar’s normal spin-down (see Section 4.1.2). Dedicated, long-integration and regular observations would be required to obtain an accurate description of this pulsar’s rotation. Finally, the last unsolved source, PSR J1928$+$1725, is a RRAT that could only be timed through its single pulses. All timing observations were conducted at AO. A number of bright, narrow single pulses were detected but at a very irregular rate, with clusters of pulses being emitted within a few 15-minutes integrations and no detectable emission in several data sets, leaving large time gaps between detections. We extracted topocentric TOAs from each single pulse and attempted to achieve phase connection but were unsuccessful. Similarly to PSR J1858$+$0239, we suspect that PSR J1928$+$1725 experiences large torque variability. We discuss in more detail the properties of this RRAT in Section 4.4 and the evidence for the possible changes in spin-down rate.

We estimated the average pulsed flux densities at 1400 MHz, $S_{1400}$, by calibrating the ALFA discovery data using the radiometer equation (Dewey et al., 1985). The ALFA receiver provides a nominal bandwidth of 323 MHz, however effective bandwidths typically range between 260 MHz and 300 MHz following RFI excision. The latter, narrower bandwidths were used in the calibration procedure. Sky temperatures at the pulsar positions were estimated by extrapolating the 408-MHz all-sky map from Remazeilles et al. (2015) to 1400 MHz, assuming a spectral index¹⁴¹⁴14Defined as $S_{\nu}\propto\nu^{\alpha}$, where $S_{\nu}$ is the flux density at frequency $\nu$ and $\alpha$ is the spectral index. of –2.7 (Remazeilles et al., 2015). System temperatures of the ALFA receiver typically varied between 28 K and 32 K. In estimating the flux density, we assumed a system temperature of 30 K. We used a value of 9 K Jy^-1 for the gain of the central ALFA beam, and scaled the gain of the outer six beams to be 79% that of the central beam (Cordes et al., 2006). The data analysis procedure discards any rotation-independent radio flux from the pulsar.

Our average pulsed flux density measurements are reported in Table 4. Considerable systematic uncertainties affect our measurements, arising notably from fluctuations in the system temperature, reductions in the effective gain due to variations in the receiver response and positional offsets between the true position of the pulsars (i.e., timing positions) and the beam center positions. These uncertainties are included in the error estimates reported in Table 4. The same approach was used to estimate the pulsed flux density of the 23 additional pulsar discoveries reported in Table 2 of Section 2.1.1.

We note here that the pulse profiles shown in Figure 1, most of which were created from data collected with the 800-MHz L-Wide receiver at Arecibo, are not the profiles that were used to estimate the average pulsed flux densities. Due to the difference in the spectral response of the ALFA and L-Wide receivers and because pulsars generally have power-law spectra and their pulse profiles evolve with observing frequency, the profiles in Figure 1 cannot be directly normalized to the flux densities reported in Table 4.

Integrated pulse profiles were examined to identify asymmetric broadening that would be indicative of frequency-dependent scattering by the turbulent interstellar medium. We attempted to quantify ISM scattering in each source by fitting a temporal pulse-broadening function to the observed pulse profile in different frequency bands. In our model, we assume that the radio waves are scattered isotropically by a single thin scattering screen. The pulse broadening function associated with this assumption takes the form $\tau_{sc}^{-1}e^{-1/\tau_{sc}}$, where $\tau_{sc}$ is the characteristic scattering time which scales with observing frequency as $\tau_{sc}(\nu)\,\propto\,\nu^{-\alpha_{sc}}$.

In the latter expression, $\alpha_{sc}$ is the scattering spectral index. We simultaneously fit the subbanded pulse profiles as a single-component Gaussian convolved with the scattering broadening function, and compared the goodness of the fit to that of a Gaussian mixture that includes up to three components. For several of our low-S/N pulsars, we did not fit subbanded profiles but instead profiles integrated over the entire observing band.

Only four pulsars display pulse shapes that are best described by the scatter-broadening model. We show in Figure 4 the fits to the subbanded pulse profiles, and best-fit scattering spectral indices $\alpha_{sc}$ and broadening timescales $\tau_{sc}$ at 1 GHz are provided in Table 5 along with the scattering timescales predicted by the NE2001 (Cordes & Lazio, 2002) Galactic electron density model. Another model that has been widely used is the YMW16 (Yao et al., 2017) model, but unlike the NE2001 model, it does not use scattering as a modeling parameter. Instead, it estimates $\tau_{sc}$ at 1 GHz for a given DM value based on the empirical scaling between scattering timescale and DM obtained by Krishnakumar et al. (2015). We include those estimates in Table 5 as well. We note significant discrepancies (by up to two orders in magnitude) between the NE2001 predictions and our measurements. These inconsistencies could be attributed to unmodeled foreground structures such as H II regions. New distance and scattering measurements along various lines of sight are thus valuable to construct more complete models of the Galactic electron density in the future.

PSRs J1911$+$0925 and J1924$+$1713 have scattering indices consistent with both an isotropic scattering mechanism ($\alpha_{sc}$ = 4, Cronyn 1970) and Kolmogorov turbulence in a cold plasma ($\alpha_{sc}$ = 4.4, Lee & Jokipii 1976; Rickett 1977), but PSRs J1907$+$0833 ($\alpha_{sc}$ = 3.0) and J1913$+$11025 ($\alpha_{sc}$ = 1.6) have flatter scattering spectra than the aforementioned theoretical models. Anomalous scattering arising from structures in the ISM (Cordes & Lazio, 2001; Rickett et al., 2009) or anisotropic scattering mechanisms (Stinebring et al., 2001; Tuntsov et al., 2013) could explain the lower $\alpha_{sc}$ values. PSRs J1907$+$0833 and J1913$+$11025 also have large DM-estimated distances, both exceeding 8 kpc. As such, their low $\alpha_{sc}$ could also be explained by the presence of multiple scattering screens along their long lines of sight. Dedicated observations would be required to probe the nature of the scattering observed in these pulsars – our flux- and band-limited data sets do not allow us to distinguish between models.

Figure 5 shows the distribution of $\tau_{sc}$ at 1 GHz as a function of DM values. For comparison, we additionally plot measurements reported in the ATNF catalog (version 1.65) for other known pulsars. The best-fit solution to the empirical scaling between $\tau_{sc}$ and DM from Krishnakumar et al. (2015) at 327 MHz, scaled to a frequency of 1 GHz using the average value of our best-fit $\alpha_{sc}$ measurements, is also shown in Figure 5 (solid line). We see that our sources, which have high DMs and large $\tau_{sc}$, are consistent with the scaling relation from Krishnakumar et al. (2015) for an average $\alpha_{sc}$ of 3.35.

We also note that two low-DM sources, PSRs J1926$+$1614 (DM = 24.0 pc cm^-3) and J1939$+$2609 (DM = 47.3 pc cm^-3) displayed obvious signs of scintillation features in their spectra. From epoch to epoch, we observed fluctuations in the flux densities of both pulsars ranging from roughly 20 to 200 $\mu$Jy. Measurement of the scintillation bandwidth, and hence the scattering delay, for these and other low-DM PALFA sources could place complementary constraints to our pulsar broadening measurements for high-DM pulsars. This analysis is however beyond the scope of this paper.

Mode changing (or switching) is a type of discontinuous transition in the radio emission where the average pulse profile abruptly switches between two or more quasi-stable states (Backer, 1970a). It is a broadband phenomenon (Bartel et al., 1982) that occurs on variable timescales. Changes in the radio beam emission pattern, and hence in the observed pulse profile, are believed to be the result of a global redistribution in the magnetosphere currents and/or magnetic fields. Mode changing has been recognized in roughly two dozens pulsars thus far (see e.g., Wang et al. 2007; Lyne et al. 2010; Ng et al. 2020). Below we describe the emission of three pulsars that displayed mode-changing behavior.

First reported by Backer (1970b), the nulling phenomenon is a sudden cessation of detectable pulsed emission that lasts for one or more pulsar cycles. RRATs and intermittent pulsars (discussed in the next sections) are extreme manifestations of pulsar nulling. The triggering mechanism responsible for nulling remains largely unknown, but it is believed to be intimately related with mode changing (e.g., van Leeuwen et al. 2002; Redman et al. 2005; Wang et al. 2007). In general terms, most interpretations of the absence of emission in nulling pulsars invoke processes occurring in the magnetosphere, for example a loss of plasma conditions required for coherent emission (Filippenko & Radhakrishnan, 1982), or intense time-varying pulse modulations where the radio flux density drops below the detection threshold (e.g., Esamdin et al. 2005).

Nulling behavior has been reported in approximately 10% of the radio pulsar population, but the proportion of pulsars experiencing nulling could be much larger since most known pulsars are not monitored regularly for the purpose of identifying and studying nulling behavior. For instance, through regular monitoring with CHIME/Pulsar, Ng et al. (2020) recently reported on the first detections of nulls from bright pulsars discovered decades ago, suggesting that nulling is a phenomenon much more common than previously thought. On the other hand, misinterpreting undetectable emission in data from a given (sensitivity-limited) instrument as a null, as opposed to an extrinsic reduction in the radio flux density (e.g., due to a reduction in intrinsic luminosity or due to scintillation), results in an overestimation of the size of this sub-population.

Among the pulsars presented in this work, we identify 12 pulsars that displayed discernible nulls. We show examples of the emission intensity as a function of time in Arecibo observations for each pulsar in Figure 10. Because of the low flux densities of these pulsars, data had to be folded into subintegrations ranging from 2 to $\sim$20 pulsar cycles in length. Only the RRAT, PSR J1928+1725 (whose properties are described separately in Section 4.4), has individual pulses with signal-to-noise ratios sufficiently high that averaging over multiple rotation is not required.

Dedicated observations with long integrations are preferable for characterizing nulling properties such as nulling fractions and average duration of nulling episodes. This is especially important for pulsars showing long nulling episodes. Such observations were not carried out for this work; all the pulsars here were observed through the same follow-up program at AO for timing purposes, with scans duration ranging from 5 to 15 minutes. Thus, we do not attempt to determine nulling parameters, except for RRAT J1928$+$1725 which we discuss in Section 4.4.

While the distribution of nulling fractions in intermittent pulsars can be similar to that in nulling pulsars (e.g., Gajjar et al. 2012), the cessation of pulsed emission in intermittent pulsars is seen on timescales lasting from days to years (e.g., Kramer et al. 2006; Lyne et al. 2017a), orders of magnitude longer than nulling pulsars. Below we discuss the properties of two intermittent pulsars identified in our sample: PSRs J1855$+$0626 and J1925$+$2513.

PSR J1928$+$1725 ($P\,=\,289.8\,$ms) is a RRAT that emits bright, heavily clustered pulses followed by long (hundreds of rotations) nulls. Other than the single pulses themselves, no underlying signal is visible when folding the data. We analyzed the distribution of wait times between consecutive pulses, and find that this distribution is inconsistent with a Poisson process (see Figure 12). The longest active phase we observe lasts 29 consecutive pulsar cycles, while the longest inactive phase lasts 1662 cycles ($\sim$482 s). Active phases are also rare – we detect neither single-pulse or periodic emission in 3/17 of our 15-min observations of PSR J1928$+$1725. The average detection rate is 78 pulses per hour.

Despite having detected numerous bright pulses, the long time gaps between observations during which pulses are detected makes it difficult to maintain timing coherence over a few epochs, and thus we have been unable to solve this RRAT. PSR J1928$+$1725 could be similar to the young and energetic ($\tau_{c}\,=\,865$ kyr and $\dot{E}\,=\,4.6\times 10^{34}$ erg s^-1) 125-ms PSR J1554$-$5209 (Keane et al., 2010), one of the very few RRATs having a spin period $<\,500\,$ms for which $\dot{P}$ has been measured¹⁵¹⁵15Based on the ATNF catalog, Manchester et al. (2005). Keane et al. (2011) conducted a timing analysis on PSR J1554$-$5209 along with a dozen more RRATs with longer periods and lower $\dot{E}$, and noted that PSR J1554$-$5209 displays a far larger scatter in its timing residuals than the other, lower-$\dot{E}$ RRATs. Significant jitter in pulse phase similar to PSR J1554$-$5209 could explain why we have been unable to derive a coherent timing ephemeris for PSR J1928$+$1725. During our phase connection attempt, it also appeared that the pulsar may be suffering from significant glitch activity, although we cannot confirm this. Of the 34 known RRATs with reported spin-down rate measurements (Manchester et al., 2005), PSR J1819$-$1458 (McLaughlin et al., 2006) is the only one for which glitch activity has been reported (Bhattacharyya et al., 2018). While the former has a much longer rotation period ($P$ = 4.263 s) than PSR J1928$+$1725, it has a relatively small characteristic age of 120 kyr. If PSR J1928$+$1725 is indeed young and energetic, the presence of glitches could very well explain why we have been unable to solve it.

Follow-up observations of PSR J1928$+$1725 were carried out with CHIME/Pulsar from 2020 October to 2021 April, but no pulses were detected. Establishing phase coherence would require regular and long-integration observations with an instrument whose sensitivity is similar to that of Arecibo.

Rotational instabilities in pulsars are generally explained by either timing noise or glitches. Whereas timing noise arises from random spin fluctuations over long timescales (e.g., one of our pulsars, PSR J0608$+$1635, displays intense timing noise; see residuals in Figure 3), a glitch is a discrete transition in the pulsar rotational state marked by an abrupt increase in spin frequency, typically with frequency jump magnitudes between $10^{-6}\,\nu$ and $10^{-10}\,\nu$, and often accompanied by a decrease in frequency derivative $\dot{\nu}$ (Espinoza et al., 2011). Glitches are believed to be caused by erratic transfers of angular momentum from the superfluid inside the star to the more slowly rotating (and cooling) crust (Anderson & Itoh, 1975; Ruderman et al., 1998). Hence their study represents a unique opportunity to gain insight into the internal structure of neutron stars.

Eight glitches of moderate and small magnitude were detected in five of the 72 pulsars presented here during the timing program. Figure 13 illustrates the glitch signatures in timing residuals when they are not included in the timing model, signatures that are characterized by the sudden onset of a steady decrease towards negative residuals. We are able to maintain phase coherence over the time gap between observations around the glitch epochs when using the best-fit glitch parameters listed in Table 6 in our timing models. Residuals corresponding to the ephemerides that contain the best-fit glitch parameters are the ones shown in Figure 3. As a result of sparse observation cadences and/or high levels of timing noise and/or large TOA uncertainties, we are unable to detect the changes in the spin-down rate ($\Delta\dot{\nu}$) during recovery for three glitch events.

The three youngest pulsars presented in this work, PSRs J1910$+$1026 ($\tau_{c}$ = 33 kyr), J1931$+$1817 ($\tau_{c}$ = 35 kyr), and J1915$+$1150 ($\tau_{c}$ = 116 kyr), have all exhibited glitch activity. This is not surprising, as young pulsars are known to display higher glitch activity (e.g., McKenna & Lyne 1990; Lyne et al. 2000; Janssen & Stappers 2006; Espinoza et al. 2011).

The largest glitch we observed occurred in PSR J1910$+$1026, the youngest pulsar in our set, with a fractional step in spin frequency of $78.4\pm 0.8\times 10^{-9}\,\nu$. Due to its low flux density (58 $\mu$Jy at 1.4 GHz), timing data were solely collected at AO and the pulsar was observed roughly twice a month for $\sim$ 1.5 years. The sparse observations combined with the relatively large timing residuals ($\sim\,1\,$ms) makes the identification of smaller glitches difficult. For the aforementioned reasons, the glitch activity of PSR J1910$+$1026 is not well constrained and could be much higher than what our dataset reveals.

On the other hand, PSRs J1931$+$1817, J1939$+$2609 and J1954$+$2529 were bright enough to be followed up at JBO and CHIME/Pulsar, and benefited from more regular timing observations over a longer time span. Over the course of our 5-yr follow-up campaign of PSR J1931$+$1817, we detected four glitches ranging from 3.3$\times 10^{-10}$ to $3.7\times 10^{-8}\,\nu$ in size. A small glitch ($\Delta\nu/\nu=6.0\times 10^{-10}$) was observed in PSR J1954$+$2529, a relatively old ($\tau_{c}\,=\,$11.7 Myr) non-recycled pulsar in a binary system which we discuss in Section 4.7.

Even more surprising is that PSR J1939$+$2609, an old pulsar with characteristic age $\tau_{c}\,\sim\,1.4\,$Gyr, suffered a glitch, albeit of small magnitude 3.3$\times 10^{-10}\,\nu$. Only two other glitching pulsars with ages $>\,100\,$Myr are known in the Galactic field: PSRs B1913$+$16 with $P\,=\,59\,$ms and $\tau_{c}\,=\,109\,$Myr (Hulse & Taylor, 1975) and J0613$-$0200, a 3.1-ms MSP with $\tau_{c}\,=\,5\,$Gyr (Lorimer et al., 1995). Similarly to PSR J1939$+$2609, only one glitch of small magnitude has been detected in each system, with $\Delta\nu/\nu=3.7\times 10^{-11}$ for PSR B1913$+$16 (Weisberg et al., 2010; Weisberg & Huang, 2016) and 2.5$\times 10^{-11}$ for PSR J0613$-$0200 (McKee et al., 2016).

We note however that for PSRs J1939$+$2609, J1954$+$2529 and J1931$+$1817 (first glitch; MJD 57645), we cannot positively rule out the possibility that the features in the timing residuals are caused by timing noise rather than glitch activity. While the latter explanation minimizes the post-fit RMS residuals in these systems and is thus currently preferred over the former, future modelling of their long-term timing behavior will be needed to be certain that those are real glitches and not part of large timing noise features.

Among the 69 pulsars with coherent timing solutions, only two are young ($\tau_{c}<10^{5}$ yr) objects: PSRs J1910$+$1026 ($P\,=\,531.5\,$ms) and J1931$+$1817 ($P\,=\,234.1\,$ms) with characteristic ages of 33 and 35 kyr, respectively. Pulsars of such young age are often found in supernova remnants (SNR). Despite the lack of proper motion measurements – important to confirm any potential pulsar/remnant association – we searched for coincident SNR by cross-matching the pulsar positions against Green’s SNR Catalog (Green, 2019). Accounting for potential angular offsets that could arise from large post-supernova tangential velocities, we search over a conservative region of 0.5^∘ radius around each pulsar. No remnant coincident with PSR J1910$+$1026 was identified, which is not surprising given its very large inferred distance ($\gtrsim$ 15 kpc). One catalogued object, SNR G53.41+0.03 (Anderson et al., 2017), is found 28^′ away from the timing position of PSR J1931$+$1817. Driessen et al. (2018) carried out a deep radio search for a pulsar associated with G53.41+0.03, covering a 10^′ region around the center of the SNR (the estimated size of the remnant), but did not detect pulsations. They also investigated and characterized G53.41$+$0.03 using multi-wavelength data, and estimated its distance at roughly 7.5 kpc and its age at $\sim$ 1000 to 8000 yrs. Given the 28^′ angular offset and adopting a distance of 7.5 kpc and the most conservative age of 8000 yrs, an unreasonably high spatial velocity ( $>\,12\times 10^{3}$ km s^-1) would be required for PSR J1931$+$1817 to have traveled from its birth location near the center of the remnant to its current position. We therefore rule out SNR G53.41$+$0.03 as being associated with PSR J1931$+$1817.

Due to their large spin-down power $\dot{E}$, young pulsars sometimes have high-energy counterparts. However, we did not find any gamma-ray point sources or pulsations in Fermi LAT data for these pulsars (analysis described in Section 5). This is also not surprising given the large distances of the pulsars and the corresponding “heuristic” energy fluxes $\sqrt{\dot{E}}/d^{2}$ of $\sim 3\times 10^{15}$ and $\sim 5\times 10^{15}$ (erg s^-1)^1/2 kpc² for PSRs J1910$+$1026 and J1931$+$1817, respectively, below the approximate LAT threshold of $10^{16}$ (erg s^-1)^1/2 kpc² (Abdo et al., 2013; Smith et al., 2019).

Of all the new pulsars presented here, only one, PSR J1954+2529, is a member of a binary system. It has an orbital period of 82.7 days and, interestingly, a system eccentricity of 0.11. The system has been monitored with Arecibo and CHIME, but mostly with the Lovell telescope. PSR J1954+2529 is also among our five pulsars that displayed a glitch. The Keplerian parameters of the pulsar’s orbit and some derived quantities are presented in Table 7. Included in our timing model is the relativistic advance of the angle of periastron, $\dot{\omega}$, which in GR is a function of the total mass of the system, $M_{\rm tot}$. Unfortunately our timing data do not yet provide a significant measurement (see upper limit in Table 7), and the 3-$\sigma$ upper limit on$\dot{\omega}$ dot doesn’t yield a meaningful constraint on $M_{\rm tot}$.

The pulsar does not appear to be recycled: with $P\,=\,0.931\,\rm s$, a characteristic age of about 12 Myr and a magnetic field at the surface of about $10^{12}\,$G, it is near the center of the cloud of “young” pulsars in the $P$-$\dot{P}$ diagram (see Figure 6). Such middle-aged radio pulsars in binary systems are rare but interesting from the point of view of stellar evolution.

In what follows, we discuss three possibilities for the formation and nature of this system. At the moment, we have no data to decide in favor of any hypothesis, but the last hypothesis is currently preferred based on theoretical arguments from stellar evolution theory.