A kinematic study of central compact objects and their host supernova remnants

Context. Central compact objects (CCOs) are a peculiar class of neutron stars, primarily encountered close to the center of young supernova remnants (SNRs) and characterized by thermal X-ray emission. Aims. Our goal is to perform a systematic study of the proper motion of all known CCOs with appropriate data available. In addition, we aim to measure the expansion of three SNRs within our sample to obtain a direct handle on their kinematics and age. Methods. We analyze multiple archival Chandra data sets, consisting of HRC and ACIS observations separated by temporal baselines between 8 and 15 years. In order to correct for systematic astrometric uncertainties, we establish a reference frame using X-ray detected sources in Gaia DR2, to provide accurate proper motion estimates for our target CCOs. Complementarily, we use our coaligned data sets to trace the expansion of three SNRs by directly measuring the spatial offset of various filaments and ejecta clumps between different epochs. Results. In total, we present new proper motion measurements for six CCOs, among which we do not find any indication of a hypervelocity object. We tentatively identify direct signatures of expansion for the SNRs G15.9+0.2 and Kes 79, at estimated significance of $2.5\sigma$ and $2\sigma$, respectively. Moreover, we confirm recent results by Borkowski et al., measuring the rapid expansion of G350.1$-$0.3 at almost $6000\,{\rm km\,s^{-1}}$, which places its maximal age at $600-700$ years. The observed expansion, combined with the rather small proper motion of its CCO, implies the need for a very inhomogeneous circumstellar medium to explain the highly asymmetric appearance of the SNR. Finally, for the SNR RX J1713.7$-$3946, we combine previously published expansion measurements with our measurement of the CCO's proper motion to obtain a constraining upper limit of $1700$ years on the system's age.


Introduction
It is generally accepted that collapsing massive stars, i.e. corecollapse 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 which allows to set a lower limit 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 a 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 e-mail: mmayer@mpe.mpg.de et al. (2005) and Verbunt et al. (2017). Such works establish the picture of NSs being very high-speed objects, with a mean threedimensional velocity around 400 km s −1 , and reliably measured projected velocities at least as high as ∼ 800 km s −1 .
In contrast to radio pulsars, for NSs exhibiting at best weak radio and optical emission, such as magnetars, X-ray-dim isolated neutron stars, or central compact objects in supernova remnants (CCOs), such measurements are much more challenging. While all three of these neutron star classes are characterized by primarily thermal X-ray emission, CCOs are particularly noteworthy for being exclusively young ( 10 4 yr) and nonvariable X-ray sources, located close to the center of supernova remnants. 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 neutron star atmosphere models with luminosities 10 33 erg s −1 , 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. 2008Mignani et al. , 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 1 of known CCOs in Table 1. Out of ten known CCOs, only for those in Kes 79, Puppis A, and PKS 1209−51/52, X-ray pulsations have been detected. These pulsations at 100 − 500 ms are likely attributable to the rotational modulation of the emission from hotspots on the neutron star surface (Gotthelf et al. 2010(Gotthelf et al. , 2005Zavlin 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, lead 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 apparent nearby SNR but in a similar location of (P,Ṗ)-space, are found (Luo et al. 2015). It has therefore been suggested that the phenomenology of CCOs is caused by magnetic field burial due to 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 ∼ 10 4 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 −1 , at 68% confidence. This corresponds to a projected velocity of (763 ± 73) km s −1 for a distance of 2 kpc, which potentially places this CCO at the fast end of the NS velocity distribution. 2 In contrast, for the CCO in PKS 1209−51/52, Halpern & Gotthelf (2015) measured a much smaller proper motion of (15 ± 7) mas yr −1 , which corresponds to a projected velocity below 180 km s −1 for a distance of 2 kpc. Only weak constraints could 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) showed 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 a 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 a 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) explored the relationship between the proper motions of 18 young neutron stars (four of which are CCOs), and the morphology of their host SNRs. Several of their proper motion measurements were obtained directly from X-ray imaging data, with others (including those for all four CCOs) inferred indirectly. They did not find any correlation between the magnitude of SNR asymmetry and the projected velocity of the neutron stars. However, they found 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) confirmed this observation. In addition, the latter work did 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 a 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. 3 Tracing the expansion of a SNR is in principle possible at all wavelengths which 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 north-western 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 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.
CCOs together with their parent SNRs present an ideal target for combining kinematic information from both the neutron star 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.

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 4 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 which 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). 5 This leaves 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 the appendix 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 which contain the CCO, taken with different detectors. This makes a 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.

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).

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. 6 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 −5 . 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 eliminates 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 7 (Gaia Collaboration et al. 2016, 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 which we deemed unreliable, for example because of large positional errors in X-rays or doubts in their correct optical identification. If 6 https://cxc.harvard.edu/cal/ASPECT/celmon/ 7 https://www.cosmos.esa.int/web/gaia/dr2 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 the appendix in Table A.2, where we also indicate which of the sources we decided to exclude in subsequent steps.

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

PSF modelling and fitting
The first step is 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 8 and the MARX software package 9 (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). 10 For each source in each observation, we obtained a profile of the logarithmic likelihood L i j (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).

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 allows 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 are unknown: 11 8 http://cxc.harvard.edu/chart/index.html 9 http://space.mit.edu/CXC/MARX/ 10 http://cxc.harvard.edu/sherpa/ 11 Of course, other parameters such as source count rates and background levels are also unknown. However, for our purpose, the only globally interesting parameters are related to source positions.   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.
-The astrometric solution of the CCO with the set of free parameters {µ α , µ δ , α 0 , δ 0 }, 12 describing its proper motion and its position at a reference epoch t 0 , in right ascension (α) and declination (δ), respectively. -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 make 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 13 calibrator source positions are effectively known prior to the observation. Therefore, for each of our targets, there is a set Λ = {µ α , µ δ , α 0 , δ 0 , ∆x 1 , ∆y 1 , ..., r N , θ N } of 4 + 4 × N parameters (with N the number of observations included), which are a priori unknown.
In order to constrain their values, we calculate the logarithmic likelihood for a given set of parameters, L tot (Λ), in the following manner: For a given astrometric solution {µ α , µ δ , α 0 , δ 0 }, we compute the true position of the CCO (x i0 , y i0 ), at the epoch 12 We follow the convention that µ α =α cos(δ), so that the proper motion is directly proportional to the object's physical velocity. Positive µ α corresponds to motion from west to east. 13 In the following, we refer to a position calibrated to the ICRS as "true" position, while a "measured" position describes the position determined in the coordinate system of a particular Chandra observation. of observation i according to: where t i is the time of observation i. The fixed parameter t 0 describes the reference epoch, for which we always chose the time of the most recent observation in our respective data set. The true positions (x i j , y i j ) of the calibrator sources are computed analogously, using their known Gaia astrometric solutions. These true positions are then converted to measured positions (x i j , y i j ), to account for a coordinate system shift modelled by the parameters {∆x i , ∆y i , r i , θ i }: The likelihood of a given measured position is evaluated by interpolating the logarithmic likelihood grid L i j (x, y), which we compute 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 Λ can then be straightforwardly calculated: Article number, page 5 of 30  proofs: manuscript no. main 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 14 software (Buchner 2021), which employs the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016(Buchner , 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 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 i j , y i j ). The distribution of these locations was then compared to the measured uncertainty contours, or equivalently the contours of L i j (x, y). This allowed us to check whether the results given by our simultaneous fit are 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.

Correction for effects of Solar motion and Galactic rotation
Ultimately, we are attempting to extract physically meaningful quantities, such as neutron star 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 co-rotation with the disk was expressed as two components of proper motion in equatorial coordinates, (µ LSR α , µ LSR δ ). From its distance d, proper motion, and LSR, it is then possible to determine the peculiar (projected) velocity of an object, v * proj as  Fig. 4) are indicated using a grayscale 2D histogram. The grid in each plot is spaced by 0.3 in right ascension and declination. 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).
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 (µ LSR α , µ LSR δ ) 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 µ * tot is biased toward large values. This is because, in the transformation to µ * tot , one essentially integrates over concentric circles of increasing radius, introducing an effective prior ∝ µ * tot , and always leading to a peak of the distribution at non-zero µ * tot . Therefore, in order to not overestimate the projected velocity of our CCOs, we reweighted the posterior distributions by a factor ∝ 1/µ * tot . The effect of this correction decreases with increasing significance of non-zero proper motion.

Expansion measurement of SNRs
When studying the dynamics of a young core-collapse SNR, a complementary tool to the analysis of the neutron star's proper motion is the measurement of the SNR's expansion. It provides Article number, page 6 of 30 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 has the advantage that, for observations taken close to each other in time, the position of the CCO (which is often the brightest source) is forced to be identical. This practically provides 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 & Mandel 2003) to interactively define rectangular regions containing emission features which we found suitable to trace the expansion of the SNR, for example sharp filaments or individual clumps. Since we are 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. 15 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 χ 2 ) 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, 15 Note that even for potential misalignment with the true direction of motion as large as 20 • , the measured expansion speed would only deviate by 1 − cos(20 • ) = 6% from the true value. within the rectangular region of each feature. For both epochs, we created N = 1000 randomized realizations of the observed data set {x * 1,i , ..., x * N,i }, via sampling with replacement from the list of observed events in the region. Thereby, the statistical uncertainty of the emission profile is approximated by the intrinsic scatter of our resampled data sets.
From each resampled data set x * n,i , we extracted onedimensional 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 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: 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 = µ exp d, where d corresponds to the (assumed) distance to the SNR.
The last step of estimating the global minimum of f is 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 such 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.
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 latetime 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 found 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, χ 2 = K k=1 ξ 2 k /∆ξ 2 k , was found to be χ 2 = 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 obtained χ 2 = 14.3. The expected 68% central interval of a χ 2 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.

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 north-western "blowout" region with only weak apparent emission, likely due to reduced density of the circumstellar medium ). Its distance is only constrained to a relatively uncertain range of values ( Table B.1 for constraints on all parameters). In the 2D 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 1D 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).
and accordingly, the age shows correlated uncertainties of factor 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 west or south-west. 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 −1 for an assumed age of 3500 years.
Unfortunately, the available data at the early epoch are split into three non-consecutive 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 −1 , 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.
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 −1 toward west seems realistic given the SNR geometry. It would imply a rather large projected velocity of the NS around 700 km s −1 when scaling to a distance of 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.
10 kpc. The formal 90% upper limit of 25 mas yr −1 (corresponding to 1200 km s −1 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, which we both avoided with our choice of regions.
In order to define a rough reference point for the expansion of the SNR, we used SAOImage ds9 (Joye & Mandel 2003) to approximate the shell of the SNR with a circle. Its center at (α, δ) = (18 h 18 m 54. s 2, −15 • 01 55 ) 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 ϑ travelled by the shock wave at the location of each region.
The results of our expansion analysis are displayed in 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 −1 . Therefore, the overall effect of this  Fig. 5), we plot the measured angular expansion speed µ exp , and the corresponding projected velocity v exp 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 which 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 is given in Table B.2. 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. We attempt to give a quantitative estimate of significance and rate of expansion now: 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 τ −1 = µ 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 τ −1 afterwards. 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 τ −1 = (4.3 ± 1.9) × 10 −4 yr −1 . This finding is moderately signif-Article number, page 9 of 30 A&A proofs: manuscript no. main µ α = − 2.6 +10.5  Table B.1 for constraints on all parameters.
icant 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 non-zero 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 2D shift of the coordinate system). For all these cases, we found that the finding of 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.

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 midinfrared 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 was inferred from H i absorption (Ranasinghe & Leahy 2018), yielding a distance of 3.5 kpc, whereas CO observations ) suggest a distance of 5.5 kpc. The age of Kes 79 has been estimated via detailed X-ray spectroscopy to lie between 4.4 and 6.7 kyr . Our measurement of the CCO's proper motion (Fig. 7) is perfectly consistent with zero, the 68% central credible intervals being (−3 +11 −10 , −3 +12 −11 ) mas yr −1 . This can be converted to a 90% upper limit on the CCO's peculiar proper motion of around 19 mas yr −1 . Equivalently, its tangential velocity is constrained at < 450 km s −1 , when scaled to a distance of 5 kpc.
While the limit on the CCO's transverse velocity is not very constraining by itself, our measurement has a noteworthy implication on the neutron star'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 −19 s s −1 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. 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 increase 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.
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 0, and several regions appearing almost unconstrained. However, the important exception is filament A as it, on its own, shows evidence for non-zero proper motion at ∼ 4σ statistical signifi-Article number, page 10 of 30  Fig. 6). The assumed distance to Kes 79 for the conversion to v exp is 5.0 kpc. The data underlying this figure is given in Table B.2. cance. 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 −1 . 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. Analogously to Sect. 4.1, we combined the measurements of individual regions to constrain the uniform expansion rate of Kes 79 at τ −1 = 1.7 +0.9 −0.8 × 10 −4 yr −1 . 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 non-zero 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 which 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 non-uniform 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.

Cas A
Cas A (G111.7-2.1) is the youngest known Galactic corecollapse 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 −1 (De-Laney 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 .
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 yielded a final measurement of the transverse velocity of (390 ± 400) km s −1 . This is equivalent to a total proper motion of (24 ± 25) mas yr −1 toward south-west, for their assumed distance of 3.4 kpc. Nevertheless, since they used different statistical methods than we do, and did not employ any PSF modelling, 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 (labelled 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 which are 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 detected 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 Figure 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 to 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 did not observe this, we concluded that all three astrometric calibrators are probably correctly identified.
The results from our proper motion measurement are displayed in Fig. 10. We find evidence for non-zero proper motion directed toward the southeast, with a best-fit value of (18 +12 −13 , −35 +17 −18 ) mas yr −1 . 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.
Converting our 2D probability distribution for (µ α , µ δ ) to a distribution for the absolute value of proper motion, we obtain a 68% central credible interval µ * tot = 35 +16 −15 mas yr −1 . This corresponds to a transverse physical velocity of (570 ± 260) km s −1 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 −1 at a position angle Article number, page 11 of 30 A&A proofs: manuscript no. main µ α = 17.6 +12.8  The SNR G330.2+1.0 is located at a distance of at least 4.9 kpc , 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 ). Analysis of shock-velocity and filament expansion points toward a rather young age of the remnant, likely below 1000 years Borkowski et al. 2018). As part of their expansion analysis, Borkowski et al. (2018) noted that no evidence for "discernible motion" of the CCO CXOU J160103.1−513353 was found, which we attempt 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 south of (−2.7 +5.3 −5.4 , −6.4 +5.5 −5.4 ) mas yr −1 . 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.
The expected proper motion of an object co-rotating 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 −1 . 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 −1 . This corresponds to a µ α = − 2.7 +5.3  Table B.1 for constraints on all parameters. transverse kick velocity below 230 km s −1 , 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 travelled by the CCO during the "lifetime" of the 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.

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. H. E. S. S. Collaboration et al. 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 apparent to consist of multiple ring-like structures which consist of thin filaments on smaller scales ). 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 Cassam-Chenaï et al. 2004). The expansion of the outermost SNR filaments has been directly measured in three quadrants Tsuji & Uchiyama 2016;Tanaka et al. 2020), yielding projected shock velocities up to 3900 km s −1 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 & Aschenbach 1996) by around 4.5 toward south-west. Since the system is located quite nearby and the SNR age is likely rather  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. 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 −1 for an age of 2000 years. There are three archival Chandra observations of the CCO, 16 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 which can be convincingly cross-matched to a Gaia counterpart. 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 −1 × ∆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 −1 .
The resulting posterior distribution can be seen in Fig. 12. We measured a proper motion of (−4 +25 −24 , −20 ± 29) mas yr −1 , 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 (ap-plying only a simple translation of the coordinate system), we also obtained 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 −1 , corresponding to a limit on the 2D space velocity of around 230 km s −1 for a distance of 1 kpc. A future observation would certainly allow to greatly reduce these large error bars on the CCO's proper motion. 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. 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 Gaensler et al. 2008). Lovchinsky et al. (2011) performed spectral analysis of a deep ACIS-S observation of the SNR and from that estimated 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) predicted a projected velocity of 1400 − 2600 km s −1 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 −1 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. 17 They provided 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 detected comparatively slow motion of the CCO toward north, at significant uncertainty. Using our independent methods, we aim 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 obtained the posterior distribution shown in Fig. 13. We find that the CCO appears to be moving northward (with large uncertainty), at (−3 ± 8, 17 +10 −9 ) mas yr −1 . From this, we obtain a 68% central credible interval for the peculiar proper motion of µ * tot = 15 +10 −9 mas yr −1 , corresponding to a transverse velocity component of 320  Table B.1 for constraints on all parameters.
sources, found a proper motion of (−5, 14) mas yr −1 , 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. 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.
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 can be 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 ϑ which 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 −1 . 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. 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 −1 , corresponding to a projected shock wave velocity close to 6000 km s −1 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 measured the fastest expansion at a very similar location (labelled "B3" by them, finding τ = 590 yr) despite a different   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. 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. 18 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 is 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 to 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 18 Mathematically, the expansion at this evolutionary stage can ideally be described by the relation between radius and age R ∝ t β , with β, the deceleration parameter, expected to range between 0.4 and 1 for the Sedov-Taylor stage (Sedov 1959;Taylor 1950) and free expansion, respectively. The true age t is therefore related to the expansion age τ = R/Ṙ via t = βτ ≤ τ. 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 have clearly detected the motion of the larger, more diffuse emission features in the north and south of the SNR (labelled M to R). The expansion rates we found 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 (labelled "NNE" and "K" by them). In these regions, they inferred 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 think the most likely origin is 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 did 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.

Proper motion and neutron star kinematics
This work constitutes the first step toward the ambitious goal of building a sample of consistently measured transverse velocities for all central compact objects. 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 non-zero 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 were 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 non-zero proper motion. For RX J1713.7−3946, we were unable to obtain any significant signature of proper motion despite the small distance to the system. For G350.1−0.3, we did measure a mildly significant non-zero value, which is Article number, page 15 of 30 A&A proofs: manuscript no. main Notes. To provide a complete four-parameter astrometric solution, we display the best-fit positions of the CCO (α 0 , δ 0 ) at a given epoch t 0 , corresponding to the latest available observation of the respective target. The values for the projected physical velocity v proj are scaled to an assumed distance d , without the inclusion of any additional errors to account for uncertainties in d . The measurements of total proper motion µ tot and projected velocity v proj from this work have been corrected for the effect of Galactic rotation. All our measurements and errors correspond to the median and 68% central interval of the underlying probability distribution. All upper limits derived in this work are at 90% confidence. 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) measured very small proper motion, despite a striking offset of the CCO from the apparent SNR center. These examples illustrate that estimates of the neutron star 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 was investigated in detail for Tycho's SNR by Williams et al. (2013), who showed 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 neutron star 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 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 = µ tot d. This can be viewed as a lower limit to the physical velocity of the neutron star, and thus serves as a proxy to the violent kick experienced by the proto-neutron star at its birth. Therefore, by comparing the fastest measured neutron star 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 −1 , 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 −1 . A likely mechanism by which large neutron star kick velocities are achieved is the "gravitational tug-boat", in which an asymmetric ejecta distribution accelerates the proto-neutron star 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 −1 can be achieved in realistic explosion scenarios (Janka 2017), none of the presently 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 19 with one-dimensional σ = 265 km s −1 as given by Hobbs et al. (2005), the probability of measuring v proj ≤ 250 km s −1 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 neutron stars. 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.

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 backwards 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) 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 & Aschenbach 1996) to illustrate the characteristic structure of this SNR.

G15.9+0.2 and Kes 79: 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 have measured 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 which 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 have demonstrated 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 ex-tend 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 ∼ 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 instance, a ∼ 60 ks observation would effect an error reduction by a factor ∼ 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 verifying if its offset location within the SNR is indeed caused by a rather violent kick, or simply due to asymmetric expansion of the SNR. 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) measured rapid motion away from the CCO, at up to ∼ 9000 km s −1 (for a distance of 5 kpc) for most parts of the SNR shell. In contrast, they detected a far smaller rate of expansion for parts of the southwest 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 south-western region causing significant deceleration, after almost free expansion for the majority of the SNR's lifetime .

G350.1−0.3: A very young and highly peculiar SNR
Contrary to all other SNRs which host a CCO, the explosion site of G350.1−0.3 which was 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.
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 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.
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 −1 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 consis-tent with the 600 − 700 years found by Borkowski et al. (2020) and deduced in our analysis.
The overall multi-wavelength 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. However, 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. 05  agate more freely, consistent with the observed radio and harder X-ray emission. Borkowski et al. (2020) came 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 was 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 h 55 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 22 cm −2 ) 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 to G350.1−0.3 would have been observable with the naked eye at all. 5.2.4. RX J1713.7−3946: Evidence for non-uniform 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 & Aschenbach 1996).
The expansion of the SNR has been independently measured along three sectors of its shell, in the south-east, south-west and north-west, in three individual studies Tsuji & Uchiyama 2016;Tanaka et al. 2020). Interestingly, these works all showed that the maximal expansion speed of the shock-   , north-west (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). The bottom panel shows the derived likelihood and corresponding 90% confidence upper limits for the undecelerated expansion age τ, inferred for the respective filaments.
wave is remarkably uniform, on the order of 4000 km s −1 (or 800 mas yr −1 ), in all directions. 20 With our determination of the expansion center of the SNR, we can now 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 south-west (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 −1 , 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 north-west (Tsuji & Uchiyama 2016) and the south-east ) rims yields higher formal limits of τ < 2520 yr and τ < 3240 yr, respectively, as illustrated in Fig. 18. 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 signif-20 All three studies measured the motion of several clumps or filaments. However, here, we focus on the highest respective velocities, since these are expected to be closest to the undecelerated speed of the shock wave in that particular region.
Article number, page 19 of 30 A&A proofs: manuscript no. main icant non-zero motion of the CCO would bias our computation in any direction. Therefore, if we assume all input expansion velocities and errors to be reliable, we can infer that the shock wave must have been decelerated non-uniformly, leading to by far the largest present-day expansion rate inferred for the south-west 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 recent age estimates have been in the range 1500 − 2300 yr Tsuji & Uchiyama 2016), depending on the exact model assumed for ejecta and ISM density profiles. 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 south-west shell to be expanding almost freely.
However, it is somewhat puzzling why the shock wave would show the least decelerated expansion along the south-western 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 −1 (Okuno et al. 2018), around 1000 km s −1 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 south-east, whereas it has likely relatively recently struck the dense cavity wall in the west Fukui et al. 2012). Thus, one would naively expect to measure the least decelerated expansion of the shock wave in the south-east 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 south-west of the SNR by Tanaka et al. (2020), and is therefore to be taken with caution. In order to confirm these constraints, future follow-up 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.

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 have been performed (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 −1 , 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 neutron star 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 non-symmetric 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 error bars on the proper motion by a factor ∼ 2.5. The trajectory of the neutron star, 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 neutron star 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 neutron star 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 neutron star kick mechanism, in which the proto-neutron star experiences the "gravitational tug" of slow and massive ejecta during the first seconds of the explosion (Wongwathanarat et al. 2013).

Summary
In this work we have analyzed several sets of archival Chandra observations with the aim of deducing proper motion measurements consistently for all known central compact objects 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 Article number, page 20 of 30 by PSF morphology, aspect reconstruction uncertainties, or the effects of Galactic rotation.
In total, we have presented six new measurements of neutron star proper motion, four of which had not been directly targeted prior to this work. Thus, we have approximately doubled the existing sample of CCOs with quantitatively constrained kinematics. We have found moderately significant non-zero proper motion only for two objects in our sample, the neutron stars in Cas A and G350.1−0.3, with no evidence for any hypervelocity CCOs. In contrast, we have 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 −1 at distances of 5 kpc and 1 kpc, respectively.
Complementarily, we have 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 have applied our methods to comparatively shallow data sets to systematically search for signatures of expansion. For both objects, we found 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 have qualitatively confirmed, 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 −1 . 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 have directly confirmed 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 have 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 have shown 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 have combined several published expansion measurements of RX J1713.7−3946 with our estimate of the supernova explosion site. We found that the current expansion rate significantly differs between individual regions. In particular, we have shown that part of the south-western 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 have determined 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 neutron star 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 neutron star 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. Notes. The table shows the quantities underlying Figs. 6, 9 and 15. To ease comparison between the systems, we provide the resulting constraints on the expansion rate τ −1 for all three SNRs. In addition, we provide its inverse, the free expansion age τ, only for G350.1−0.3, since the expansion rate has to be significantly different from zero to obtain a sensible value for individual regions.