© M. G. F. Mayer and W. Becker 2021

1. Introduction

It is generally accepted that collapsing massive stars (i.e., core-collapse supernovae) leave behind a compact remnant, a neutron star (NS) or black hole. A natural consequence of asymmetries in core-collapse supernovae is a strong kick acting on the compact remnant (e.g., Wongwathanarat et al. 2013), whose kinematics are thus directly connected to explosion properties. Proper motion studies of NSs are therefore a powerful tool that allows a lower limit to be set on the kinetic energy and momentum transferred to the NS during the explosion of the progenitor. Beyond that, depending on its age, the proper motion of an NS may reveal its origin or even the exact supernova explosion site, which can be used to robustly infer the age of the NS and/or associated supernova remnant (SNR).

For radio pulsars, which make up the vast majority of known NSs, proper motion can be measured precisely for a comparatively large number of sources, using precise positions from pulsar timing or very-long-baseline interferometry. This allows for population studies of NS kinematics, as for example by Hobbs et al. (2005) and Verbunt et al. (2017). Such works establish the picture of NSs being very high-speed objects, with a mean three-dimensional velocity around 400 km s⁻¹, and reliably measured projected velocities at least as high as ∼ 800 km s⁻¹.

In contrast to radio pulsars, for NSs exhibiting at best weak radio and optical emission – such as magnetars, X-ray-dim isolated NSs, or central compact objects (CCOs) in SNRs – such measurements are much more challenging. While all three of these NS classes are characterized by primarily thermal X-ray emission, CCOs are particularly noteworthy for being exclusively young (≲10⁴ yr) and nonvariable X-ray sources, located close to the center of SNRs. Unlike radio pulsars, they show no signs of interaction with their surroundings, such as extended nebular or jet-like emission. Furthermore, their spectra can be described entirely using blackbody or NS atmosphere models with luminosities ≳10³³ erg s⁻¹, without any unambiguous evidence for nonthermal magnetospheric emission (see Becker 2009; Gotthelf et al. 2013; De Luca 2017). Despite intensive searches (e.g., Mignani et al. 2008, 2019), no point-like or diffuse counterparts to any CCO at radio, infrared, or optical wavelengths have been found with present-day instrumentation. This places tight upper limits on the presence of potential binary companions and indicates no need for additional components to explain the observed spectral energy distribution.

We give an overview¹ of known CCOs in Table 1. Out of the ten known CCOs, X-ray pulsations have only been detected for those in Kes 79, Puppis A, and PKS 1209−51/52. These pulsations at 100 − 500  ms are likely attributable to the rotational modulation of the emission from hotspots on the NS surface (Gotthelf et al. 2010, 2005; Zavlin et al. 2000). Analysis of their spin-down points toward CCOs having a very weak magnetic dipole field. This, in combination with their steady emitting behavior, led to their designation as “anti-magnetars” (Gotthelf et al. 2013; Halpern & Gotthelf 2010a).

An important issue is that, observationally, few orphaned CCOs, meaning older pulsars without an apparent nearby SNR but in a similar location of (P, Ṗ)–space, have been found (Luo et al. 2015). It has therefore been suggested that the phenomenology of CCOs is caused by magnetic field burial due to the fallback accretion of material shortly after the supernova explosion. The magnetic field is then believed to resurface and evolve into a macroscopic dipolar field on timescales of ∼10⁴ yr, placing the object in a different region of parameter space (e.g., Bogdanov 2014; Luo et al. 2015; Gourgouliatos et al. 2020).

Since they are exclusively detected in X-rays, the proper motion of CCOs can currently only be measured directly using data from the Chandra X-ray Observatory, which presently is the only X-ray telescope providing sufficient spatial resolution for the task. Prior to this work, this has been attempted for five CCOs, with very diverse results: Most prominently, such measurements have been performed multiple times for RX J0822−4300, located in the SNR Puppis A (Hui & Becker 2006; Winkler & Petre 2007; Becker et al. 2012; Gotthelf et al. 2013; Mayer et al. 2020). Owing to the increasing temporal baseline, their precision has improved strongly over the years, the most recent value by Mayer et al. (2020) being (80.4 ± 7.7) mas yr⁻¹, at 68% confidence. This corresponds to a projected velocity of (763 ± 73) km s⁻¹ for a distance of 2 kpc, which potentially places this CCO at the fast end of the NS velocity distribution². In contrast, for the CCO in PKS 1209−51/52, Halpern & Gotthelf (2015) measure a much smaller proper motion of (15 ± 7) mas yr⁻¹, which corresponds to a projected velocity below 180 km s⁻¹ for a distance of 2 kpc. Only weak constraints can be placed on the proper motions of the CCOs in Cas A (DeLaney & Satterfield 2013) and G266.2−1.2 (Vela Jr., Mignani et al. 2019), owing to a paucity of X-ray bright field sources for astrometric frame registration. Finally, very recently, Borkowski et al. (2020) have shown that the CCO in G350.1−0.3 seems to exhibit only modest proper motion, contrary to expectations (see Sect. 4.6).

Alternatively to a direct measurement, it can be feasible to infer the proper motion of an NS indirectly, if its origin (i.e., the explosion site) and the age of the SNR are known sufficiently precisely (e.g., Lovchinsky et al. 2011; Fesen et al. 2006). However, if the estimate for the explosion site is solely based on the apparent present-day morphology of the SNR, this method is much more prone to systematic biases than the direct one. The main reason is that the apparent geometric center of an SNR does not necessarily coincide with its true explosion site. This was first shown by Dohm-Palmer & Jones (1996), who demonstrated numerically that SNRs expanding into interstellar medium with a density gradient can show significant offsets between their morphological center and their true explosion site, despite maintaining a relatively circular shape.

A previous study by Holland-Ashford et al. (2017) explores the relationship between the proper motions of 18 young NSs (four of which are CCOs), and the morphology of their host SNRs. Several of their proper motion measurements are obtained directly from X-ray imaging data, with others (including those for all four CCOs) being inferred indirectly. They do not find any correlation between the magnitude of SNR asymmetry and the projected velocity of the NSs. However, they find that an anticorrelation exists between the direction of proper motion and the direction of brightest ejecta emission, a proxy to highest ejecta mass. A similar study based on six systems by Katsuda et al. (2018) confirms this observation. In addition, the latter work does find a positive correlation between the degree of SNR asymmetry and the NS kick velocity. These findings support a hydrodynamic kick mechanism, where the NS reaches its extreme velocity due to asymmetric ejection of matter during the explosion, rather than via a neutrino-induced channel.

The most immediate handle on the age of an SNR is the direct measurement of its expansion. This approach yields a robust upper limit on the age since any physically justifiable assumption about the expansion history of the SNR would result in an age estimate lower than that for free expansion³. Tracing the expansion of an SNR is in principle possible at all wavelengths that allow to resolve the relevant emission features at sufficiently high resolution, such as the radio, optical, and X-ray regimes. Typically, the optical regime yields the most precise age constraints since astrometric uncertainties are very small. Moreover, the optical emission usually originates from high density ejecta clumps, which are subject to little deceleration by the circumstellar medium. A popular example is Cas A, whose explosion date was quite precisely determined by Thorstensen et al. (2001) and Fesen et al. (2006) to around the year 1670.

However, many young SNRs are not visible at optical wavelengths, often due to Galactic absorption. In this case, X-ray expansion measurements are a powerful tool to constrain age and internal kinematics, despite the lower spatial resolution and higher deceleration of the observed features. Examples for the viability of this method include the confirmation of the nature of G1.9+0.3 as the youngest Galactic SNR (Borkowski et al. 2014), the analysis of the propagation of the northwestern rim of the “Vela Jr” SNR (G266.1−1.2, Allen et al. 2015), or multiple recent efforts to trace the shock wave propagation of the TeV SNR RX J1713.7−3946 (Acero et al. 2017; Tsuji & Uchiyama 2016; Tanaka et al. 2020). Measurements of precisely dated very young SNRs (such as Cas A, Tycho’s, or Kepler’s SNR) show that expansion as traced in X-rays typically appears somewhat decelerated compared to the expectation for free expansion (Patnaude & Fesen 2009; Hughes 2000; Vink 2008). This highlights that SNR expansion measurements generally provide only an upper limit on their true age.

Central compact objects together with their parent SNRs present an ideal target for combining kinematic information from both the NS and the supernova shock wave. Since such systems are exclusively young, one can expect to be able to measure large expansion velocities, corresponding to a significant fraction of their undecelerated value. In addition, a precisely measured proper motion of the CCO provides accurate constraints on the supernova explosion site, the exact origin of the SNR shock wave. Combining these constraints promises unbiased results on the SNR’s kinematics and age as both the starting point and the present-day location of the shock wave are known.

The goal of this work is to establish a complete and internally consistent sample of proper motion measurements for CCOs and of expansion measurements for their host SNRs. We aim to provide direct constraints on the motion of all currently known CCOs (with appropriate data available) for which this task has not been previously performed. The ultimate target of this effort is to search for indications of violent supernova kicks, manifesting themselves in large velocities, and to constrain the location of the true center of targeted SNRs, with respect to their present-day morphology. For the three SNRs for which an expansion measurement has not yet been carried out at the time of writing, we perform an independent expansion study using our combined data set.

Our paper is organized as follows: We outline our data selection and reduction in Sect. 2 and describe our analysis strategy in Sect. 3. We then give a detailed overview over the results for each individual CCO/SNR in Sect. 4. Finally, we discuss physical consequences of the compiled sample of proper motion and expansion measurements in Sect. 5, and summarize our results in Sect. 6. Unless stated otherwise, all errors given in this paper are computed at 68%, and all upper limits at 90% confidence.

2. Data selection and reduction

Our initial list of potential targets for analysis consisted of the SNRs listed in Table 1. From this list, we excluded Puppis A, PKS 1209−51/52 (due to already precisely measured proper motion of the CCO) and Vela Jr. (due to the established lack of astrometric calibrators, Mignani et al. 2019). Even though the proper motion of the CCO in Cas A has already been constrained directly (DeLaney & Satterfield 2013), we decided to include it in our target list, in order to apply our own methodology to the data.

We searched the Chandra Data Archive⁴ for existing imaging observations covering the remaining potential targets. We took into account data sets from both Chandra imaging detectors, the High Resolution Camera (HRC, Murray et al. 2000) and the Advanced CCD Imaging Spectrometer (ACIS, Garmire et al. 2003). Whenever possible, we preferred all our observations to be taken with the same detector. The archival data were required to span a minimum temporal baseline of five years in order to qualify for a proper motion measurement. Suitable sets of observations that fulfill these criteria were found for six targets, the SNRs G15.9+0.2, Kes 79 (G33.6+0.1), Cas A, G330.2+1.0, RX J1713.7−3946 (G347.3−0.5) and G350.1−0.3. For G353.6−0.7, we found only a single suitable observation in imaging mode (taken in 2008), not permitting further analysis for our purpose. Out of the six target SNRs, the expansion of three (Cas A, G330.2+1.0, RX J1713.7−3946) has been measured directly in past studies (Fesen et al. 2006; Borkowski et al. 2018; Acero et al. 2017; Tsuji & Uchiyama 2016; Tanaka et al. 2020)⁵. This left us with the SNRs G15.9+0.2, Kes 79 and G350.1−0.3 as targets for an expansion measurement.

A journal of all observations used in our analysis is given in Table A.1. For some of our targets, such as G330.2+1.0 and G350.1−0.3, there are multiple observations with very deep coverage, promising high-precision results. In contrast, for RX J1713.7−3946, we found only two shallow observations that contain the CCO, taken with different detectors. This made performing a meaningful measurement more challenging and prone to systematic effects. In all cases treated in our sample, the target is only covered at two epochs that are spaced far enough apart in time to perform a proper motion measurement, despite most having multiple closely spaced observations at one of the epochs.

We reprocessed all archival data using the Chandra Interactive Analysis of Observations (CIAO, Fruscione et al. 2006) task chandra_repro with standard settings, to create level 2 event lists calibrated according to current standards. For ACIS data, we enabled the option to use the Energy-Dependent Subpixel Event Repositioning (EDSER) algorithm (Li et al. 2004), to be able to exploit on-axis data at optimal spatial resolution. We used CIAO version 4.12 and the calibration database (CALDB) version 4.9.0 throughout the work described in this paper. After the reprocessing step, we screened all data for soft proton flares, excluding the few affected time intervals with background levels 3σ above the respective quiescent average, yielding cleaned and calibrated event lists on which we based our subsequent analysis.

3. Methods

Our main analysis strategy consisted of several steps: First, we searched for serendipitous field sources for astrometric calibration (Sect. 3.1), which we then used to align the coordinate frames of the observations and simultaneously fit for the CCO’s proper motion (Sect. 3.2). For the three target SNRs of our expansion study, we exploited the astrometrically aligned data set to quantify the radial motion of key features of the SNR (Sect. 3.3).

3.1. Astrometric calibration sources

While the underlying principle of a proper motion analysis is rather straightforward, practically performing such a measurement in X-rays poses some challenges: Owing to uncertainties in aspect reconstruction, absolute source positions in Chandra images are only accurate to around 0.8″ at 90% confidence⁶. In order to maximize the precision of the CCO’s position and proper motion, it is thus necessary to align the coordinate system of our observations to an absolute reference, the International Celestial Reference System (ICRS). For this purpose, it is optimal to use serendipitous sources in the field of view with precisely known astrometric positions. Since most of our systems had not been studied in the context of a proper motion analysis previously, we searched for the presence of such calibrator sources systematically:

From the cleaned event lists, we extracted images over the entire detector area, which we binned at 1 × 1 and 2 × 2 detector pixels for ACIS and HRC, respectively. For ACIS, we included only events in the standard broad band, 0.5 − 7.0 keV. We then employed the CIAO task wavdetect to perform wavelet-based source detection and rough astrometric localization at a detection threshold of 10⁻⁵. We chose this rather “generous” threshold since we required the presence of a source in all observations in order for it to qualify as a calibrator. This eliminated the need to strictly reject potentially spurious sources in the individual detection runs.

For each of the observations, we matched the output source detection lists to the Gaia DR2 catalog⁷ (Gaia Collaboration 2016, 2018), which offers by far the best accuracy of all astrometric catalogs available at the time of writing. Those Gaia sources with an X-ray counterpart within at most 1.5″ in each observation made up our candidate list. We then visually inspected these candidates, excluding those sources that we deemed unreliable, for example because of large positional errors in X-rays or doubts in their correct optical identification. If many potential reference sources were found across the detector, we selected a subset of reliably identified and bright sources, spread evenly across the detector.

In further analysis steps, we occasionally decided to exclude objects from our list of frame registration sources. Reasons for this include the presence of a nearby overlapping source, the disagreement of positional offsets determined from different calibrator sources, or simply the failure to provide useful constraints compared to brighter sources. We display the calibrator source positions as well as the location of the respective target CCOs within their SNRs in Fig. 1. In addition, the list of astrometric reference sources used for our analysis is given in Table A.2, where we also indicate which of the sources we decided to exclude in subsequent steps.

Fig. 1. Overview over all target SNRs, with positions of the NS and the calibrator sources (labeled as in Table A.2) indicated on exposure-corrected images of the respective SNR. We show the deepest available observation for all targets except for RX J1713.7−3946, where we used the ACIS observation (ID 5559) instead of the HRC observation, due to its lower intrinsic noise level.

3.2. Proper motion analysis of CCOs

Our method for measuring the proper motion of CCOs followed the approach used in Mayer et al. (2020) and Becker et al. (2012). The main difference to these previous works was the performance of a simultaneous fit for the object’s proper motion and the astrometric calibration of individual observations, which provided a natural way to deal with observations of varying depth (Sect. 3.2.2).

3.2.1. PSF modelling and fitting

The first step was the precise measurement of the positions of all sources, that is, of the CCO and calibrator objects, at each epoch. We closely followed the procedure outlined and discussed in depth in Mayer et al. (2020) to obtain models of the Chandra point spread function (PSF) via ray-tracing simulations, using the ChaRT (Chandra Ray Tracer) online tool⁸ and the MARX software package⁹ (Davis et al. 2012, version 5.4.0). As in Mayer et al. (2020), we fitted the respective PSF models to source image cutouts, using the spectral and image fitting program Sherpa (Freeman et al. 2001)¹⁰. For each source in each observation, we obtained a profile of the logarithmic likelihood ℒ_ij(x, y) of the source position on a finely spaced grid (typical spacing ≲20 mas) in the region around the best fit. In our notation, the variables x and y describe relative positions as measured by Chandra in a tangent plane coordinate system around the approximate location of the CCO, while we use the indices i and j to label the individual observations and objects (including j = 0 for the CCO).

3.2.2. Coordinate transformation and proper motion fit

In the next step, we corrected the astrometry of each observation for minor misalignment with the ICRS by applying a linear transformation to the Chandra coordinate system. This transformation allowed for a positional shift as well as a small stretch and rotation of the coordinate system. Prior to the analysis, two groups of “interesting” parameters were unknown¹¹:

First, the astrometric solution of the CCO with the set of free parameters {μ_α, μ_δ, α₀, δ₀}¹², describing its proper motion and its position at a reference epoch t₀, in right ascension (α) and declination (δ), respectively. Second, a set of transformation parameters {Δx_i, Δy_i, r_i, θ_i} per observation i, describing translation in the tangent plane, as well as minor scaling and rotation corrections to the coordinate frame.

We made the realistic assumption that the errors on Gaia positions of calibrator stars can be neglected when compared to the errors in measured X-ray source positions. Thus, the true¹³ calibrator source positions were effectively known prior to the observation. Therefore, for each of our targets, there was a set Λ = {μ_α, μ_δ, α₀, δ₀, Δx₁, Δy₁, …, r_N, θ_N} of 4 + 4 × N parameters (with N the number of observations included), which were a priori unknown.

In order to constrain their values, we calculated the logarithmic likelihood for a given set of parameters, ℒ_tot(Λ), in the following manner: For a given astrometric solution {μ_α, μ_δ, α₀, δ₀}, we computed the true position of the CCO , at the epoch of observation i according to:

(1)

where t_i is the time of observation i. The fixed parameter t₀ describes the reference epoch, for which we always chose the time of the most recent observation in our respective data set. The true positions of the calibrator sources were computed analogously, using their known Gaia astrometric solutions.

These true positions were then converted to measured positions (x_ij, y_ij), to account for a coordinate system shift modeled by the parameters {Δx_i, Δy_i, r_i, θ_i}:

(2)

The likelihood of a given measured position was evaluated by interpolating the logarithmic likelihood grid ℒ_ij(x, y), which we computed as described in the previous section. By summing over the contributions of all N observations and S sources, the total likelihood for a set of parameters Λ could then be straightforwardly calculated:

(3)

In order to extract the most likely proper motion solution, it is necessary to derive the posterior probability distribution for Λ, which is implied by our likelihood function in its (4 + 4 × N)–dimensional parameter space. To achieve this in practice, we used the UltraNest¹⁴ software (Buchner 2021), which employs the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016, 2019), to robustly tackle this multi-dimensional problem. We used uninformative (flat, truncated) priors on all variables since we wanted our study to be as independent as possible from external assumptions.

This method is equivalent to performing the PSF fits to all sources at all epochs simultaneously, under the constraints of NS motion at constant velocity and known positions of all calibrators. The advantage compared to a truly simultaneous implementation is that our method allows for very efficient evaluation of the likelihood function via interpolation on a precomputed grid. With our approach, sources with large measured positional uncertainties only add little weight to the fit, and thus do not “wash out” the final result. Furthermore, covariances between parameters as well as non-Gaussian errors are properly taken into account.

However, the main caveat of this approach is that it neglects the possible presence of any systematic offsets beyond the linear coordinate transformation allowed by our fit. Such offsets could be caused by subtle differences between different detectors and roll angles or by the possible misidentification of a frame registration source with the wrong astrometric counterpart. To identify such potential systematic effects, we visualized our model’s predictions in the following way: For each point in the posterior sample of our fit, we computed the expected location of each source in the Chandra coordinate system (x_ij, y_ij). The distribution of these locations was then compared to the measured uncertainty contours, or equivalently the contours of ℒ_ij(x, y). This allowed us to check whether the results given by our simultaneous fit were sensible for all objects, and to recognize potentially misidentified astrometric calibrators. An example of such a comparison, illustrating the strongly varying scales of our astrometric uncertainties, is shown in Fig. 2.

Fig. 2. Illustration of PSF-fit likelihood for the source position and the prediction by our simultaneous proper motion fit for the CCO of G15.9+0.2. Each column shows data from the individual epochs (with the observation ID displayed in the title), while each row represents one particular object. The one to five sigma contours of the PSF-fit likelihood ℒij(x, y) are shown in green, and the predicted positions obtained from our simultaneous fit (see posterior distribution in Fig. 4) are indicated using a grayscale two-dimensional histogram. In the lower left corner of each panel, we provide a reference scale of length 0.3″. We display the true source position for our calibrator stars (known from Gaia) with a blue star. The uncertainty contours for observation 16766 are by far the smallest since the corresponding exposure was much deeper than at the other epochs (at ∼90 ks).

3.2.3. Correction for effects of solar motion and Galactic rotation

Ultimately, we attempt to extract physically meaningful quantities, such as NS velocity, with reference to the “rest frame” of the SNR. Therefore, it is necessary to take into account the expected motion of the NS from Galactic rotation alone (see e.g., Hobbs et al. 2005; Halpern & Gotthelf 2015). If this effect were not considered, our estimate of the physical CCO velocity would be contaminated by the motion of the sun and the rotation of the Galactic disk.

To determine the local standard of rest (LSR) of a system, we assumed the linear Galactic rotation curve and solar velocity given in Mróz et al. (2019). For each system, the expected velocity for corotation with the disk was expressed as two components of proper motion in equatorial coordinates, . From its distance d, proper motion, and LSR, it was then possible to determine the peculiar (projected) velocity of an object, as

(4)

One complication here is that, for some of our targets, the SNR distance d is not constrained precisely. For these cases, we assumed a flat probability distribution over the interval between its lower and upper limit (see Table 1 for available constraints). The resulting distribution of was then appropriately convolved with our distribution of the proper motion of the CCO. In practice, we found that the uncertainty on the LSR is significantly smaller than the uncertainty on the proper motion of all our targets.

A further technical issue is that, due to our choice of flat priors on (μ_α, μ_δ) in Sect. 3.2.2, the resulting probability distribution for is biased toward large values. This is because, in the transformation to , one essentially integrates over concentric circles of increasing radius, introducing an effective prior , and always leading to a peak of the distribution at nonzero . Therefore, in order to not overestimate the projected velocity of our CCOs, we reweighted the posterior distributions by a factor . The effect of this correction decreases with increasing significance of nonzero proper motion.

3.3. Expansion measurement of SNRs

When studying the dynamics of a young core-collapse SNR, a complementary tool to the analysis of the NS’s proper motion is the measurement of the SNR’s expansion. It provides the most immediate constraint on the shock velocity in its shell, and, by backward extrapolation, can be used to provide a very reliable constraint on its maximum age. Therefore, we took advantage of our results from the previous section, and measured the expansion of three CCO-hosting SNRs (G15.9+0.2, Kes 79, G350.1−0.3) in the following way:

First, we used the optimal transformation parameters {Δx_i, Δy_i, r_i, θ_i} for each observation i from Sect. 3.2.2, to perform the astrometric calibration of our data set. Since we fitted for the motion of the CCO and the frame calibration simultaneously, this had the advantage that, for observations taken close to each other in time, the position of the CCO (which is often the brightest source) was forced to be identical. This practically provided an additional reliable source for coordinate frame calibration. We aligned the data from all observations to a common coordinate system, using the CIAO task wcs_update with a transformfile customized to contain the optimal transformation parameters derived in Sect. 3.2.2. All data from observations taken at the same epoch were then merged and reprojected to a common tangent plane using reproject_obs. The systematic error in the coordinate frame alignment can be conservatively estimated from the combined statistical uncertainties of {Δx_i, Δy_i}, yielding typical values between 0.1″ and 0.3″.

In order to define the features of interest for our measurement, we extracted exposure-corrected images of each SNR in the standard ACIS broad band (0.5 − 7.0 keV). We then used SAOImage ds9 (Joye et al. 2003) to interactively define rectangular regions containing emission features that we found suitable to trace the expansion of the SNR, for example sharp filaments or individual clumps. Since we were mostly interested in measuring radial motion, the orientation of these boxes was adapted to what we deemed the probable direction of motion given the shape of the feature and its location within the SNR¹⁵. In order to avoid any possible confirmation bias toward the detection of expansion, we refrained from excluding regions in retrospect, even if their motion turned out to be almost unconstrained by the available data.

At this point, it would be possible to determine the motion of an individual feature and its errors, for instance by fitting a model for the flux profile to the data of both epochs (e.g., Xi et al. 2019). However, since some clumps and filaments targeted in this work show a rather complex profile, it seems infeasible to empirically model their shape, because the required number of free parameters would be high, or one would need to hand-tune the model for each feature. Alternatively, one could directly compare the observed emission profiles via a custom likelihood function (often some variant of χ²) quantifying their deviation, in one (e.g., Tanaka et al. 2020; Tsuji & Uchiyama 2016), or two dimensions (e.g., Borkowski et al. 2020). However, after some initial testing, we found that such a likelihood function behaves quite erratically in many cases, meaning that small sub-pixel variations in the shift of the profiles lead to large “jumps” in the likelihood. This makes a direct determination of errors prone to systematics and irreproducible for varying bin sizes.

Thus, in order to avoid biased results or underestimated errors, we used a bootstrap resampling technique (Efron 1982) to estimate uncertainties for the motion of all our features: We extracted the set of X-ray event coordinates x_i (with i = 1, 2 labeling the early/late observational epoch) along the suspected direction of expansion, as well as the corresponding exposure profile, within the rectangular region of each feature. For both epochs, we created N = 1000 randomized realizations of the observed data set , via sampling with replacement from the list of observed events in the region. Thereby, the statistical uncertainty of the emission profile was approximated by the intrinsic scatter of our resampled data sets.

From each resampled data set , we extracted one-dimensional flux profiles F_i(x) describing the emission along the (expected) direction of motion. For each bin x, we estimated a simple Gaussian flux error as , where C_i(x) describes the number of events in the bin. In order to estimate the optimal shift between the epochs, we defined the function f, which quantifies the deviation between the emission profiles at early and late epochs, in dependence of a positional shift ξ and an amplitude scale factor s:

(5)

where s was introduced to account for possible flux changes of the feature.

By minimizing f for each of the N realizations, we obtained a sample approximately distributed according to the likelihood of ξ. From this distribution, we extracted the angular speed μ_exp of the emission feature given by μ_exp = ξ/Δt, with Δt being the difference between the average event times of the early and late epochs. This is related to the projected expansion velocity via v_exp = μ_expd, where d corresponds to the (assumed) distance to the SNR.

The last step of estimating the global minimum of f was however complicated by the noisy nature of our function on small spatial scales. The reason for this is that the binning of events required to obtain an image is only an approximation to the underlying continuous flux distribution of a feature, which introduces a “granularity” in f. We resolved this issue by applying a Gaussian kernel to smooth the flux histograms F_i(x) before computing f. The width of the kernel was chosen so as to not over-smooth the prominent features in each region, but at the same time eliminate most of the purely statistical bin-to-bin fluctuations. Other approaches, like smoothing only one of the profiles or applying Poisson statistics to directly compare the two unsmoothed profiles, were found to typically yield similar or larger intrinsic errors. The basic principle of our approach is visualized for an example region in G350.1−0.3 in Fig. 3.

Fig. 3. Example expansion measurement in region D of G350.1−0.3 (see Sect. 4.6). Top and middle panels: exposure corrected images at the early and late epoch, and the white rectangle marks the measurement region. Bottom panel: one-dimensional flux profiles Fi(x) at the early (red) and late (blue) epoch. The markers display the binned distributions with errors, and the solid lines reflect the flux profiles smoothed by a Gaussian kernel. The dashed red line corresponds to the profile of the early epoch, shifted and rescaled by the indicated parameters (see upper right corner), with the shaded area displaying the associated uncertainty from bootstrapping. We define the x coordinate such that it increases in radial direction, meaning from west to east in this case.

The main advantage of our method is that we only extracted the optimal value of ξ for each realization. We therefore did not need to perform any “stepping” through parameter space to obtain the errors from a given likelihood threshold as we would consider this likely to yield underestimated statistical errors. Furthermore, our method allowed us to implicitly account for statistical fluctuations in both measurement epochs, which would be neglected when treating the deeper observation as a definitive “model” without intrinsic errors.

As a simple practical test of our approach, we used two late-time data sets of G350.1−0.3 (ObsIDs 21118 and 21119) and applied our method, as if trying to measure expansion. Naturally, the true shift of features between these two observations is essentially nonexistent due to the small time difference of around one day. Therefore, one expects the results of our experiment to cluster around zero, tracing not the motion of features but the intrinsic noise of our method. We find that the overall deviation from zero of our measurements can be satisfactorily accounted for with statistical and systematic errors. More quantitatively, the sum of normalized squared deviations from zero, ,i s found to be χ² = 19.1 for K = 18 regions (i.e., “degrees of freedom”), when accounting only for statistical uncertainties. Adding in quadrature the estimated systematic error of 0.1″, we obtain χ² = 14.3. The expected 68% central interval of a χ² distribution with 18 degrees of freedom is (12.1, 23.9). Therefore, our test demonstrates that we obtain reasonable estimates for our total uncertainty, and that our overall error estimates are likely rather conservative.

4. Results

4.1. G15.9+0.2

G15.9+0.2 is a small-diameter Galactic SNR, first detected in the radio by Clark et al. (1973). In X-rays, its morphology resembles that of an incomplete, but very symmetric, circular shell. It shows a northwestern “blowout” region with only weak apparent emission, likely due to reduced density of the circumstellar medium (Sasaki et al. 2018). Its distance is only constrained to a relatively uncertain range of values (7 − 16 kpc, Tian et al. 2019), and accordingly, the age shows correlated uncertainties of a factor of two (2900 − 5700 years, Sasaki et al. 2018).

The bright CCO CXOU J181852.0−150213 is quite clearly offset by around 35″ from the apparent geometrical center of the SNR toward the west or southwest. As first noted by Reynolds et al. (2006), this seems to imply a quite large transverse velocity, or a proper motion of around 10 mas yr⁻¹ for an assumed age of 3500 years.

Unfortunately, the available data at the early epoch are split into three nonconsecutive observations of 5.0, 9.2, and 15.1 ks length, with the detector aimpoint being located around 3′ north of the CCO, which reduces the usability of this data set for our purpose. In contrast, the late observation is very deep (with an exposure around 90 ks) and aimed directly at the CCO, resulting in very small statistical errors on its late-time position.

The resulting posterior distribution from our astrometric fit can be seen in Fig. 4. The 68% central credible intervals for the proper motion are given by (−17 ± 12, −4 ± 10) mas yr⁻¹, in right ascension and declination, respectively. This seems to agree well with the expected proper motion direction given the CCO’s present-day location in the SNR.

Fig. 4. Corner plot displaying the joint and marginal posterior distributions for the parameters of the astrometric solution for the CCO in G15.9+0.2. The astrometric calibration parameters {Δxi, Δyi, ri, θi} have been marginalized over, and are not shown here in order to keep the figure uncluttered (see Table B.1 for constraints on all parameters). In the two-dimensional correlation plots, we indicate the 1σ, 2σ, and 3σ contours, meaning the smallest regions containing 39.3%, 86.5%, and 98.9% of the total probability mass, with solid lines. In the marginalized one-dimensional histograms, we indicate the median and 68% central interval of the posterior probability distribution, with dashed lines. The units of proper motion and the astrometric zero-point are mas yr−1 and arcsec, respectively. This plot and analogous figures have been created using the corner.py package (Foreman-Mackey 2016).

While this cannot be counted as an unambiguous detection of significant proper motion for the CCO, a proper motion of the order of 15 mas yr⁻¹ toward west seems realistic given the SNR geometry. It would imply a rather large projected velocity of the NS around 700 km s⁻¹ when scaling to a distance of 10 kpc. The formal 90% upper limit of 25 mas yr⁻¹ (corresponding to 1200 km s⁻¹ for that distance) on the total proper motion with respect to the LSR of G15.9+0.2 is physically unconstraining.

The almost perfectly circular shape of G15.9+0.2 allowed us to define quite intuitive regions along its outer rim for an expansion measurement. They are illustrated in Fig. 5. As can be seen, all regions trace the motion of the edge of the SNR shell, and due to the observed morphology, we would expect to observe quite similar expansion speeds along all directions. Minor problems faced in their definition were the bright star “1” superimposed on the eastern shell and a detector chip gap crossing the southernmost part of the remnant, both of which we avoided with our choice of regions.

Fig. 5. SNR G15.9+0.2 with regions used for the expansion measurement indicated. We show a smoothed, exposure-corrected image of the SNR with logarithmic scaling of the color map. We overlay the boxes from which we extracted the one-dimensional emission profiles for our expansion measurement, and label them alphabetically for easier identification and to avoid confusion with numeric labels for the astrometric calibrator stars.

In order to define a rough reference point for the expansion of the SNR, we used SAOImage ds9 (Joye et al. 2003) to approximate the shell of the SNR with a circle. Its center at corresponds to the probable location of the explosion, with an assumed 1σ–error of 15″, corresponding to around 10% of the SNR’s radius. This is an intuitive way to estimate the origin of the SNR’s expansion since the CCO is clearly offset with respect to the apparent center, and its proper motion is not constrained very precisely at this point. This reference point allowed us to obtain a rough estimate of the angular distance ϑ traveled by the shock wave at the location of each region.

The results of our expansion analysis are displayed in Fig. 6 (for a full illustration of the data at both epochs, see Fig. C.1). Due to the comparatively low exposure time available in the early epoch, the statistical error bars on the proper motion of the SNR shell are quite large. We estimate that the systematic error introduced by astrometric misalignment of the observations is at most 0.15″ (see Table B.1), corresponding to a proper motion error around 15 mas yr⁻¹. Therefore, the overall effect of this systematic error is rather small, also since there does not seem to be any bias toward detecting motion along a certain direction in our results.

Fig. 6. Constraints on the expansion of G15.9+0.2. Top: for each feature (labels as in Fig. 5), we plot the measured angular expansion speed μexp, and the corresponding projected velocity vexp for a distance of 10 kpc, against its angular distance to the SNR center ϑ. To assess their respective significance, we have used black error bars for those measurements that enclose > 95% of the total probability mass on either side of 0, and gray for the others. In addition, we indicate the expected free expansion speed for ages of 1000, 2000, and 3000 years with dashed lines. Bottom: expansion rate τ−1 as measured in each individual region. The thick black line with the gray-shaded area represents the combination of the individual measurements (or weighted average), assuming uniform expansion. All error bars are at 68% confidence. The data underlying this figure are given in Table B.2.

We obtained a quantitative estimate of significance and rate of expansion in the following way: For each region, we combined the measured expansion speed μ_exp with its distance from the SNR center ϑ to obtain an estimate of the expansion rate τ⁻¹ = μ_exp/ϑ. In order to combine the constraints from individual boxes, we made the simplest possible assumption, which is that the SNR is expanding uniformly in all directions. The large size of the measured error bars presently does not warrant more complex (albeit more realistic) models, such as an azimuthal variation of the expansion velocity. Our assumption allowed us to directly multiply the inferred likelihoods from all regions to estimate an error-weighted average expansion rate of G15.9+0.2. In this process, we marginalized over systematic uncertainties introduced by astrometric misalignment of observations and the unknown location of the SNR center. This means we recalculated μ_exp and ϑ many times for random samples drawn from possible center locations and systematic shifts, and summed up the combined distributions of τ⁻¹ afterward. Thereby, we removed the need to artificially convolve distributions.

As can be seen in the lower panel of Fig. 6, our combined best estimate for the current average rate of expansion is τ⁻¹ = (4.3 ± 1.9)  × 10⁻⁴ yr⁻¹. This finding is moderately significant at around 2.5σ, as the probability for the true value to be less than or equal to zero is found to be 0.6%. Given the almost perfectly circular shape of the SNR, and the measured error bar sizes, presently there appears to be no obvious deviation from uniform expansion. Thus, while it is infeasible to constrain the kinematics of individual portions of the shell at this point, the statistical distribution of the results from our seven regions supports nonzero expansion. During the refinement of our analysis strategy, we attempted several techniques to measure the motion of the rim, and also attempted a simpler form of astrometric alignment of the two epochs (applying only a constant two-dimensional shift of the coordinate system). For all these cases, we find that overall expansion still holds.

It may therefore be tempting to state that, assuming perfectly free expansion, our average expansion rate would correspond to an age τ ∼ 2500 yr. Of course however, given the available data and the size of statistical errors, this value constitutes only a zeroth-order estimate for the SNR’s age and is therefore to be taken with caution until the expansion of G15.9+0.2 has been independently measured.

4.2. Kes 79

Kes 79 (G33.6+0.1) is a bright mixed-morphology SNR, whose X-ray appearance is dominated by multiple filaments and shells, as well as the bright central X-ray source CXOU J185238.6+004020 (Sun et al. 2004). In the radio and mid-infrared regimes, its morphology appears similarly complex and filamentary (Giacani et al. 2009; Zhou et al. 2016). The distance of Kes 79 is not very well constrained. One of the most recent estimates, inferred from H I absorption, yields a distance of 3.5 kpc (Ranasinghe & Leahy 2018), whereas CO observations (Kuriki et al. 2018) suggest a distance of 5.5 kpc. The age of Kes 79 is estimated via detailed X-ray spectroscopy to lie between 4.4 and 6.7 kyr (Zhou et al. 2016).

Our measurement of the CCO’s proper motion (Fig. 7) is perfectly consistent with zero, the 68% central credible intervals being . This can be converted to a 90% upper limit on the CCO’s peculiar proper motion of around 19 mas yr⁻¹. Equivalently, its tangential velocity is constrained at < 450 km s⁻¹, when scaled to a distance of 5 kpc.

Fig. 7. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in Kes 79. See Table B.1 for constraints on all parameters.

While the limit on the CCO’s transverse velocity is not very constraining by itself, our measurement has a noteworthy implication on the NS’s timing properties, particularly its period derivative: the magnitude of the Shklovskii effect, an entirely kinematic contribution to the measured period derivative (Shklovskii 1970), can now be constrained to < 4 × 10⁻¹⁹ s s⁻¹ at 90% confidence (for a distance of 5 kpc). This corresponds to < 5% of the measured period derivative (Halpern & Gotthelf 2010a), which justifies to neglect this contribution to first order.

As can be seen in Fig. 8, Kes 79 exhibits several sharp filaments and shell fragments, which are in principle ideally suited for tracing the propagation of the supernova shock wave through the interstellar medium. However, our data set has several caveats: First, the exposure times of all observations are comparatively low (at 29.5, 9.8, and 10.0 ks), such that statistical fluctuations in the measured flux profiles are quite high. Second, for all three observations, the telescope pointing directions and roll angles were different. Therefore, since we attempted to avoid all the chip gaps of the ACIS-I detector, we had to exclude several emission features from our region definitions, which would otherwise have increased the measurable signal quite significantly (e.g., close to regions A, B, C, H). In order to quantify the distance between our emission features and the SNR center, we made use of our results in Sect. 5.2, establishing our estimated explosion site as (relatively uncertain) reference point.

Fig. 8. SNR Kes 79 with regions used for the expansion measurement indicated (as in Fig. 5).

The measured proper motion μ_exp of the shock wave in the individual regions is displayed in the upper panel of Fig. 9. Again, we observe that the overall picture is visually dominated by large statistical uncertainties, with most points scattered around zero, and several regions appearing almost unconstrained. However, the important exception is filament A as it, on its own, shows evidence for nonzero proper motion at ∼4σ statistical significance. In contrast, no regions show stronger deviations from zero than ∼1σ in the negative direction, which would be expected if this result was simply due to underestimated errors. We conservatively estimate the systematic astrometric error from coordinate system alignment to around 0.25″ (see Table B.1), corresponding to an error on the angular motion around 17 mas yr⁻¹. While the proximity of region A to a chip gap seems concerning for our overall interpretation, even an inspection by eye of the shock profiles in region A (see Fig. C.2) shows an obvious shift in the exposure-corrected flux profile, strengthening our confidence in the observed result.

Fig. 9. Constraints on the expansion of Kes 79 (as in Fig. 6). The assumed distance to Kes 79 for the conversion to vexp is 5.0 kpc. The data underlying this figure are given in Table B.2.

Analogously to Sect. 4.1, we combined the measurements of individual regions to constrain the uniform expansion rate of Kes 79 at . The probability for the true value to be less than or equal to zero is estimated to 1.6%. Thus, the nominal statistical significance of nonzero expansion is still low (at around 2σ). However, it should be noted that it is possible that certain features, traced for instance by the filament in region A, propagate more freely than others that could be more strongly decelerated. Therefore, the assumption of uniform expansion is likely unrealistic given the complex X-ray and radio morphology of Kes 79. Indirect signs of nonuniform expansion have actually been found in spatially resolved X-ray spectroscopy of Kes 79 by Zhou et al. (2016). If the expansion is indeed nonuniform, a combination of likelihoods is technically not valid, since one does not truly measure the same underlying quantity in different regions.

4.3. Cas A

Cas A (G111.7-2.1) is the youngest known Galactic core-collapse SNR at an age of around 350 years (e.g., Fesen et al. 2006; Alarie et al. 2014). Its X-ray morphology can be clearly separated into an inner, thermally emitting ejecta-dominated shell, and an outer forward shock region, which exhibits mainly synchrotron emission and expands at around 5000 km s⁻¹ (DeLaney et al. 2004; Patnaude & Fesen 2007, 2009). The distance to Cas A has been precisely constrained to (3.33 ± 0.10) kpc via an analysis of the expansion of filaments in the optical (Alarie et al. 2014).

The CCO CXOU J232327.9+584842 was discovered in the Chandra first-light image of Cas A (Tananbaum 1999). There exists a past proper motion analysis for the CCO performed by DeLaney & Satterfield (2013), which yields a final measurement of the transverse velocity of (390 ± 400) km s⁻¹. This is equivalent to a total proper motion of (24 ± 25) mas yr⁻¹ toward the southwest, for their assumed distance of 3.4 kpc. Nevertheless, since they used different statistical methods than we did, and did not employ any PSF modeling, we decided to reanalyze the data to produce a comparable measurement to the other objects in this paper.

Our initial search for detectable calibrator sources yielded only one source common to all five observations (labeled as “1” here), which is identical to the one found in DeLaney & Satterfield (2013). However, instead of using quasi-stationary flocculi of Cas A for image registration, we decided to stack the late-time observations to increase our chance to detect faint sources. To do this, we aligned the observations using point-like sources found by wavdetect that were highly significant in all observations (> 10σ). These correspond mostly to bright clumps of emission of the SNR, which can be regarded as effectively stationary on the relevant timescales (< 10 days). We applied the CIAO task wcs_match to find the optimal transformation between two respective epochs, and then reprojected and merged all late-time observations using the wcs_update and reproject_obs tasks.

Using this combined late-time data set, we detect two more point sources common to early and late observations, and with close positional match to Gaia stars (see Fig. 1 and Table A.2). Due to their location in the outer regions of Cas A (see Fig. 1), we made some further considerations to ensure their correct identification: We verified that their spectral nature is softer than that of the surrounding emission, using an archival 1 Ms ACIS observation (PI U. Huang, ObsIDs 4634−4639). Furthermore, we emphasize that if these sources were actually associated with ejecta clumps of the SNR rather than being foreground sources, they would be expected to show very large proper motion due to their location at the edge of the remnant shell. Such behavior would be clearly recognizable with respect to the secure source 1. Since we do not observe this, we conclude that all three astrometric calibrators are probably correctly identified.

The results of our proper motion measurement are displayed in Fig. 10. We find evidence for nonzero proper motion directed toward the southeast, with a best-fit value of . An analogous fit to the data from the late epochs taken individually yields a similar result, showing that the impact of our merging technique is rather small.

Fig. 10. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in Cas A. See Table B.1 for constraints on all parameters.

Converting our two-dimensional probability distribution for (μ_α, μ_δ) to a distribution for the absolute value of proper motion, we obtain a 68% central credible interval . This corresponds to a transverse physical velocity of (570 ± 260) km s⁻¹ at the distance to Cas A of 3.4 kpc. Here, including the effect of Galactic rotation does not considerably alter the result. Fesen et al. (2006) performed an analysis of the expansion of high-velocity optically emitting ejecta knots, yielding a precise estimate for the explosion date of Cas A (around the year 1670). Connecting the explosion site inferred in Thorstensen et al. (2001) with the present-day location of the CCO, they constrained its velocity to around 350 km s⁻¹ at a position angle 169° ± 8°. Our measurement is marginally consistent with both that indirect estimate and the analysis of DeLaney & Satterfield (2013).

4.4. G330.2+1.0

The SNR G330.2+1.0 is located at a distance of at least 4.9 kpc (McClure-Griffiths et al. 2001), and can be seen as a faint but complete shell in X-rays. Its emission is almost entirely nonthermal in nature, except for a region in the east showing thermal emission (Williams et al. 2018). Analysis of shock-velocity and filament expansion points toward a rather young age of the remnant, likely below 1000 years (Williams et al. 2018; Borkowski et al. 2018). As part of their expansion analysis, Borkowski et al. (2018) note that no evidence for “discernible motion” of the CCO CXOU J160103.1−513353 is found, which we attempted to further quantify here.

We can see the posterior distribution for the astrometric solution of the CCO in Fig. 11. We find a very small proper motion toward the south of . There are no obvious systematic discrepancies between the posterior predictions from our fit and the PSF fit contours, and we find no need to artificially add a systematic error component to the likelihood contours. Thus, due to the very long exposure times of the observations at early and late epochs, and the large number of clean astrometric calibrators, we find the small statistical error bars of our result plausible.

Fig. 11. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in G330.2+1.0. See Table B.1 for constraints on all parameters.

The expected proper motion of an object corotating with the Galactic disk at the location of G330.2+1.0 at a distance between 4.9 and 10.0 kpc is (−4.9 ± 0.8, −4.5 ± 0.6) mas yr⁻¹. Thus, the measured proper motion of the CCO is perfectly consistent with zero when effects of Galactic rotation are taken into account. We quote the corresponding 90% upper limit on the CCO’s peculiar proper motion, which lies at 9.9 mas yr⁻¹. This corresponds to a transverse kick velocity below 230 km s⁻¹, assuming a distance of 5 kpc.

This case demonstrates that, even given the “poor” spatial resolution of X-ray data when compared to the optical, effects of Galactic rotation cannot generally be neglected for proper motion measurements. Furthermore, our very tight constraints on its motion imply that the maximum angular distance traveled by the CCO during the “lifetime” of G330.2+1.0 (≲1000 years) is only around 10″. This is close to negligible compared to the SNR’s present-day radius around 5′, making CXOU J160103.1−513353 quite literally a “central” compact object.

4.5. RX J1713.7−3946

RX J1713.7−3946 (G347.3−0.5) is one of the brightest known emitters of ultra-high-energy gamma rays (e.g., Hess & Abdalla 2018). It appears quite luminous in X-rays and comparatively dim in the radio regime, with its X-ray spectrum being absolutely dominated by nonthermal synchrotron emission (Okuno et al. 2018). Morphologically, the SNR can be described as elliptically shaped, with the bright western shell appearing to consist of multiple ring-like structures, which in turn consist of thin filaments on smaller scales (Cassam-Chenaï et al. 2004). The distance to the SNR is estimated to be around 1.0 − 1.3 kpc, based on CO and H I observations of clouds interacting with the shock wave (Fukui et al. 2003; Cassam-Chenaï et al. 2004). The expansion of the outermost SNR filaments has been directly measured in three quadrants (Acero et al. 2017; Tsuji & Uchiyama 2016; Tanaka et al. 2020), yielding projected shock velocities up to 3900 km s⁻¹ and implying an age for RX J1713.7−3946 between 1500 and 2300 years.

The position of the CCO 1WGA J1713.4−3949 is slightly offset from the apparent center of the SNR (Pfeffermann et al. 1996) by around 4.5′ toward the southwest. Since the system is located quite nearby and the SNR age is likely rather low, it therefore seems reasonable to expect significant proper motion. For instance, assuming the geometric center to correspond to the explosion site, we would expect motion of around 130 mas yr⁻¹ for an age of 2000 years.

There are three archival Chandra observations of the CCO¹⁶, one from 2005 using ACIS-I, one from 2013 using HRC-I, and one from 2015 (ID 15967) using ACIS-I. The latter observation was carried out using only a subarray of the detector, and with an offset aimpoint. Therefore, it is not suitable for our analysis, leaving us with only two moderately deep observations.

Due to the intrinsically high noise level in the HRC observation, it proved very difficult to find adequate reference sources that could be convincingly cross-matched to a Gaia counterpart. Thus, in order to avoid having to rely on potentially erroneous optical associations, we decided to modify our approach, measuring relative source offsets between our two observations directly instead of measuring absolute positions. To achieve this, for each source, we convolved the PSF-fit likelihood contours of the late epoch with the “mirrored” contours of the early epoch (meaning we set x → −x, y → −y). This effectively subtracts the two positions from each other, resulting in a probability distribution for the astrometric offset between the two observations for a given object. To account for the unknown possible proper motion of our calibrators, we convolved their offset distributions with a Gaussian of width σ = 10 mas yr⁻¹ × Δt, where Δt = 7.9 yr corresponds to the time difference between the two epochs. We then performed our astrometric fit with only {μ_α, μ_δ} and a single set of frame registration parameters {Δx, Δy, r, θ} as free parameters. The advantage of this method is that it does not require precisely knowing the true source positions, and “only” assumes their motions to not be much larger than 10 mas yr⁻¹.

The resulting posterior distribution can be seen in Fig. 12. We measure a proper motion of , which is consistent with zero within its quite large statistical uncertainties. We verified this result by checking that for both of our calibrators taken as the sole astrometric reference (applying only a simple translation of the coordinate system), we also obtain proper motion consistent with zero. This demonstrates that our result is not overly biased by our method or choice of calibration sources. We estimate the 90% credible upper limit on the peculiar proper motion to be around 48 mas yr⁻¹, corresponding to a limit on the two-dimensional space velocity of around 230 km s⁻¹ for a distance of 1 kpc. A future observation would certainly allow these large error bars on the CCO’s proper motion to be greatly reduced. Given the issues we faced, it would probably be advantageous to use the ACIS detector for such a task as it is generally better suited for the detection of faint sources than the HRC.

Fig. 12. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in G347.3−0.5. For this target, a relative astrometric frame registration method was applied as described in the text. We therefore measured no absolute positions, but only relative source displacements, which correspond to the components of proper motion. See Table B.1 for constraints on all parameters.

4.6. G350.1−0.3

G350.1−0.3 exhibits a very unusual morphology at multiple wavelengths: In the radio, it shows a distorted and elongated shape (Gaensler et al. 2008), and in X-rays, a very bright irregularly shaped clump in the east dominates the overall morphology with only weak emission otherwise. The most likely explanation for the bright clumpy X-ray emission is the interaction of hot, metal-rich ejecta with a molecular cloud, based on which a distance of 4.5 kpc seems likely (Bitran et al. 1997; Gaensler et al. 2008). Lovchinsky et al. (2011) performed spectral analysis of a deep ACIS-S observation of the SNR and from that estimate an age of 600 − 1200 years for the SNR.

The bright X-ray source XMMU J172054.5−372652 is located 3′ west of the dominant emission region, quite offset from the apparent center of the remnant. Assuming the origin of the CCO to be located at the apparent center of emission of the SNR, Lovchinsky et al. (2011) predict a projected velocity of 1400 − 2600 km s⁻¹ for the NS. This would require an extremely strong kick to have acted on the NS during the supernova explosion, and the implied proper motion of 65 − 130 mas yr⁻¹ should be easily detectable for us.

During the preparation of this work, Borkowski et al. (2020) published a study on the expansion and age of G350.1−0.3¹⁷. They provide a detailed analysis of the motion of many individual emission clumps as well as a brief measurement of the CCO’s proper motion. Their main results include the measurement of rapid expansion of the SNR, constraining its age to be at most around 600 yr. Moreover, they detect comparatively slow motion of the CCO toward the north, at significant uncertainty. Using our independent methods, we aimed to confirm their measurement of the SNR’s expansion and obtain quantitative constraints on the proper motion of the CCO.

The total exposure of the data set at both early and late epochs is exquisite, at around 100 and 200 ks, respectively. The data are therefore ideal for obtaining precise constraints on the motion of both the CCO and the SNR ejecta in the nine-year time span between the two epochs.

From our proper-motion fit, where we treated all five late-time observations independently, we obtain the posterior distribution shown in Fig. 13. We find that the CCO appears to be moving northward (with large uncertainty), at . From this, we obtain a 68% central credible interval for the peculiar proper motion of , corresponding to a transverse velocity component of at a distance of 4.5 kpc. Borkowski et al. (2020), using positions from Gaussian fits to the CCO and field sources, find a proper motion of (−5, 14) mas yr⁻¹, which is consistent with our finding of slow motion approximately due north. Both measured values are certainly much lower in magnitude (and quite likely different in direction) than predicted from the SNR geometry in Lovchinsky et al. (2011). The implications of this will be discussed in Sect. 5.2.3.

Fig. 13. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution of the CCO in G350.1-0.3. See Table B.1 for constraints on all parameters.

As a first test for the expansion of G350.1−0.3, we performed a simple image subtraction of the early-time from the late-time image, to check for obvious changes between the two epochs. The lower panel of Fig. 14 displays the differential between the merged, exposure-corrected and smoothed images. There is very obvious evidence for substantial eastward motion of the bright ejecta clump, which appears more significant toward its outer edge. In addition, the image nicely confirms the measured proper motion of the CCO, with its position showing a small but visible offset between the two epochs. This image is a very useful qualitative confirmation of the SNR’s rapid expansion, independent of the exact method used in quantitative analysis.

Fig. 14. SNR G350.1−0.3 with regions used for the expansion measurement indicated. Top panel: as in Fig. 5. Bottom panel: difference image between the early and late epochs. We use a symmetric logarithmic color scale to highlight features on small scales, with blue corresponding to the late-epoch emission being brighter, and red being the opposite. For the main image, we smoothed the data with a Gaussian kernel of 1″, while in the inset we show the region of the bright clump smoothed with only 0.5″.

As can be seen in the upper panel of Fig. 14, the complex morphology and bright emission of the eastern clump allowed us to define numerous small-scale features of emission. Their propagation was measured in order to infer the internal kinematics and global expansion behavior of the ejecta. In contrast, for the much fainter emission in the north and south of the SNR, we defined regions of diffuse emission on much larger scales, which we hoped would allow us to trace the interaction of the shock wave with the interstellar medium.

For each feature, we used the indirectly inferred location of the explosion site from Sect. 5.2 to estimate the angular distance ϑ that has been traversed since the supernova. The measured expansion speeds μ_exp for all 18 regions can be seen in the upper panel of Fig. 15 (for a full illustration of the data at both epochs, see Fig. C.3). We estimate systematic uncertainties due to astrometric frame registration to be smaller than 0.1″, or equivalently 11 mas yr⁻¹. The proper motion of most features is much better constrained than for those in Sects. 4.1 and 4.2. Therefore, we combined our probability distributions for ϑ and μ_exp directly to obtain the free expansion age of the SNR, defined as τ = ϑ/μ_exp. The resulting constraints on τ are displayed in the lower panel of Fig. 15.

Fig. 15. Expansion and age estimate for G350.1−0.3. Top panel: as in Fig. 6, with an assumed distance of 4.5 kpc. Bottom: constraints on the expansion age τ from each individual region, with error bars at 68% confidence. The data underlying this figure are given in Table B.2.

As we can see, there is spectacular evidence for significant expansion in almost all regions investigated. The maximal observed expansion speed is around 250 mas yr⁻¹, corresponding to a projected shock wave velocity close to 6000 km s⁻¹ for a distance of 4.5 kpc. In combination with its angular distance from the SNR center, we find the lowest expansion age for region E, at τ = (580 ± 50) yr, including systematic errors. These findings are in great agreement with those from Borkowski et al. (2020), who measure the fastest expansion at a very similar location (labeled “B3” by them, finding τ = 590 yr) despite a different region shape. Among the “runners-up” in terms of likely expansion rate are the regions A, B, F, M, P, all of which show values consistent with τ ∼ 700 yr. However, none of these regions independently confirms that τ < 700 yr. In principle, the highest undecelerated expansion age τ can be considered a hard upper limit on the true age¹⁸. However, it seems important to ask why region E would indicate expansion at a higher rate than regions B or C. Its location does not appear to be “special” in the sense that it is not located close to the edge of the bright interacting ejecta clump. Naively, one would expect the least decelerated expansion there, if the bent shape of the clump was caused by the direct interaction of the shock wave with an obstacle. We thus choose to place a more conservative upper limit of 700 yr on the true age of G350.1−0.3, which nonetheless confirms it as one of the three youngest known Galactic core-collapse SNRs (see Borkowski et al. 2020).

The internal kinematics of the bright eastern clump are found to vary quite significantly with location: In particular, its southern elongated feature shows a significant velocity gradient, with the regions I and J being decelerated more strongly than the southern regions K and L. Similarly, the easternmost regions B and C display quite rapid motion when compared with the rest of the clump, especially the nearby region D. Furthermore, upon closer inspection of flux profiles at early and late epochs, it becomes evident that some features in the clump, for example regions J and K, have evolved in shape and decreased in brightness (see Fig. C.3) over the time span of nine years. The general trend seems to be that the X-ray emission in these regions becomes less “clumpy” and more diffuse as the SNR evolves. All these findings trace the interaction of the ejecta with an X-ray dark obstacle (Gaensler et al. 2008), which appears to cause the unique outward bent shape of the clump. Clearly, the interaction is the strongest close to the center of the clump, leading to the lowest measured expansion rates there, and higher velocities for its eastern and southern parts. In this scenario, if parts of the shock wave encounter a reduced ambient density after having passed by the obstacle, morphological changes of associated emission features could be a natural consequence.

Finally, we clearly detect the motion of the larger, more diffuse emission features in the north and south of the SNR (labeled M to R). The expansion rates we find here are a lot less certain, but on average appear comparable to or slightly smaller, than those in the bright emission clump. Here, there are some subtle differences between our findings and those of Borkowski et al. (2020), for instance, for the faint regions M and O (labeled “NNE” and “K” by them). In these regions, they infer rapid, almost undecelerated expansion, corresponding to expansion ages < 700 yr, with quite small statistical errors. While our measurement for region M is formally consistent with theirs, we find much larger statistical errors in both regions M and O, making our constraints on the local degree of deceleration of the shock wave rather weak. These differences in uncertainty could be caused by different shapes or orientations of the measurement regions. However, we consider the most likely origin to be the difference in analysis methods, with our one-dimensional resampling technique yielding more conservative errors than a direct comparison of two-dimensional flux images.

Overall, we do not find significant evidence for spatially varying degrees of deceleration within the faint diffuse emission regions. The fact that every single region of the SNR shows (more or less) significant signatures of expansion away from a common center clearly proves that the unusual morphology of the SNR is not caused by the superposition of two unrelated objects, as was conjectured by Gaensler et al. (2008). Instead, G350.1−0.3 is indeed a single and highly peculiar SNR.

5. Discussion

5.1. Proper motion and NS kinematics

This work constitutes the first step toward the ambitious goal of building a sample of consistently measured transverse velocities for all CCOs. At the present time, only four members of the class (the CCOs in Puppis A, Cas A, G350.1−0.3 and PKS 1209−51/52) have a measured nonzero proper motion. However, we hope that future observations will aid in reducing error bars and obtaining more precise constraints for those CCOs for which only upper limits on the proper motion could be established here.

We display an overview of the astrometric solutions of all known CCOs, derived in this work and previous studies, in Table 2. Our work almost doubles the number of CCOs with quantitative proper-motion measurements from five to nine. For several systems in this work, the results of our proper motion measurements are markedly different from reasonable expectations. For instance, for both G350.1−0.3 and RX J1713.7−3946, the present-day location of the CCO within the SNR suggests measurable nonzero proper motion. For RX J1713.7−3946, we are unable to obtain any significant signature of proper motion despite the small distance to the system. For G350.1−0.3, we do measure a mildly significant nonzero value, which is however considerably smaller than what the SNR morphology would suggest. This is similar to the case of PKS 1209−51/52, where Halpern & Gotthelf (2015) measure very small proper motion, despite a striking offset of the CCO from the apparent SNR center. These examples illustrate that estimates of the NS kick purely based on SNR geometry have to be interpreted with care, and in most cases cannot replace a direct measurement. As demonstrated by Dohm-Palmer & Jones (1996), an inhomogeneous circumstellar medium can easily result in an offset between the SNR’s explosion site and its apparent center, as the shock wave experiences different degrees of deceleration in different directions, leading to a distortion of the remnant’s morphology. This effect is investigated in detail for Tycho’s SNR by Williams et al. (2013), who show that the observed density and shock velocity variations along its outer rim imply a significant offset between the apparent and true SNR center, by up to 20% of its radius. Therefore, even if a central NS appears clearly offset from its host’s morphological center, such an offset need not be due to its proper motion, but may also be caused – solely or to some fraction – by the SNR’s distorted shape.

The two most important factors limiting the precision of the measurements in this work are the availability of reliable astrometric calibration sources and their photon statistics. For instance, the uncertainties on CCO proper motion in G15.9+0.2 and Kes 79 are primarily due to the (comparatively) short exposures at the former and latter of the two epochs, respectively, which lead to a small number of photons available for precise astrometric localization. In contrast, for RX J1713.7−3946, we were hindered by the lack of observed X-ray sources with reliable astrometric counterparts in the HRC observation. We argue that the main reason for this is the smaller effective area and higher intrinsic background rate of the Chandra HRC when compared to the ACIS instrument. In principle, the design of the HRC as a microchannel plate instrument (Murray et al. 2000) is ideal for astrometric measurements. However, such measurements can only be performed in an absolute manner in the presence of frame calibration sources, for which the ACIS usually possesses better detection prospects, unless their X-ray emission is very soft.

By combining the measured total proper motion μ_tot with an estimate for the distance d to the SNR, one can constrain the projected velocity of the CCO in the plane of the sky v_proj = μ_totd. This can be viewed as a lower limit to the physical velocity of the NS and thus serves as a proxy to the violent kick experienced by the proto-NS at its birth. Therefore, by comparing the fastest measured NS velocities to theoretical considerations and numerical simulations of core-collapse supernovae, one can attempt to constrain the kick mechanism and the degree of explosion asymmetry.

Currently, the CCO with the largest securely measured velocity is RX J0822−4300 in Puppis A, with a value of v_proj = (763 ± 73) km s⁻¹, scaled to a distance of 2 kpc. As is shown in Table 2, none of the other known CCOs (except possibly those in G15.9+0.2 and Vela Jr. for which sufficient precision has not yet been achieved in a direct measurement) are likely to show a projected velocity in excess of ∼ 800 km s⁻¹. A likely mechanism by which large NS kick velocities are achieved is the “gravitational tug-boat”, in which an asymmetric ejecta distribution accelerates the proto-NS by exerting a gravitational pull in the direction of slow and massive ejecta clumps (Wongwathanarat et al. 2013). As it has been shown that kick velocities above 1000 km s⁻¹ can be achieved in realistic explosion scenarios (Janka 2017), none of the measured CCO velocities are in serious conflict with theoretical expectations.

At the present time, the observed distribution of CCO velocities shows no obvious departure from that for radio pulsars. For instance, assuming a Maxwellian velocity distribution¹⁹ with one-dimensional σ = 265 km s⁻¹ as given by Hobbs et al. (2005), the probability of measuring v_proj ≤ 250 km s⁻¹ is around 35%. Therefore, observing three such cases in our sample is unsurprising. However, the current sample of measured CCO velocities is neither large nor precise enough for definite conclusions on the velocity distribution of CCOs to be drawn.

Most likely, the task of performing a statistically meaningful comparison between kick velocities of radio pulsars and CCOs will require a lot of time and observational effort to complete (including the detection of new CCOs). However, it may ultimately help in shedding light on the question of what is the fundamental difference driving the phenomenological diversity of young NSs. At the present time, it is unclear how the difference in magnetic field strengths between CCOs, rotation-powered pulsars, and magnetars is related to the conditions before, during, or after the supernova explosion.

5.2. SNR expansion and explosion sites

Our proper motion measurements provide a constraint on the CCO’s position over time, back to the explosion date of the supernova. This provides additional input to the physical interpretation of the present-day expansion and morphology of the SNR as the location of its true center can be estimated. With this in mind, we inferred explosion sites of our six SNRs by extrapolating our posterior distribution of the CCO’s proper motion backward from its present-day location, to the assumed supernova explosion time. We then obtained the smallest two-dimensional regions containing 39.3% and 86.5% of the total probability mass, corresponding to 1σ and 2σ constraints on the explosion location, which can be regarded as the true center of the SNR. For this purpose, we assumed the following input ages, along the lines of available constraints (Table 1): 4000 yr (G15.9+0.2), 5000 yr (Kes 79), 340 yr (Cas A), 1000 yr (G330.2+1.0), 2000 yr (RX J1713.7−3946), and 700 yr (G350.1−0.3).

We compare the explosion sites to the present-day morphology of our six SNRs in Fig. 16. The extent of RX J1713.7−3946 is much larger than the Chandra field of view, which is why we used an archival ROSAT PSPC observation (Pfeffermann et al. 1996) to illustrate the characteristic structure of this SNR.

Fig. 16. Constraints on the explosion sites of the six SNRs targeted in this paper. We show exposure-corrected images of the SNRs (with logarithmic intensity scaling), overlaid with 1σ and 2σ contours for the most likely origin of the CCO. The apparent geometrical centers of G15.9+0.2 (estimated in this work), RX J1713.7−3946 (Pfeffermann et al. 1996) and G350.1−0.3 (Lovchinsky et al. 2011) are indicated with white crosses. For Cas A, we plot the precisely constrained expansion center from Thorstensen et al. (2001) in an inset to make it more easily distinguishable. All images were created from the Chandra data used in this work, except for RX J1713.7−3946, where we used an archival ROSAT observation.

5.2.1. G15.9+0.2 and Kes 79: Possible direct evidence for expansion

It is immediately obvious that, for some SNRs, our proper motion measurements are physically much less constraining than for others. For instance, the extent of the uncertainty contours of the explosion sites of G15.9+0.2 and Kes 79 is almost comparable to the extent of the two SNRs themselves. The reason for this is the comparatively shallow exposure of part of the used data sets, in combination with the small intrinsic angular size and large age of the SNRs.

Nevertheless, we measure possible direct signatures of expansion for both remnants, which at this point appear slightly more significant for G15.9+0.2. If the true age of G15.9+0.2 were indeed lower than the 4000 yr that we assumed for estimating the explosion site, its uncertainty contours would “shrink” closer toward the CCO. The resulting constraints would likely still be consistent with what we estimate as the geometric center of the SNR.

For both objects, we demonstrate the feasibility of an exploratory proper motion and expansion study, whose constraints could be significantly improved with future Chandra follow-up observations. For instance, a dedicated observation of Kes 79 around the year 2025 with a depth of ∼30 ks would extend the measurement baseline to ∼25 yr in combination with the 2001 observation. This would be expected to reduce the errors on proper motion and expansion by a factor of ∼2. For G15.9+0.2, given that the very deep late-time observation was carried out in 2015, a followup observation appears to be feasible only from around 2025. For example, a ∼60 ks observation would effect an error reduction by a factor of ∼1.9. In combination, reduced error bars on the motion of both the CCO and the SNR shell would provide insights into the kinematics and temporal evolution of filaments at the shock front as well as the SNR age. Moreover, a more accurate measurement of the proper motion of the CCO in G15.9+0.2 would further constrain its kick velocity. This would allow us to verify if its offset location within the SNR is indeed caused by a rather violent kick or is simply due to asymmetric expansion of the SNR.

5.2.2. G330.2+1.0: Rapid shockwave expansion away from a stationary CCO

In contrast to G15.9+0.2 and Kes 79, we provide stringent upper limits on the proper motion of the CCO of G330.2+1.0. Thus, its present-day location can be regarded as an almost perfect approximation to the SNR expansion center. The most obvious application of this is a comparison with the present-day expansion of the SNR. Borkowski et al. (2018) measure rapid motion away from the CCO, at up to ∼ 9000 km s⁻¹ (for a distance of 5 kpc) for most parts of the SNR shell. In contrast, they detect a far smaller rate of expansion for parts of the southwestern rim, as well as significant temporal decline in its brightness. This, together with the relatively symmetric present-day morphology of the SNR, points toward a quite recent collision of the blast wave with a dense cloud in the southwestern region, causing significant deceleration after almost free expansion for the majority of the SNR’s lifetime (Borkowski et al. 2018).

5.2.3. G350.1−0.3: A very young and highly peculiar SNR

Contrary to all other SNRs that host a CCO, the explosion site of G350.1−0.3 deduced in this work is found to be strongly offset from the remnant’s emission center, but matches the position of the CCO closely (see Fig. 16). Consequently, the striking asymmetry of the SNR is further enhanced by our inferred explosion site: only the rapidly propagating eastern “half” of the remnant appears to be visible at all, while the supposed half west of its true center, if present, remains almost entirely undetected. In any case, conservation of momentum requires the presence of a large amount of mass (either in the form of a compact object or ejecta) moving toward the west.

It is interesting to observe that, when ignoring the presence of the bright eastern clump, the remaining part of the SNR appears to form an almost circular incomplete shell, centered approximately on our inferred explosion site (indicated in Fig. 17). In this context, the bright eastern clump is required to be connected to fast ejecta as it is located far outside this potential shell, and shows the highest physical expansion velocities. As significant expansion is detected in all investigated regions, including the more diffuse northern and southern parts of the shell, it seems conceivable that the initial explosion imprinted a quadrupolar asymmetry on the ejecta distribution.

Fig. 17. Potential diffuse shell of G350.1−0.3. We show the exposure-corrected image of G350.1−0.3 and its explosion site as in Fig. 16, overlaid with the incomplete circular shell tracing the SNR’s diffuse emission (dashed line).

We consider a possible explanation for the missing western SNR shell to be a large density gradient in the surrounding interstellar medium (Gaensler et al. 2008). The alternative, a dipolar asymmetry in the amount of ejected material, would necessarily have imprinted a very strong westward kick on the CCO. Given our measurement of its modest northward proper motion, this seems rather unlikely. The detection of rapid eastward expansion of the shockwave at up to 5000 − 6000 km s⁻¹ here and in Borkowski et al. (2020) only strengthens this conclusion, since it increases the momentum attributed to the bright eastern clump and thus the implied recoil. Furthermore, the expansion measurement leads to a precisely constrained explosion site, due to the now certain very low age of the SNR. Therefore, it can be stated with high confidence that the circumstellar material around G350.1−0.3 is highly inhomogeneous: It likely exhibits a density gradient from east to west and a significant localized density enhancement in the east. This eastern clump causes, as we have shown, a local velocity gradient in the expanding ejecta, which manifests itself in the fascinating outward-bent shape of the interacting shock wave.

Given our observations, we can infer that the true extent of the remnant is necessarily greater than apparent in X-rays, with a radius of at least 3.5′ (the approximate distance between the bright eastern ejecta clumps and the current location of the CCO), larger than the 2.5′ estimated in Lovchinsky et al. (2011). Since the remnant is larger, but also expanding much faster, than predicted, their age estimate of 600 − 1200 years is still consistent with the 600 − 700 years found by Borkowski et al. (2020) and deduced in our analysis.

The overall multiwavelength morphology of G350.1−0.3 is quite rich, and can aid in understanding its peculiar appearance in X-rays. As shown by Lovchinsky et al. (2011), the bright eastern clump has an obvious counterpart in the mid-infrared (24 μm), likely due to dust shocked by the interaction of the supernova ejecta with a molecular cloud. In its fainter regions, the 24 μm image of the SNR shows features corresponding to the southern part of the diffuse X-ray shell, but little apparent emission in the north. In contrast, as displayed by Gaensler et al. (2008), the SNR as viewed in the radio domain (at 4.8 GHz) exhibits a horizontal “bar,” coincident with the northern part of the shell seen in X-rays. In turn, there appears to be no radio emission originating from the southern region. Keeping these points in mind, it is interesting to observe that the X-ray emission in these two regions also shows spectral differences, with the north emitting harder and the south softer radiation (see Fig. 1 in Borkowski et al. 2020). We interpret these findings as a further signature of inhomogeneous circumstellar medium, where denser circumstellar material in the south causes stronger interaction with the supernova blast wave, and thus leads to infrared emission from shocked or heated dust. Less dense material in the north interacts weakly with the shock wave, allowing it to propagate more freely, consistent with the observed radio and harder X-ray emission. Borkowski et al. (2020) come to a related conclusion, demonstrating that the northern part of the SNR shows little evidence for the presence of ejecta, and its emission is thus likely attributable to the propagation of the supernova blast wave through a comparatively thin medium.

Given the youth of G350.1−0.3, we checked the catalog of “guest stars” by Stephenson & Green (2009) for any potential historical counterpart. There is no obvious match to our SNR, except possibly one reported guest star in AD 1437, which would imply an age ∼580 yr. The specified position of this historic event is (α, δ) = (16^h55^m, −38° ), which places it about five degrees from the position of G350.1−0.3 (Stephenson & Green 2009). This spatial discrepancy and the possible nova explanation preferred by the authors argue against an association of this guest star with G350.1−0.3. In fact, the observed hydrogen column density N_H ∼ 4 × 10²² cm⁻² (Lovchinsky et al. 2011) implies a visual extinction by around A_V ∼ 22 mag, according to the relation of Predehl & Schmitt (1995). Combining this with an assumed peak absolute magnitude M_V ∼ −18 for a core-collapse supernova at a distance of 4.5 kpc results in an observed apparent magnitude m_V ∼ 18. It therefore seems very unlikely that the supernova associated with G350.1−0.3 would have been observable with the naked eye at all.

5.2.4. RX J1713.7−3946: Evidence for nonuniform expansion history

For RX J1713.7−3946, we have determined that the supernova explosion site is likely located within ∼2′ of the present-day location of the CCO, which is somewhat at odds with the apparent geometrical center (Pfeffermann et al. 1996). The expansion of the SNR has been independently measured along three sectors of its shell, in the southeast, southwest, and northwest in three individual studies (Acero et al. 2017; Tsuji & Uchiyama 2016; Tanaka et al. 2020). Interestingly, these works all show that the maximal expansion speed of the shockwave is remarkably uniform, on the order of 4000 km s⁻¹ (or 800 mas yr⁻¹), in all directions²⁰.

With our determination of the expansion center of the SNR, we can provide updated constraints on the age of RX J1713.7−3946, by combining the location of its true center with the proper motion of the fastest “blobs” in the three regions, in analogy to Sect. 4.6: Comparing our inferred explosion site with the estimated position of “blob A” in the southwest (Tanaka et al. 2020), we find the angular distance covered by the shock wave up to today to be around ϑ = 1290″ ± 50″. Together with its measured proper motion of μ_exp = (810 ± 30) mas yr⁻¹, this yields an undecelerated expansion age of τ = (1590 ± 80) yr. Since it is likely that the shock wave has been somewhat decelerated during its expansion history, we quote only the corresponding 90% upper limit on the expansion age of the SNR at τ < 1700 yr. An analogous investigation using the reported values and errors for the northwest (Tsuji & Uchiyama 2016) and the southeast (Acero et al. 2017) rims yields higher formal limits of τ < 2520 yr and τ < 3240 yr, respectively, as illustrated in Fig. 18. We should note here that assuming a different age as input to the computation of the SNR’s explosion site would slightly alter the size of the uncertainties on its expansion age. However, it would have little effect on the overall result since only significant nonzero motion of the CCO would bias our computation in any direction.

Fig. 18. Expansion and age estimate for RX J1713.7−3946. Top panel: measured maximal proper motion μexp of filaments in the expanding shell versus their angular distance ϑ from the explosion site. We indicate values from the literature for the southeast (Acero et al. 2017), northwest (Tsuji & Uchiyama 2016), and southwest (Tanaka et al. 2020) regions, respectively. For the proper motion, we display the error bars as given in the respective paper, scaled to 68% confidence and including systematic errors for Acero et al. (2017). Bottom panel: derived likelihood and corresponding 90% confidence upper limits for the undecelerated expansion age τ, inferred for the respective filaments.

Therefore, if we assume all input expansion velocities and errors to be reliable, we can infer that the shock wave must have been decelerated nonuniformly, leading to the by far largest present-day expansion rate inferred for the southwest filament. Generally, it can be considered unlikely that any physical mechanism could have substantially accelerated the blast wave in the past, while it is very probable for the shock wave to have experienced some degree of deceleration due to past interaction with the surrounding medium. Therefore, we choose to quote the smallest out of the three measurements of τ, corresponding to a true age of < 1700 yr. This is a quite constraining result, as other recent age estimates are in the range 1500 − 2300 yr (Acero et al. 2017; Tsuji & Uchiyama 2016), depending on the exact model assumed for the density profiles of ejecta and the interstellar medium. In particular, our value is still consistent with the tentative association of RX J1713.7−3946 with the historical guest star, a possible supernova, observed by Chinese astronomers in the year 393 (Wang et al. 1997), if we assume parts of the southwest shell to be expanding almost freely.

However, it is somewhat puzzling that the shock wave seems to show the least decelerated expansion along the southwestern direction, in particular for blob A. The observed X-ray emission from blob A is much softer than its surroundings with a power law spectral index Γ = 2.74 ± 0.07 (Tanaka et al. 2020). Using the simplest assumptions, this implies a shock velocity of only ∼ 2800 km s⁻¹ (Okuno et al. 2018), around 1000 km s⁻¹ smaller than the measured value. Additionally, X-ray and molecular gas observations imply that the shock wave is expanding in a cavity blown by the progenitor in the southeast, whereas it has likely relatively recently struck the dense cavity wall in the west (Cassam-Chenaï et al. 2004; Fukui et al. 2012). Thus, one would naively expect to measure the least decelerated expansion of the shock wave in the southeast of RX J1713.7−3946. We therefore emphasize that our age estimate is largely based on the published proper motion of only a single blob in the southwest of the SNR by Tanaka et al. (2020), and is therefore to be taken with caution. In order to confirm these constraints, future followup observations and/or independent reanalyses are needed to verify both the proper motion of the CCO and that of blob A.

In this context, we consider again the double ring-like morphology of RX J1713.7−3946 in X-rays, which is most striking for the western shell (see Fig. 16). As pointed out by Cassam-Chenaï et al. (2004), both the inner and outer apparent ring exhibit a deformed elliptical structure. This structure appears nearly continuous around the entire SNR, and for the inner part, is matched well by the SNR’s radio morphology. It is interesting to observe that the southern and eastern portions of this inner ring appear to be located quite close to our derived explosion site, as could be expected from the reverse shock. However, due to the overall weakness of thermal X-ray emission (Katsuda et al. 2015) and the radio detection of the inner ring, such an interpretation seems quite unlikely for this feature. Alternatively, one could also imagine that the complicated morphology stems from the superposition of two unrelated SNRs. This however can be considered unlikely due to the observed dominance of nonthermal X-ray emission across the entire remnant (e.g., Okuno et al. 2018) and the observed expansion in multiple regions. The most likely scenario is that the SNR’s peculiar morphology is caused by the projection of different parts of the shock wave, expanding into an anisotropic wind-blown cavity from the progenitor star. Locally varying collision times with the cavity wall are also a reasonable explanation for the observed differences in the relative expansion rate of the SNR.

5.2.5. Cas A: Adding independent constraints to optical expansion measurements

Finally, for Cas A, the situation is quite unique: While X-ray measurements of the propagation of the supernova shock wave do exist (e.g., Vink et al. 1998; Patnaude & Fesen 2009), their general target was to trace the exact location and propagation of the forward shock along the SNR shell. These studies were able to demonstrate that the forward shock of Cas A travels at a speed around (4200 − 5200) km s⁻¹, implying significant deceleration of the shockwave compared to free expansion. Direct expansion and age measurements for this infant SNR can be performed with much higher precision by tracing the motion of clumpy filaments in optical line emission. Prominent examples are the works by Fesen et al. (2006) and Thorstensen et al. (2001), who determined the explosion date of Cas A to around the year 1670.

Our constraint on the explosion site via NS proper motion is presently far less precise than the constraints from optical filament expansion. However, the accuracy of the latter result depends on the presence of systematic effects, for example from nonsymmetric deceleration of ejecta knots. If we imagine being able to use a potential future observation to extend our temporal baseline from around 10 to more than 20 years, we can envisage greatly reduced error bars on our X-ray measurement of the CCO’s proper motion. For instance, we estimate that a further HRC observation of similar depth as the 1999 data set (∼50 ks) would reduce the uncertainty on the proper motion by a factor of ∼2.5. The trajectory of the NS, which is almost perpendicular to the motion of the fastest optical knots (Fesen et al. 2006), could then serve as a complementary and completely independent ingredient to measuring the age and location of the expansion center of Cas A. The main advantage is that, due to its enormous density, the NS can be regarded as a “bullet,” moving through its surroundings at effectively constant velocity. This makes its proper motion a useful tool to search for systematics in the motion of the fastest optical ejecta knots, possibly improving the accuracy of the age estimate for Cas A.

As pointed out by Fesen et al. (2006), it is interesting to observe that the NS in Cas A appears to be moving in a direction unrelated to that of the prominent “jets” of high-velocity ejecta. Instead, the CCO moves into the direction roughly opposite the brightest ejecta emission (Holland-Ashford et al. 2017). This is likely a manifestation of the NS kick mechanism, in which the proto-NS experiences the “gravitational tug” of slow and massive ejecta during the first seconds of the explosion (Wongwathanarat et al. 2013).

6. Summary

In this work, we analyzed several sets of archival Chandra observations with the aim of deducing proper motion measurements consistently for all known CCOs with suitable data. Important elements of our analysis were the systematic search for serendipitous sources for astrometric frame calibration, the fitting of models of the Chandra PSF to the data to obtain unbiased source positions, and the alignment of measured positions with a calibrated reference frame. The final results were obtained by a simultaneous fit to the data from all observational epochs, to provide accurate estimates of the proper motion of the target. In this way, we attempted to exclude biases by PSF morphology, aspect reconstruction uncertainties, or the effects of Galactic rotation.

In total, we present six new measurements of NS proper motion, four of which had not been directly targeted prior to this work. Thus, we approximately double the existing sample of CCOs with quantitatively constrained kinematics. We find moderately significant nonzero proper motion only for two objects in our sample, the NSs in Cas A and G350.1−0.3, with no evidence for any hypervelocity CCOs. In contrast, we set comparatively tight upper limits on the motion of several objects in our sample. For instance, the transverse physical velocities of the CCOs in G330.2+1.0 and RX J1713.7−3946 are constrained at < 230 km s⁻¹ at distances of 5 kpc and 1 kpc, respectively.

Complementarily, we used our astrometrically coaligned data sets to directly measure the expansion of three SNRs, by comparing the one-dimensional emission profiles of characteristic features at different epochs. For the SNRs G15.9+0.2 and Kes 79, we applied our methods to comparatively shallow data sets to systematically search for signatures of expansion. For both objects, we find tentative evidence for the expansion of fragments of the SNR shell away from its center, at 2 − 2.5σ formal statistical significance. However, at the present time, the suboptimal exposure times and pointing configurations of the archival data sets prohibit a more detailed analysis of internal kinematics or an exact age estimate.

For the much younger SNR G350.1−0.3, we qualitatively confirm, using our independent methods, recently published results by Borkowski et al. (2020), most importantly the SNR’s rapid expansion at close to 6000 km s⁻¹. Conservative interpretation of our measurement implies an upper limit on the SNR’s age of 700 years, which places G350.1−0.3 among the three youngest Galactic core-collapse SNRs. Additionally, we directly confirm that kinematics within the SNR vary on small scales within its bright eastern clump, likely due to the strong interaction of ejecta with a molecular cloud. Furthermore, expansion is visible also for much fainter emission features, proving that the peculiar observed morphology of G350.1−0.3 is in fact not due to chance superposition of two SNRs.

We used our CCO proper motion measurements to constrain the exact explosion locations of the six SNRs in our sample. For instance, for the young SNR G330.2+1−0, we show that its supernova occurred within only 10″ of the present-day location of the CCO. A very interesting finding is the confirmation of a striking asymmetry in the SNR G350.1−0.3, whose explosion site is located far west of its apparent geometric center. This finding has severe implications on the properties of the SNR, requiring an ultra-asymmetric explosion or strong density inhomogeneities in the circumstellar environment of one of the youngest and most bizarre Galactic SNRs.

Finally, we combined several published expansion measurements of RX J1713.7−3946 with our estimate of the supernova explosion site. We find that the current expansion rate significantly differs between individual regions. In particular, we show that part of the southwestern shock wave has likely experienced much weaker past deceleration than other regions, since otherwise its present-day velocity could not be reconciled with its origin. Using the most stringent resulting constraint, we determine an upper limit of 1700 years on the age of RX J1713.7−3946.

Some of our findings, as any multi-epoch astrometric analysis, would certainly profit from future Chandra follow-up observations, ideally similar to the 19-year campaign to measure the proper motion of the CCO in Puppis A (see Mayer et al. 2020). For instance, a future Chandra observation to be taken a few years from today would allow us to constrain the age and NS kick velocity of G15.9+0.2 and Kes 79, via an analysis of the kinematics of the CCO and the SNR shell. Additionally, a more precise trajectory of the NS could be used to indirectly constrain deceleration of optical ejecta filaments in Cas A, identifying potential systematic errors on its age estimate. Similarly, combining X-ray measurements of expansion for RX J1713.7−3946 with a more precise measurement of the proper motion of its CCO may help in supporting or rejecting its association with the historical SN 393. We hope that some of these exciting questions can be addressed during the lifetime of Chandra, as it is the only instrument allowing for the required high precision in X-ray astrometry – now, as in the foreseeable future.

Acknowledgments

We would like to thank Johannes Buchner for helpful discussions on statistical problems. We are grateful to the anonymous referee for the provided constructive criticism and useful suggestions. We appreciate the reliable support by the Chandra helpdesk at several stages of our analysis. MGFM acknowledges support by the International Max-Planck Research School on Astrophysics at the Ludwig-Maximilians University (IMPRS). We acknowledge the use of the Chandra data archive and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO and Sherpa. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). Further, we acknowledge the use of the Python packages Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), and NumPy (Harris et al. 2020).

References

Appendix A: Observation log and astrometric calibrators

Appendix B: Full results of proper motion and expansion measurements

Appendix C: Illustration of SNR expansion measurements in individual regions

Fig. C.1. Direct comparison of flux profiles in different regions of G15.9+0.2. Each panel is labeled with the letter corresponding to the respective region. Top and middle sub-panels: exposure corrected images at the early and late epoch, and the white rectangle outlines the measurement region. Bottom sub-panel: one-dimensional flux profiles Fi(x) at the early (red) and late (blue) epoch. The x coordinate increases outward along the radial direction of the SNR. The markers display the binned distribution, and the solid line reflects the flux profile smoothed with a Gaussian kernel. The dashed red line corresponds to the profile of the early epoch, shifted and rescaled by the parameters indicated in the upper right corner, with the shaded area displaying the associated uncertainty.

Fig. C.2. Direct comparison of flux profiles in different regions of Kes 79 (as in Fig. C.1).

Fig. C.3. Direct comparison of flux profiles of G350.1−0.3 (as in Fig. C.1).

Fig. C.3. continued.

All Tables

All Figures

	Fig. 1. Overview over all target SNRs, with positions of the NS and the calibrator sources (labeled as in Table A.2) indicated on exposure-corrected images of the respective SNR. We show the deepest available observation for all targets except for RX J1713.7−3946, where we used the ACIS observation (ID 5559) instead of the HRC observation, due to its lower intrinsic noise level.
In the text

Fig. 2. Illustration of PSF-fit likelihood for the source position and the prediction by our simultaneous proper motion fit for the CCO of G15.9+0.2. Each column shows data from the individual epochs (with the observation ID displayed in the title), while each row represents one particular object. The one to five sigma contours of the PSF-fit likelihood ℒij(x, y) are shown in green, and the predicted positions obtained from our simultaneous fit (see posterior distribution in Fig. 4) are indicated using a grayscale two-dimensional histogram. In the lower left corner of each panel, we provide a reference scale of length 0.3″. We display the true source position for our calibrator stars (known from Gaia) with a blue star. The uncertainty contours for observation 16766 are by far the smallest since the corresponding exposure was much deeper than at the other epochs (at ∼90 ks).

In the text

Fig. 3. Example expansion measurement in region D of G350.1−0.3 (see Sect. 4.6). Top and middle panels: exposure corrected images at the early and late epoch, and the white rectangle marks the measurement region. Bottom panel: one-dimensional flux profiles Fi(x) at the early (red) and late (blue) epoch. The markers display the binned distributions with errors, and the solid lines reflect the flux profiles smoothed by a Gaussian kernel. The dashed red line corresponds to the profile of the early epoch, shifted and rescaled by the indicated parameters (see upper right corner), with the shaded area displaying the associated uncertainty from bootstrapping. We define the x coordinate such that it increases in radial direction, meaning from west to east in this case.

In the text

Fig. 4. Corner plot displaying the joint and marginal posterior distributions for the parameters of the astrometric solution for the CCO in G15.9+0.2. The astrometric calibration parameters {Δxi, Δyi, ri, θi} have been marginalized over, and are not shown here in order to keep the figure uncluttered (see Table B.1 for constraints on all parameters). In the two-dimensional correlation plots, we indicate the 1σ, 2σ, and 3σ contours, meaning the smallest regions containing 39.3%, 86.5%, and 98.9% of the total probability mass, with solid lines. In the marginalized one-dimensional histograms, we indicate the median and 68% central interval of the posterior probability distribution, with dashed lines. The units of proper motion and the astrometric zero-point are mas yr−1 and arcsec, respectively. This plot and analogous figures have been created using the corner.py package (Foreman-Mackey 2016).

In the text

Fig. 5. SNR G15.9+0.2 with regions used for the expansion measurement indicated. We show a smoothed, exposure-corrected image of the SNR with logarithmic scaling of the color map. We overlay the boxes from which we extracted the one-dimensional emission profiles for our expansion measurement, and label them alphabetically for easier identification and to avoid confusion with numeric labels for the astrometric calibrator stars.

In the text

Fig. 6. Constraints on the expansion of G15.9+0.2. Top: for each feature (labels as in Fig. 5), we plot the measured angular expansion speed μexp, and the corresponding projected velocity vexp for a distance of 10 kpc, against its angular distance to the SNR center ϑ. To assess their respective significance, we have used black error bars for those measurements that enclose > 95% of the total probability mass on either side of 0, and gray for the others. In addition, we indicate the expected free expansion speed for ages of 1000, 2000, and 3000 years with dashed lines. Bottom: expansion rate τ−1 as measured in each individual region. The thick black line with the gray-shaded area represents the combination of the individual measurements (or weighted average), assuming uniform expansion. All error bars are at 68% confidence. The data underlying this figure are given in Table B.2.

In the text

	Fig. 7. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in Kes 79. See Table B.1 for constraints on all parameters.
In the text

	Fig. 8. SNR Kes 79 with regions used for the expansion measurement indicated (as in Fig. 5).
In the text

	Fig. 9. Constraints on the expansion of Kes 79 (as in Fig. 6). The assumed distance to Kes 79 for the conversion to vexp is 5.0 kpc. The data underlying this figure are given in Table B.2.
In the text

	Fig. 10. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in Cas A. See Table B.1 for constraints on all parameters.
In the text

	Fig. 11. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in G330.2+1.0. See Table B.1 for constraints on all parameters.
In the text

Fig. 12. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution for the CCO in G347.3−0.5. For this target, a relative astrometric frame registration method was applied as described in the text. We therefore measured no absolute positions, but only relative source displacements, which correspond to the components of proper motion. See Table B.1 for constraints on all parameters.

In the text

	Fig. 13. Corner plot (as in Fig. 4) showing the posterior distribution of the astrometric solution of the CCO in G350.1-0.3. See Table B.1 for constraints on all parameters.
In the text

Fig. 14. SNR G350.1−0.3 with regions used for the expansion measurement indicated. Top panel: as in Fig. 5. Bottom panel: difference image between the early and late epochs. We use a symmetric logarithmic color scale to highlight features on small scales, with blue corresponding to the late-epoch emission being brighter, and red being the opposite. For the main image, we smoothed the data with a Gaussian kernel of 1″, while in the inset we show the region of the bright clump smoothed with only 0.5″.

In the text

	Fig. 15. Expansion and age estimate for G350.1−0.3. Top panel: as in Fig. 6, with an assumed distance of 4.5 kpc. Bottom: constraints on the expansion age τ from each individual region, with error bars at 68% confidence. The data underlying this figure are given in Table B.2.
In the text

Fig. 16. Constraints on the explosion sites of the six SNRs targeted in this paper. We show exposure-corrected images of the SNRs (with logarithmic intensity scaling), overlaid with 1σ and 2σ contours for the most likely origin of the CCO. The apparent geometrical centers of G15.9+0.2 (estimated in this work), RX J1713.7−3946 (Pfeffermann et al. 1996) and G350.1−0.3 (Lovchinsky et al. 2011) are indicated with white crosses. For Cas A, we plot the precisely constrained expansion center from Thorstensen et al. (2001) in an inset to make it more easily distinguishable. All images were created from the Chandra data used in this work, except for RX J1713.7−3946, where we used an archival ROSAT observation.

In the text

	Fig. 17. Potential diffuse shell of G350.1−0.3. We show the exposure-corrected image of G350.1−0.3 and its explosion site as in Fig. 16, overlaid with the incomplete circular shell tracing the SNR’s diffuse emission (dashed line).
In the text

Fig. 18. Expansion and age estimate for RX J1713.7−3946. Top panel: measured maximal proper motion μexp of filaments in the expanding shell versus their angular distance ϑ from the explosion site. We indicate values from the literature for the southeast (Acero et al. 2017), northwest (Tsuji & Uchiyama 2016), and southwest (Tanaka et al. 2020) regions, respectively. For the proper motion, we display the error bars as given in the respective paper, scaled to 68% confidence and including systematic errors for Acero et al. (2017). Bottom panel: derived likelihood and corresponding 90% confidence upper limits for the undecelerated expansion age τ, inferred for the respective filaments.

In the text

Fig. C.1. Direct comparison of flux profiles in different regions of G15.9+0.2. Each panel is labeled with the letter corresponding to the respective region. Top and middle sub-panels: exposure corrected images at the early and late epoch, and the white rectangle outlines the measurement region. Bottom sub-panel: one-dimensional flux profiles Fi(x) at the early (red) and late (blue) epoch. The x coordinate increases outward along the radial direction of the SNR. The markers display the binned distribution, and the solid line reflects the flux profile smoothed with a Gaussian kernel. The dashed red line corresponds to the profile of the early epoch, shifted and rescaled by the parameters indicated in the upper right corner, with the shaded area displaying the associated uncertainty.

In the text

	Fig. C.2. Direct comparison of flux profiles in different regions of Kes 79 (as in Fig. C.1).
In the text

	Fig. C.3. Direct comparison of flux profiles of G350.1−0.3 (as in Fig. C.1).
In the text

	Fig. C.3. continued.
In the text