The Tail of PSR J0002+6216 and the Supernova Remnant CTB 1

TL;DR

This study combines radio and gamma-ray observations to measure the proper motion of pulsar PSR J0002+6216, linking it to the supernova remnant CTB 1 and providing insights into pulsar birth and supernova asymmetries.

Contribution

It presents the first combined radio and gamma-ray analysis of PSR J0002+6216, measuring its proper motion and associating it with CTB 1, revealing high pulsar velocity and implications for supernova dynamics.

Findings

Pulsar is at the tip of a bow-shock nebula pointing to CTB 1

Proper motion measured at 115±33 mas/yr, indicating a velocity of 1100 km/s

Pulsar likely originated from the supernova that created CTB 1

Abstract

We have carried out VLA imaging and a Fermi timing analysis of the 115 ms gamma-ray and radio pulsar PSR J0002+6216. We found that the pulsar lies at the apex of a narrowly collimated cometary-like 7 arcmin tail of non-thermal radio emission which we identify as a bow-shock pulsar wind nebula. The tail of the nebula points back toward the geometric center of the supernova remnant CTB 1 (G116.9+0.2) 28 arcmin away, at a position angle $\theta_\mu=113^\circ$. We measure a proper motion with 2.9$\sigma$ significance from a Fermi timing analysis giving $\mu$=115$\pm$33 mas yr$^{-1}$ and $\theta_\mu=121^\circ\pm{13}^\circ$, corresponding to a large transverse pulsar velocity of 1100 km s$^{-1}$ at a distance of 2 kpc. This proper motion is of the right magnitude and direction to support the claim that PSR J0002+6216 was born from the same supernova that produced CTB 1. We explore the implications for pulsar birth periods, asymmetric supernova explosions, and mechanisms for pulsar natal kick velocities.

Tables

Table 1. Table 1: Median and 1- σ 𝜎 \sigma confidence intervals for Fermi -LAT Timing Model

Parameter	Value
Right Ascension ($\alpha$, J2000)	00h02m58$\stackrel{s}{.}$14(1)
Declination ($\delta$, J2000)	+62°16′09$\stackrel{\prime \prime }{.}$52(8)
Proper Motion in R.A. (${\mu }_{\alpha }\cos \delta$, mas yr-1)	97$\pm$ 32
Proper Motion in Decl. (${\mu }_{\delta }$, mas yr-1)	-57$\pm$27
Proper Motion Magnitude (mas yr-1)	115 $\pm$ 33
Proper Motion Position Angle (degrees from North)	121 $\pm$ 13
Epoch of position (MJD)	56500.0
Timescale	TDB
Solar System Ephemeris	DE421

Full text

pulsars: individual (PSR  J0002+6216), supernova remnants, ISM: individual objects(CTB 1), proper motions

††facilities: VLA, Fermi††software: CASA

1 Introduction

Deviations from spherical symmetry appear to be an essential ingredient in successful core collapse supernovae (SN) explosions, and the observational signatures of this asymmetry include the morphology and kinematics of supernova remnant (SNR) ejecta and the natal kick velocities of pulsars (e.g., Holland-Ashford et al., 2017). The high-velocity outliers in the pulsar velocity ($V_{PSR}$) distribution are of special interest because they provide stringent tests of neutron star kick mechanisms (Janka, 2017). Characterizing this tail in the kick distribution is a timely topic since, for example, it can affect the fraction of binary neutron star systems that remain bound and their distribution within the host galaxy (Berger, 2014; Vigna-Gómez et al., 2018). There are currently only a small number of pulsars that have well enough measured proper motions and distances to robustly claim $V_{PSR}$ in excess of 1000 km s*-1* (Hobbs et al., 2005; Shternin et al., 2017). Similar high velocities have been inferred from the offset of pulsars from the center of SNRs (Frail & Kulkarni, 1991; Frail et al., 1994) but the burden of proof is high, and in at least two cases follow up scintillation or proper motion measurements do not confirm high velocities (Nicastro et al., 1996; Zeiger et al., 2008).

Recently it was proposed that PSR J0002+6216 and the SNR CTB 1 may be physically associated, based on their angular proximity on the sky and roughly similar distances (Zyuzin et al., 2018). The 115 ms $\gamma$-ray and radio pulsar PSR J0002+6216 is one of the newest additions to the Fermi sample (Clark et al., 2017; Wu et al., 2018). It has a large spin-down energy $\dot{{E}}$=1.53 $\times 10^{35}$ erg s*-1* and a dipole magnetic field $B=0.8\times 10^{12}$ G. The SNR CTB 1 (G116.9+0.2) has been extensively studied at all wavelengths. In optical H$\alpha$ and non-thermal radio it shows a well-defined circular shell of radius $\theta_{s}$=17.8′ (Landecker et al., 1982), while the X-rays are centrally concentrated, making it one of a small number of mixed morphology remnants (Craig et al., 1997; Lazendic & Slane, 2006; Pannuti et al., 2010; Katsuragawa et al., 2018). Given the angular offset (28′) and estimates for the age (104 yr) and distance (2 kpc; see §4), a physical SNR-PSR association would require $V_{PSR}>$1000 km s*-1*. Accordingly, we have carried out new VLA imaging observations and we have re-analyzed Fermi timing data. Our data support the hypothesis that PSR J0002+6216 is a high-velocity pulsar that originated from the same SN that produced the SNR CTB 1.

2 Imaging Observations

We observed a field toward PSR J0002+6216 with the Karl G. Jansky Very Large Array (VLA) as part of project 17B-384 using Director’s Discretionary Time (Perley et al., 2009). Observations were carried out on 2017 August 19 with the VLA in the C configuration. We observed in the 1-2 GHz frequency range (L band) using the standard Wideband Interferometric Digital Architecture (WIDAR) correlator setup for continuum observing with 16 spectral windows, 64x1-MHz wide channels each to get 1 GHz of total bandwidth centered on 1.52 GHz. Data were saved in 5 s integrations and the total time on source was 63 min. The radio source J2350+6440 was used as a phase calibrator, while 3C 48 (J0137+331) was used as both the bandpass and flux density calibrator.

The correlated visibilities were calibrated using CASA 5.4.0-32 together with the automated VLA pipeline version 41722 (McMullin et al., 2007; Kent et al., 2018)111https://science.nrao.edu/facilities/vla/data-processing/pipeline. The pipeline calibrated measurement set had its weights reinitialized prior to imaging according to integration time and bandwidth. The imaging was performed using a customized version of CASA 5.3.0-123 for development of imaging algorithms by NRAO’s Algorithms Research and Development group. For imaging, the wideband AWProjection gridding algorithm with conjugate beam models was used in combination with multi-term multi-frequency synthesis (mtmfs) and multiscale clean algorithm (Rau & Cornwell, 2011; Bhatnagar et al., 2013; Rau et al., 2016). The “Briggs” weighting scheme (Briggs D., 1995, PhD Thesis, New Mexico Institute of Mining and Technology) was used with a robust setting of 0, corresponding to a compromise between uniform and natural weighting. For w-projection, 128 projection planes were used. The ray-traced primary antenna beam pattern was rotated with 5.0∘ steps with corresponding parallactic angle. The imaged area covered a region of 3.4${}^{\circ}\times$3.4∘ around the R.A./Dec. pointing of 00h02m41.88s and +62*∘18′2.2′′, that is well beyond the half-power point and the first null of the primary beam of the VLA. The rms noise in the final image near PSR J0002+6216 was 31 $\mu$Jy beam-1* with a synthesized beam of approximately 11*′′. This is about a factor of two above the expected theoretical thermal noise limit of 17$\mu$Jy beam-1*, not taking into account an increase of antenna temperature due to observing in the Galactic plane.

3 Timing Analysis

We constructed a coherent pulse timing model for PSR J0002$+$6216 using data from the Fermi Large Area Telescope (LAT, Atwood et al., 2009). We selected Pass 8 (P8R3 Atwood et al., 2013; Bruel et al., 2018) events collected from 2008 August 4 to 2018 November 12 with reconstructed position $<$2*∘* from the pulsar position and with reconstructed energy in the interval 100 MeV to 30 GeV. To improve sensitivity, we used a preliminary version of the 4FGL sky model222https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/ to compute the probability that each event originates from the pulsar rather than the background (“photon weights”, Kerr, 2011).

Because the integration time required to measure the pulse phase (“time of arrival”) is about one month, the time signature of proper motion is only marginally resolved. Thus we opted to use the PINT pulsar timing package333https://github.com/nanograv/PINT to perform an unbinned analysis. Moreover, even for large proper motions, the expected timing deviations are small and on the edge of detectability. Maximizing the sensitivity demands the sharpest possible model for the pulse profile, but determining the ideal pulse profile is not possible without knowing the timing solution in the first place. There is thus tension between sensitivity and bias for the proper motion detection, and we take great care to break this circular relationship by adopting the method we describe below.

To determine a good, sharp, starting model of the pulse profile, we optimized both the timing parameters ($\nu$, $\dot{\nu}$, and position) and a series of 12 Fourier components representing timing noise and other unmodelled signals (e.g. proper motion) by maximizing the $Z^{2}_{10}$ statistic (Buccheri et al., 1983), which is sensitive to the power in the first 10 Fourier components of the pulse profile. We believe this approach allows sufficient freedom in the timing model to produce the sharpest possible profile, but it avoids overfitting by restricting the harmonic content of the pulse profile444$\gamma$-ray Pulsars with 10 years of Fermi data typically require 40 or more Fourier modes to adequately model.. To the resulting phases we fit an analytic model, $f(\phi)$, comprising six wrapped Gaussians, shown in in Figure 1.

Next, we studied the timing noise by analyzing a series of models of increasing complexity, specifically including from 2 to 12 Fourier modes with frequencies of 1/Tobs, 2/Tobs, etc. We used emcee (Foreman-Mackey et al., 2013) to explore the posterior distribution of the likelihood $\log\mathcal{L}=\sum_{i}\log w_{i}f(\phi_{i},\lambda)+(1-w_{i})$, with the $w_{i}$ the photon weights and $\phi_{i}$ the rotational phases as determined by PINT for a timing model with parameters $\lambda$. (We adopt uniform priors on all parameters.) We find that the likelihood improves substantially with the addition of 3 Fourier modes, but that the best-fit power in higher modes drops immediately to a white noise floor. In tandem, the maximum likelihood only improves marginally with these additional degrees of freedom. Specifically, we observe a typical increase in the best-fit log likelihood of $\sim$2 for each additional Fourier component. Because these components satisfy Wilks’ Theorem (Wilks, 1938), such increases are formally insignificant. Similar conclusions arise from, e.g., the Akaike information criterion (Akaike, 1973). Thus, we adopt a three-Fourier component model as our baseline for further analysis. We note that the overall distribution of power in the Fourier modes is more complicated than the power-law relation seen in many other young pulsars (e.g., Kerr et al., 2015), showing an excess of power on few-year time scales that might originate from an unmodelled glitch recovery preceding Fermi observations.

We next added degrees of freedom for the proper motion to the timing model and studied the posterior over the range of timing noise models. For our preferred three-component model, the maximum likelihood improves by 5.7, which by Wilks’ Theorem has a chance probability of 0.0033, or 2.9$\sigma$ significance. From the samples of the posterior distribution for each model we computed the 16%, 50%, and 84% quantiles for the magnitude and position angle of the proper motion, and the measurement results are displayed in Figure 2. For our preferred model, we find a total proper motion of 115$\pm$33 mas yr*-1* at a position angle of 121°$\pm$13°. The samples for the coordinate proper motions are correlated ($-0.24$) such that more positive values of $\mu_{\alpha}$ prefer more negative values of $\mu_{\delta}$, leading to a narrower distribution of position angle than the one-dimensional marginalized distributions indicate. By analyzing the variation of the measurements over the timing noise models, we estimate a systematic uncertainty in the position angle of about 10°, while the magnitude is largely unaffected in comparison to the statistical uncertainty.

Finally, we assessed the effect of the assumed pulse profile model on the proper motion measurement by performing Monte Carlo Markov chain (MCMC) runs with the 18 parameters (3 for each Gaussian) of the model allowed to vary, thus marginalizing over these nuisance parameters. Examples of realizations of the template parameters are shown in Figure 1. We find similar results: the maximum log likelihood is improved by 5.1, with a chance probability of 0.006 (2.8$\sigma$), and a total proper motion of 112$\pm$39 mas yr*-1* at a position angle of 125°$\pm$19°.

These are large values for a pulsar proper motion. In the most recent ATNF Pulsar Catalog (Manchester et al., 2005), only 2% of the sample have measured proper motions as large, and most of these pulsars are at distances $<$1 kpc. We will explore the implications of this result in §5.

4 Results

Figure 3 shows our 20-cm (1.5 GHz) continuum image toward PSR J0002+6216. The pulsar sits at the tip of an elongated cometary tail with Galactic coordinates $l,b$=117.327*∘,$-0.074^{\circ}$. The observed emission is dominated by the cometary tail. The angular extent of this tail is at least 7′ and it remains uniformly bright and highly collimated along most of its length. The feature is unresolved (width$<$13′′) for the first half of its length, and is only marginally resolved (width$\sim$17′′-20′′) for the second half. There is extended emission (width$>$60′′*) that projects beyond the 7′ tail that may be associated with the cometary feature. However, the rms noise is 20% higher in this region and thus it is possible that this emission is associated with the nearby SNR CTB 1 instead.

The total flux density $S_{\nu}$ of this 7′ feature at 1.5 GHz is 14$\pm$0.9 mJy. In order to measure a spectral index, where $S_{\nu}\propto\nu^{\alpha}$, we looked for radio images at different frequencies from archival surveys (Rengelink et al., 1997; Intema et al., 2017), and pointed observations of CTB 1 (e.g., Dickel & Willis, 1980) but none of these had the requite sensitivity, resolution or field-of-view. We attempted an in-band spectral index measurement by splitting our 1 GHz band into two halves, centered at 1.25 GHz and 1.75 GHz, and imaging and deconvolving separately. The resulting value $\alpha=-0.98\pm 0.31$ is not very precise but it is sufficient to suggest that the emission is non-thermal in origin. We imaged the field with full Stokes parameters, but the rms noise is too high to put meaningful limits on the degree of polarization ($<$24%).

We also looked for evidence of the tail-like feature at other wavelengths. There is an X-ray point source whose position is consistent with PSR J0002+6216 (Wu et al., 2018; Zyuzin et al., 2018). The cometary feature is not visible in ROSAT images (Craig et al., 1997), but it lies outside of the field-of-view of deeper pointings made by the narrow-field instruments of ASCA, Chandra and Suzaku (Lazendic & Slane, 2006; Pannuti et al., 2010; Katsuragawa et al., 2018). Nor is there a structure with this morphology in deep H$\alpha$, [N ii], [S ii], or [O iii] images (Fesen et al., 1997).

Based on this morphology and the non-thermal spectral index, we identify G 117.33$-$0.07 as a sub-class of supersonic or bow-shock pulsar wind nebulae (PWNe; Kothes, 2017; Kargaltsev et al., 2017; Slane, 2017). These PWNe are typically found after the pulsar has escaped the high pressure environs of their parent SNR. The bow shock is formed as the relativistic pulsar wind is shocked and confined by ram pressure due to the high space velocity of the pulsar through the interstellar medium (ISM). The shocked particles and magnetic energy are swept backward where they emit broadband, non-thermal synchrotron trailing along in a “tail”.

The PWN is detected close to the noise threshold in Stokes I images at 1.42 GHz from the Canadian Galactic Plane Survey (CGPS; Kothes et al., 2006). We show this CGPS image of CTB 1 in Figure 4. While the resolution is only 1′, the CGPS image includes both single-dish and interferometric data, and thus it faithfully reproduces angular structure on all scales up the resolution limit. The tail in this CGPS image extends from PSR J0002+6216 unbroken 11′ to the southeastern edge of the SNR CTB 1. A line from the pulsar along the tail appears to point back to the center of the remnant 28′$\pm$1′ away, suggesting a possible common origin. To test this hypothesis we measured the positions of 10 radio peaks along the 7′ feature and we fit a linear least squares solution to a line of these data. The extrapolation of this 7′ line passes within 5*′′* and 11*′′* of the geometric center of CTB 1 independently derived by Landecker et al. ($\alpha_{1950}=23^{\rm h}45^{\rm m}45\fs$, $\delta_{1950}=62^{\circ}10.5^{\prime}$; 1982) and Kothes et al. ($\alpha_{2000}=23^{\rm h}59^{\rm m}16\fs$, $\delta_{2000}=62^{\circ}27^{\prime}$; 2006), respectively. This offset is considerably smaller than the $\pm$1′ position uncertainty of the geometric center of CTB 1. Note that we are assuming here that the geometric center is a good proxy for the origin of the supernova event, which is not necessarily the case (Gvaramadze, 2004).

A further test of the association is to look for proper motion of the pulsar. The age of CTB 1, derived from X-rays, assuming Sedov evolution, is 9-13.3 kyr (Craig et al., 1997; Lazendic & Slane, 2006), while similar values of 7.5 kyr (Hailey & Craig, 1994) are obtained from optical data. We adopt a value of 10 kyr, recognizing that this value is uncertain by about 20%. With the 28′ PSR-SNR offset555Zyuzin et al. (2018) use a PSR-SNR angular offset of 17′ rather than 28′, apparently in error. we predict the proper motion of PSR J0002+6216 to be 168$\pm$35 mas yr*-1* in the southeast direction with a position angle $\theta_{\mu}=113^{\circ}$ (measured from north through east). This is in good agreement with pulsar timing proper motion measurement of $\mu=115\pm 33$ mas yr*-1* at a position angle $\theta_{\mu}=121\pm{13}^{\circ}$. Moreover, we can rule out a proper motion with magnitude $<$63 mas yr*-1* with 95% confidence. If we constrain the proper motion position angle to be $113\pm 3^{\circ}$, we find $\mu=121\pm 30$ mas yr*-1* and the 95% confidence range increases to $>$73 mas yr*-1*. We conclude that PSR J0002$+$6216 passed through the center of CTB 1 between 10 and 20 kyr ago.

5 Discussion

In §4 we have shown that PSR J0002+6216 lies at the “head” of a cometary-shaped feature whose tail points back to the geometric center of the SNR CTB 1. We argue that this feature is a bow-shock PWN shaped by the pulsar’s supersonic motion through the ISM, and that the PSR and SNR may have a common origin. From a pulsar timing analysis we measure a proper motion that agrees in magnitude and direction to the predicted value derived from the PSR-SNR angular offset and the SNR age. In the following subsections we will explore the implications of each of these main results.

5.1 The velocity and distance to PSR J0002+6216

Our measured proper motion from the timing analysis (§3) of $\mu=115\pm$ 33 mas yr*-1* corresponds to a transverse pulsar velocity $V_{PSR}$=550$\times{d_{kpc}}$ km s -1, where $d_{kpc}$ is the pulsar distance in kpc. The distance to PSR J0002+6216 can be estimated from its dispersion measure DM=218.6 pc cm*-3* (Wu et al., 2018). The two models for the Galactic distribution of ionized gas give $d$=7.9 kpc (Cordes & Lazio, 2002) and $d$=6.4 kpc (Yao et al., 2017), respectively, leading to a suspiciously large $V_{PSR}$=3500-4300 km s -1. Both Wu et al. (2018) and Zyuzin et al. (2018) have also noted that the DM distance would also require that the $\gamma$-ray efficiency of PSR J0002+6216 exceed unity. They derive a pulsar distance from an empirical relationship between the $\gamma$-ray luminosity and $\dot{{E}}$ of 2.0 kpc and 2.3 kpc, respectively. At $d$=2.3 kpc our measured proper motion gives $V_{PSR}$=1260$\pm$ 360 km s -1. While still large, there are other pulsars with well-measured proper motions with comparable values of $V_{PSR}$ including PSRs B1727$-$47, B2011+38, and B2224+65 (Hobbs et al., 2005; Shternin et al., 2017).

The foregoing suggests that the DM-based distance to PSR J0002+6216 may be an overestimate and that there is an excess source of free electrons beyond that assumed in the Galactic models (Cordes & Lazio, 2002; Yao et al., 2017). Following Kulkarni et al. (2014) we can attempt to derive the size, $L$, and distance, $d$, to this excess contribution $DM_{e}$ using three observational constraints: the diffuse H$\alpha$ emission, the known ionizing stars, and the DM distribution of local pulsars. At an angular resolution of 1∘ the Wisconsin H$\alpha$ Mapper Northern Sky Survey (WHAM; Haffner et al., 2003) measures an integrated H$\alpha$ intensity of 28.5 Rayleigh, or EM=65 pc cm*-6*, where we have used the conversion of 1 Rayleigh corresponding to EM=2.25 pc cm*-6* for an ionized gas with electron temperature of 8000 K. We can obtain a second estimate of EM on arcminute scales from the deep H$\alpha$ image of Fesen et al. (1997). Along the southeastern edge of the nearby CTB 1 they measure an extinction-corrected H$\alpha$ surface brightness 6.6$\times 10^{-16}$ erg cm*-2* s*-1* arcsec*-1* or 2.8$\times 10^{-5}$ erg cm*-2* s*-1* sr*-1*, or EM=260 pc cm*-6*, where we have used the conversion of 1 Rayleigh=2.41$\times 10^{-7}$ erg cm*-2* s*-1* sr*-1* (Kulkarni et al., 2014). From the same deep H$\alpha$ images we estimate that the diffuse H$\alpha$ emission toward PSR J0002+6216 is about 10$\times$ fainter, or EM=30 pc cm*-6*. Thus any nebular source of excess electrons is constrained to have EM$\leq$30-65 pc cm*-6*, where EM=DM${}_{e}^{2}$/$L_{pc}$ and $DM_{e}$=n_{e}\times$$L_{pc}.

It is straightforward to account for this excess contribution if the pulsar lies at a distance $d\gtrsim$3.4 kpc. At that distance the excess DM required in the Yao et al. (2017) model over the observed DM for PSR J0002+6216 is DM${}_{e}\simeq$50 pc cm*-3* and EM=2500/$L_{pc}$. $DM_{e}$ is about a factor of two higher for the Cordes & Lazio (2002) model. Neither of these models account for the fact that PSR J0002+6216 passes within 8′ of the line of sight of HD 225160, an 8th magnitude O8 blue supergiant with a parallax distance from Gaia DR2 of 3.4$\pm$0.4 kpc (Gaia Collaboration et al., 2018). Such a star will ionize an extended region and be surrounded by a Strömgren sphere whose size $R_{s}$ for a uniform density $n_{e}$ is given by $R_{s}$=$n_{e}^{-2/3}U$. The excitation parameter $U$ (in units of pc cm*-2*) is calculated from stellar atmospheric models and conveniently expresses the ionizing flux as a function of spectral type and spectral class (e.g., Panagia, 1973). A line-of-sight that intersects an H ii region with impact parameter $R$ will see an excess dispersion measure $DM_{e}$=2$R_{s}{n}_{e}\sqrt{1-({R}/{R}_{s})^{2}}$ where $L=2\sqrt{{R}_{s}^{2}-{R}^{2}}$ (Prentice & Ter Haar, 1969). For ISM number densities $n_{e}<30$ cm*-3* the line-of-sight of the pulsar ($R$=8 pc, i.e., 8′ at 3.4 kpc) always passes through the Strömgren sphere of HD 225160. The derived value of $DM_{e}$ is relatively insensitive to ne with $DM_{e}$=70-150 pc cm*-3* for ne=0.1 to 1.0 cm*-3* and the corresponding emission measures are 10-150 pc cm*-6*. Thus an H ii region ionized by HD 225160 can easily explain the required $DM_{e}$ and is consistent with the EM constraints from H$\alpha$.

At smaller distances a third constraint comes into play since the angular size of the nebula providing $DM_{e}$ cannot be so large as to affect the DMs of known pulsars in the vicinity. Within a 6∘ radius of PSR J0002+6216 there are 12 other pulsars. The majority (9) have DMs in the range of 100-125 pc cm*-3* with an average distance of 2.4 kpc, while the remaining have DMs of about 200 pc cm*-3*. PSR B2351+61 with DM=94 pc cm*-3* is closest on the sky to PSR J0002+6216, which at 1.1∘ provides a constraint on $L/d$. For example, at $d$=2.3 kpc the maximum dimension the nebular region can have in the plane of the sky is $L$=44 pc. At this distance the Yao et al. (2017) model predicts DM${}_{e}\simeq$130 pc cm*-3* for PSR J0002+6216 and EM=1.7$\times{10^{4}}/L$. To satisfy the deep H$\alpha$ EM constraints (EM$\leq$30-65 pc cm*-6*) any putative nebula would need to be significantly elongated (10:1) or there must be large extinction toward PSR J0002+6216, reducing the observed EM. We examined 3D extinction models in this direction (e.g., Sale et al., 2014). Within 12′ of PSR J0002+6216 the extinction (defined at $\lambda$=549.5 nm) has a value A*∘=1.3 mag by $d$=2.3 kpc and nearly doubles by $d$=3.5 kpc. The distribution of dust is quite patchy; within a 1∘ radius $A_{\circ}$ at 1 kpc varies from 0.43 mag to 1.38 mag, while at 3 kpc the extinction range is 1.2 to 2.4 mag. This is consistent with the observed column density and optical extinction toward the nearby SNR CTB 1, varying from $A_{V}$=1 to 2 mag (Fesen et al., 1997; Katsuragawa et al., 2018). It becomes exceedingly difficult, however, to satisfy all constraints at smaller distances. At $d$=1 kpc, for example, PSR J0002+6216 is less an outlier in the pulsar velocity distribution, but the local DM constraints give $L\leq$20 pc and the Yao et al. (2017) model predicts DM${}_{e}\simeq$200 pc cm-3* and EM=4$\times{10^{4}}/L$. To satisfy the EM constraints would require a narrow tube of nearly 1000 pc in length and width 20 pc, a highly unlikely configuration.

Summarizing, we find that the DM-distance gives unreasonable values for $V_{PSR}$ and requires that the pulsar $\gamma$-ray efficiency exceeds unity. We argue that there is evidence that there are additional source of free electrons along the line of sight. The nearest distance that is consistent with the existing observational constraints is $d$=2 kpc. This is the same distance derived from an empirical relation between between the $\gamma$-ray luminosity and $\dot{{E}}$. For the remaining discussion and derivations we will adopt $d$=2 kpc ($d_{2}\equiv$1) for PSR J0002+6216, recognizing that this estimate is uncertain.

5.2 The PWN of PSR J0002+6216

Within the sub-class of bow-shock or supersonic PWNe there are a diverse range of morphologies (see list in Kargaltsev et al., 2017), likely driven by small-scale variations in the density and magnetic field of the ISM, and/or the different orientations of the pulsar magnetospheric spin axis with respect to the velocity vector (Bucciantini, 2014; Barkov & Lyutikov, 2018; Toropina et al., 2018). The PWN that most resembles G 117.33$-$0.07 is G 315.78$-$0.23 around PSR J1437$-$5959, a moderately aged pulsar with a period of 61.7 ms and a spin-down energy $\dot{{E}}$=1.5 $\times 10^{36}$ erg s*-1* (Camilo et al., 2009; Ng et al., 2012). Both PWNe have the same angular extent with narrow heads and a narrow, uniformly bright tail. Here we will follow Ng et al. (2012) in calculating some standard PWN properties for G 117.33$-$0.07 and compare them to G 315.78$-$0.23. Our estimates should be viewed as preliminary since we have only a single radio image (1-2 GHz) in which G 117.33$-$0.07 is only marginally resolved (Figure 3) and the synchrotron spectrum is poorly constrained (§4).

The size of the bow shock region $r_{s}$ can be estimated by equating the ram pressure of the fast-moving pulsar $\rho{V}^{2}_{PSR}$ to the wind pressure of the pulsar $\dot{{E}}$/4$\pi r^{2}_{s}$c (see Frail et al., 1996). Since the source is unresolved we assume an isotropic wind with $\dot{{E}}$=1.53 $\times 10^{35}$ erg s*-1* for PSR J0002+6216 and adopt ISM number densities of 0.1-1 cm*-3*, which bracket the density range derived from HI and X-ray observations (e.g., Yar-Uyaniker et al., 2004; Lazendic & Slane, 2006). We have assumed a 10% He mass in the ISM when calculating $\rho$. From this we derive $r_{s}$=(1-4)$\times 10^{-3}d_{2}^{-1}$ pc. This is identical to the small standoff distance of $r_{s}\simeq 2.4\times 10^{-3}$ pc derived for G 315.78$-$0.23, but since G 117.33$-$0.07 is 4x closer it might be possible to resolve its bow shock (0.1-0.4*′′) with future radio observations. An alternative way to express this pressure balance is to write it in terms of the ISM pressure (Kargaltsev et al., 2008), which is generally more robustly determined than the ISM density. Now $\rho{V}^{2}_{PSR}$ is re-written as $\gamma{\mathcal{M}^{2}}P_{ISM}$, where $\gamma\equiv 5/3$ is the ISM adiabatic index, $\mathcal{M}$ is the Mach number, and $P_{ISM}$ is a typical thermal ISM pressure of 10-12* dyn cm*-2*. The estimated shock Mach numbers are $\mathcal{M}\simeq{200}$ for both G 315.78$-$0.23, and G 117.33$-$0.07. These large Mach numbers and the narrowness of the PWN tails are suggestive of a high pulsar velocity in both systems.

While the head of G 117.33$-$0.07 is unresolved with our 12*′′* beam, the linear dimensions of the 7′ tail are $l$=4.1$d_{2}$ pc and the half-width of the end of the tail is $h$=0.09$d_{2}$ pc. This implies that the PWN has expanded by about a factor of 45 from the tip of the head (rs) to the end of the tail ($h$). At $d$=8 kpc the corresponding dimensions of G 315.78$-$0.23 are a tail length of $l$=20 pc and half-width of $h$=1.2 pc. Provided that the distance of 8 kpc is correct, G 315.78$-$0.23 has the longest radio tail among the known bow-shock PWNe, while G 117.33$-$0.07 has an average tail length (Kargaltsev et al., 2008; Ng et al., 2012; Kargaltsev et al., 2017). We approximate the volume of G 117.33$-$0.07 using a half-cone of radius $r_{s}$ plus a frustum of a right circular cone with length $l$ and radii $r_{s}$ and $h$, giving a volume $V=10^{54}d_{2}^{3}$ cm*-3*.

The slope of the radio spectrum of G 117.33$-$0.07 of $\alpha=-0.98\pm 0.31$ (§4), is atypical of radio PWNe and is more similar to the slopes of X-ray PWNe (Kothes, 2017; Slane, 2017). We caution however that $\alpha$ was estimated within a narrow 1 to 2 GHz band and needs to be better measured over a wider frequency range. For now we will calculate PWN parameters with both $\alpha=-1$ and a more typical radio PWN of $\alpha=-0.3$. The radio luminosity $L_{R}$ is obtained by integrating from 10 MHz to 100 GHz, giving values of 2.7$\times{10}^{30}d_{2}^{2}$ erg s*-1* ($\alpha=-0.3$) and 9.3$\times{10}^{29}d_{2}^{2}$ erg s*-1* ($\alpha=-1$). The radio efficiency $L_{R}$/$\dot{{E}}$ is $\sim$10*-5*, typical of other radio PWNe (Frail & Scharringhausen, 1997; Gaensler et al., 2000).

The synchrotron spectrum can be used to derive useful PWN parameters such as the strength and energy density of the magnetic field, the energy density of the relativistic particles, and the lifetime of the radiating electrons. The standard expressions are given in Pacholczyk (1970) but we use equations 8-10 from Ng et al. (2012) in which the magnetic field B, the synchrotron lifetime tsyn and the relativistic gas and field pressure in the tail $P_{tail}$ are conveniently expressed in terms of the lesser-known wind parameters such as the ion to electron energy density ratio $\eta$ and magnetic to particle (electrons and ions) energy density ratio $k_{m}$. If we adopt their values of $\eta$=0 (i.e. no ions) and $k_{m}$=0.1 (electron-dominated PWN) we derive $B$=30$d_{2}^{-2/7}$ $\mu$G, $t_{syn}$=1$d_{2}^{3/7}$ Myr (at 10 GHz), and $P_{tail}$=1.6$\times{10}^{-10}{d_{2}^{-4/7}}$ dyn cm*-2*. Remarkably, apart from a distance scaling, these values for G 117.33$-$0.07 are identical to the values derived for G 315.78$-$0.23. We have integrated the synchrotron spectrum in the range 107 Hz to 1013 Hz, but the values are relatively insensitive to the outer frequency range and they are not sensitive to which value of $\alpha$ we use.

G 117.33$-$0.07 and G 315.78$-$0.23 appear to be near morphological twins of each other, although the tail of G 117.33$-$0.07 is five times smaller at our adopted distance of 2 kpc. Both PWNe are powered by high-energy pulsars moving at high velocity away from a parent SNR (see §5.3). Their PWNe are characterized by a small shock standoff distance and a large Mach number ($\mathcal{M}\simeq 200$) with long synchrotron-emitting tails which are likely dominated by relativistic particles rather than B-fields. In the case of G 315.78$-$0.23 the magnetic field is detected and is seen to be aligned with the tail. For G 117.33$-$0.07 we lack such data. Future broadband radio continuum and polarimetric observations will fully resolve the tail and measure the magnetic field orientation, and together with X-rays will better constrain the spectrum and thus the parameters derived above. Likewise, higher angular resolution observations of the head of G 117.33$-$0.07 would also show whether the PWN morphology is being shaped by the geometry of the pulsar magnetopshere and/or the pulsar proper motion.

The age of the radio PWN is obtained simply by dividing its angular size with the pulsar proper motion giving $t_{R}$=3600 yr. As expected, $t_{R}$ is much smaller than the 1 Myr radiative lifetime of the radio-emitting electrons but at X-ray energies (assuming B=30 $\mu$G from above) $t_{syn}\simeq$200 yr (1 keV). Thus we might expect to see a smaller X-ray PWN scaled approximately by the ratio of 7′($t_{syn}$/$t_{R})\simeq$20*′′, or even longer if the streaming velocity of the post-shock wind is a significant fraction of the speed of light. The prospects for the detection of an X-ray PWN appear good since for typical X-ray PWN efficiencies of 10-3* the X-ray luminosity $L_{X}\simeq{10}^{32}$ erg s*-1* (Kargaltsev et al., 2008), and G 117.33$-$0.07 is relatively nearby at 2 kpc with only modest gas column densities expected $\sim$1021 cm*-2*.

5.3 The PSR J0002+6216 SNR CTB 1 association

All neutron stars, given a substantial kick at the time of birth, will eventually escape their parent SNR on a timescale $\tau_{esc}=({E}_{\circ}/{n}_{\circ})^{1/3}$ $V^{-5/3}_{PSR}$ that is only weakly dependent on the explosion kinetic energy E∘ and the ISM density n∘, but is sensitive to the magnitude $V_{PSR}$ (van der Swaluw et al., 2003; Bykov et al., 2017). Our measured proper motion of PSR J0002+6216 is of the right magnitude and direction to support the claim that it was born from the same SN that produced the SNR CTB 1 (Zyuzin et al., 2018). Likewise, the direction and the morphology of the 7-11′ (4.1-6.4 pc) tail of the PWN suggest a physical connection between a high-velocity pulsar and its SNR. While promising, a secure PSR-SNR association also requires good agreement of distances and ages (Kaspi, 1996).

The distance to CTB 1 can be reliably estimated since it is claimed to be part of a much large star-forming complex in the Perseus arm (Fich, 1986). On smaller scales CTB 1 is embedded in a neutral hydrogen (H i) hole and is interacting with neutral gas along its bright radio continuum edges at a local standard of rest velocity of $-30$ km s*-1* (Landecker et al., 1982). This value agrees remarkably well with the mean velocity derived from optical spectroscopy (Hailey & Craig, 1994, and references therein). Thus we adopt the kinematic distance of 2.0$\pm$0.4 kpc from Landecker et al. (1982), but we note that Yar-Uyaniker et al. (2004) find the same H i structures and velocities but argue instead that CTB 1 is a blue-shifted Local arm object at a distance of 1.6 kpc. Larger kinematic distances quoted for CTB 1 use older IAU Galactic rotation parameters (Reid et al., 2014; Wenger et al., 2018), or neglect to account for the well-known non-circular motions toward the Perseus arm in this direction (Choi et al., 2014; Sakai et al., 2018).

In §5.1 we argued that the nominal DM distance to PSR J0002+6216 was an overestimate and we showed that there were additional sources of free electrons along the line of sight. The nearest distance that PSR J0002+6216 could be without violating several observational constraints is $d$=2 kpc, although we can not rule out $d\simeq$3.5 kpc. If PSR J0002+6216 lies at 2 kpc its distance would agree with the kinematic distance of SNR CTB 1. We note that this is the same distance as the adjacent PSR B2334+61 (DM=58.41 pc cm*-3*) SNR G 114.3+0.3 association (Kulkarni et al., 1993) 3∘ away, and it would place both PSR-SNR associations in the same star-forming complex within the Perseus arm.

The characteristic age of the pulsar $\tau_{c}$=306 kyr greatly exceeds the SNR age of $\sim$10 kyr by a factor of 30. If $\tau_{c}$ is the correct age, our proper motion value would suggest that the pulsar was born far away (11∘) and that the PWN and CTB 1 are a chance alignment. Ng et al. (2012) argue that the chance alignment of a PWN tail pointing back to the geometric center of an SNR is exceedingly unlikely. These two ages could be reconciled if PSR J0002+6216 was born with an initial period (P∘) close to its current period P=115 ms (Camilo et al., 1994; Zyuzin et al., 2018). If so, PSR J0002+6216 would be one of a growing number of young PSRs in SNRs (including PSR J1437$-$5959; Camilo et al., 2009) for which P∘/P $\not\ll$1 (Popov & Turolla, 2012), casting doubt on the reliability of $\tau_{c}$ for assessing PSR-SNR associations.

Comparing independently derived pulsar velocities provides another consistency check on the association. At a distance of 2$\pm{0.4}$ kpc and age 10$\pm{0.2}$ kyrs, the 28′$\pm$1′ PSR-SNR offset predicts $V_{PSR}$=1600$\pm$450 km s -1, where the uncertainty is calculated from the quadrature sum of the errors on age, distance and angular offset. From our measured proper motion and the 2 kpc distance we derive an age-independent $V_{PSR}$=1100$\pm$315 km s*-1*. A third method for estimating the transverse velocity of the pulsar is given by Frail et al. (1996) where the velocity of the pulsar can be written as $V_{PSR}$=$V_{s}\,\beta$/$c_{\circ}$, where $V_{s}$ is the shock velocity of the SNR, $\beta$ is the fractional offset of the pulsar from the center of the remnant normalized by the remnant radius, and $c_{\circ}$ is a constant equal to 2/5 for a remnant in the Sedov phase. This method has the advantage of being independent of both the age and the distance of the association, but it does weakly depend on the evolutionary state ($c_{\circ}$) of the remnant. A lower bound of $V_{s}\gg$ 100 km s*-1* comes from optical line ratios such as [O iii]/H$\beta$ (Fesen et al., 1997). Thermal plasma is seen from CTB 1 with X-ray temperatures of 0.2 keV to 0.28 keV, which corresponds to Vs=410-480 km s*-1* (see eqn. 2 of Lazendic & Slane, 2006), and thus $V_{PSR}$=1780$\pm$150 km s -1. The consistency of these three $V_{PSR}$ estimates does bolster the claim of a PSR-SNR association.

5.4 PSR J0002+6216 and pulsar kick mechanisms

One of the more robust conclusions from this work is that PSR J0002+6216 is a high-velocity pulsar. Most pulsars have average transverse velocities of order 250 km s*-1* (Hobbs et al., 2005; Verbunt et al., 2017), so PSR J0002+6216 is a rare outlier with $V_{PSR}>$1000 km s*-1*. While it has long been known that pulsars receive a substantial kick at birth, the debate about the mechanism(s) is still an active research topic. Numerous kick mechanisms have been proposed including binary disruption, asymmetric neutrino emission, jets that accelerate the pulsar, and hydrodynamic instabilities (Iben & Tutukov, 1996; Lai et al., 2001; Scheck et al., 2006; Ng & Romani, 2007; Wongwathanarat et al., 2013; Kochanek et al., 2018). A high-velocity pulsar like PSR J0002+6216 poses a strong challenge to these models. From the magnitude of the velocity alone we can rule out binary disruption (Iben & Tutukov, 1996; Fryer et al., 1998). Some jet models predict large scale morphological distortions in the SNR along the axis defined by the pulsar’s proper motion (Bear & Soker, 2018). This is not seen in the case of SNR CTB 1 which is remarkably circularly symmetric in both optical H$\alpha$ and non-thermal radio.

A general prediction for several natal kick models, that has some observational support (Johnston et al., 2007; Ng & Romani, 2007), is an alignment between the pulsar’s rotation axis and the direction of the pulsar velocity (Spruit & Phinney, 1998). Wu et al. (2018) have fit the pulse profiles of PSR J0002+6216 to constrain some of the magnetospheric parameters using outer gap (OG) and slot gap emission models (TPC). They measure the angle between the rotation axis and the line of sight $\xi_{TPC}=54^{\circ}\pm 2\arcdeg$ and $\xi_{OG}=58\arcdeg^{+25}_{-1}$, and the angle between the rotation axis and the magnetic axis $\alpha_{PC}=64\arcdeg^{+3}_{-2}$ and $\alpha_{OG}=69\arcdeg^{+8}_{-1}$. These estimates are highly model dependent and are subject to systematic errors of 10∘ or larger, so that no strong conclusions should be drawn from these values.

The angle $\xi$ is particularly interesting since for those PWNe with a toroidal wind it describes the angle that the axis of the torus makes with respect to the plane of the sky (Helfand et al., 2001; Ng & Romani, 2004), and it can have a strong impact on the morphology of bow-shock PWNe (Barkov & Lyutikov, 2018). Unfortunately, while we know $\xi$, we do not know its position angle on the plane of the sky to compare with $\theta_{\mu}$, the direction of the pulsar’s motion. This angle is derived from polarization measurements of the radio pulse profile and PSR J0002+6216 is too faint to make such measurements in the near term. However, if the velocity and rotation axis are aligned then we can say that the 3D velocity would be $(\cos\xi)^{-1}\simeq 1.7\times$ larger than the transverse value $V_{PSR}$. Another trend expected from these finite duration kick models is that the fastest moving pulsars should have the longest birth periods (Figure 12 of Ng & Romani, 2007). Interestingly for our high-velocity PSR J0002+6216, we agree with Zyuzin et al. (2018) that if the PSR-SNR association is real, the pulsar birth period P${}_{\circ}\simeq$P=115 ms (§5.3).

Recent work on young ejecta-dominated SNRs has established a link between pulsar kick velocities and SNR asymmetries, as measured via X-ray morphologies (Katsuda et al., 2018). These results support the hydrodynamic instability kick model in which the SN recoil expels the pulsar in one direction and the newly-synthesized heavy elements from the core in the opposite direction (Janka, 2017). The high velocity of PSR J0002+6216 could pose a particularly strong test of this mechanism. To illustrate this we note that the kinetic energy of a 1.5 M⊙ neutron star moving with $V_{PSR}$=1100 km s*-1* is $E_{k}=0.2\times{10}^{50}$ erg. This is likely a lower limit since we have only the transverse 2D velocity of the pulsar. Several different analyses have all concluded that CTB 1 is a low-energy SNR with an explosion kinetic energy $E_{\circ}\simeq{10}^{50}$ erg (Craig et al., 1997; Fesen et al., 1997; Lazendic & Slane, 2006). As $E_{k}$ is a significant fraction of E∘ it may be difficult for a low-energy SN like CTB 1 to impart a substantial natal kick to PSR J0002+6216, but if the mechanism does work it would predict strong asymmetries.

Is there a signature from the natal kick given to PSR J0002+6216 in either the large-scale X-ray morphology of CTB 1 or the distribution of its ejecta? The X-ray emission from CTB 1 is concentrated within the radio shell with its X-ray centroid shifted to the northwest (e.g., Pannuti et al., 2010). This asymmetry could be due to the natal recoil, but since CTB 1 is a middle-aged SNR, it is more likely that the morphology is dominated by its interaction with the ISM (Lazendic & Slane, 2006). There is ample evidence from the H i kinematics that CTB 1 is strongly interacting with the surrounding gas (Landecker et al., 1982), and from X-ray, H$\alpha$ and H i velocity images we see that the northeast rim of the remnant has broken out of this cavity and is expanding into a lower density region (Craig et al., 1997; Fesen et al., 1997; Yar-Uyaniker et al., 2004). Thus we find no evidence to suggest that a kick signature is seen in the morphology of CTB 1.

Despite CTB 1 being a mixed morphology remnant that is strongly shaped by the ISM, the analysis of X-ray spectra shows evidence for enhanced heavy element abundances (Lazendic & Slane, 2006). The abundance ratios of these ejecta are consistent with CTB 1 being an oxygen-rich SNR produced in a core-collapse SN with a progenitor mass of 13-15 M⊙ (Pannuti et al., 2010). In a deep (82 ks), narrow-field Suzaku pointing of the center of CTB 1, Katsuragawa et al. (2018) carried out a spectral analysis of three regions, finding enhanced iron abundances in one region (region E) that they interpret as originating from asymmetric SN ejecta. This region encompasses a large solid angle of the SNR, so it is not possible to say at this time if the centroid of this ejecta asymmetry is directed away from the direction of motion of PSR J0002+6216. Future searches for heavy element asymmetries should be made along an axis defined by the PWN tail and the geometric center of the SNR.

6 Summary and Conclusions

In this paper we have presented an analysis of VLA 20-cm radio continuum observations and Fermi timing data toward the $\gamma$-ray and radio pulsar PSR J0002+6216. We discovered a bow-shock PWN G 117.33$-$0.07 with a long (4 pc), narrowly collimated non-thermal tail. We estimated the physical properties of G 117.33$-$0.07 and found that it was strikingly similar to another bow-shock PWN G 315.78$-$0.23 (Ng et al., 2012). They are both shaped by the large spin-down energy and high (inferred) pulsar velocities, leading to high Mach number shocks with small, compact heads and long synchrotron-emitting tails. In the case of G 117.33$-$0.07 we have measured the pulsar proper motion and confirmed the high velocity.

We found that the tail of G 117.33$-$0.07 points back toward the geometric center of the SNR CTB 1, suggesting a physical association that is similar to one claimed for the pulsar PSR J1437$-$5959, its PWN G 315.78$-$0.23 and the SNR G 315.9$-$0.0 (aka The Frying Pan; Ng et al., 2012). Our proper motion measurement is of the right magnitude and direction to support the claim that PSR J0002+6216 is a high-velocity pulsar that has escaped its parent SNR CTB 1. Nonetheless, problems remain and the association, while plausible, is not yet demonstrated. The DM-based distance for  J0002+6216 of  7 kpc is inconsistent with the CTB 1 distance, and its characteristic age of 306 kyr greatly exceeds the SNR age of 10 kyr. In §5 we have shown how this age and distance for PSR J0002+6216 are likely gross overestimates, but it remains the case that a secure PSR-SNR association requires agreement of independently measured distances and ages.

A more robust inference, independently arrived at from the pulsar proper motion and the PWN properties, is that PSR J0002+6216 is a high-velocity pulsar ($V_{PSR}>$1000 km s*-1*). This conclusion does not require the PSR-SNR association. We have adopted a conservative distance of 2 kpc, while a larger distance as inferred from the DM would only serve to increase $V_{PSR}$. We looked at different kick mechanisms for the origin of this very high pulsar velocity. While we lack the data to test whether the rotation axis of the pulsar is aligned with the velocity vector, we note that the large birth period estimated for a high-velocity pulsar like PSR J0002+6216 is a general outcome of some impulsive kick models. There is evidence that the heavy elements from CTB 1 are ejected asymmetrically, in support of the hydrodynamic kick model. However, a reanalysis of the Suzaku data could test whether the ejecta are preferentially found with direction opposite the pulsar’s motion.

The main limitation of this present work is the lack of a well-determined pulsar distance. A radio pulsar parallax measurement would settle the distance debate while also confirming the proper motion obtained by pulsar timing. The analysis of the PWN properties would be substantially improved with more radio data spanning a larger frequency range, and with greater angular resolution. Spectrally resolved X-ray and radio images could be used to help understand why the tail is so uniformly bright over most its length by looking for signatures of in situ particle acceleration of the shocked wind or evidence of synchrotron cooling. Of more immediate importance would be an independent measure of the pulsar velocity and the PWN from deep H$\alpha$ imaging and spectroscopy (e.g., Romani et al., 2010). In this regard PSR J0002+6216 has the same distance and large velocity as PSR B2224+65, known for its prominent H$\alpha$ “Guitar” nebula (Chatterjee & Cordes, 2004). Scaling relations from Brownsberger & Romani (2014) suggest that the putative bow shock nebula would be compact ($\theta\propto\dot{{E}}^{1/2}/d^{2}$) but bright ($f_{{\rm H}\alpha}\propto\dot{{E}}$/$d^{7/2}$), provided that the larger $\dot{E}$ of PSR J0002+6216 (1.53$\times 10^{35}$ erg s*-1* vs 1.2 $\times 10^{33}$ erg s*-1*) does not ionize the H i in the surrounding ISM (Landecker et al., 1982; Yar-Uyaniker et al., 2004).

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. This research has made use the High Energy SNR catalog of Ferrand & Safi-Harb (2012)666http://www.physics.umanitoba.ca/snr/SNRcat/. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The research presented in this paper has used data from the Canadian Galactic Plane Survey, a Canadian project with international partners, supported by the Natural Sciences and Engineering Resources Council. Work at NRL is supported by NASA. DAF thanks Paul Demorest, Shri Kulkarni, Sterl Phinney and Roger Romani for useful discussions. MK thanks David Smith for pointing out the connection with HD 225160. We thank the anonymous referee for useful comments that improved the quality of this manuscript.