Three-dimensional Velocity Diagnostics to Constrain the Type Ia Origin of Tycho's Supernova Remnant

While various methods have been proposed to disentangle the progenitor system for Type Ia supernovae, their origin is still unclear. A circumstellar environment is key to distinguishing between the double-degenerate and single-degenerate (SD) scenarios since a dense wind cavity is expected only in the case of the SD system. We perform spatially resolved X-ray spectroscopy of Tycho’s supernova remnant (SNR) with XMM-Newton and reveal the three-dimensional velocity structure of the expanding shock-heated ejecta measured from Doppler-broadened lines of intermediate-mass elements. Obtained velocity profiles are fairly consistent with those expected from a uniformly expanding ejecta model near the center, whereas we discover a rapid deceleration (∼4000 to ∼1000 km s−1) near the edge of the remnant in almost every direction. The result strongly supports the presence of a dense wall entirely surrounding the remnant, which is confirmed also by our hydrodynamical simulation. We thus conclude that Tycho’s SNR is likely of SD origin. Our new method will be useful for understanding progenitor systems of Type Ia SNRs in the era of high-angular/energy-resolution X-ray astronomy with microcalorimeters.


INTRODUCTION
The origin of Type Ia supernovae (SNe) has not yet been clarified.While Type Ia SN explosion is known to be caused by a white dwarf (WD) in a binary system, two scenarios are under debate: the double-degenerate (DD; Iben & Tutukov 1984;Webbink 1984) or singledegenerate (SD; Whelan & Iben 1973).Since only SD systems require a stellar companion, a detection of the surviving star could be a smoking gun for this scenario.
Many previous optical surveys, however, failed to detect such stars in supernova remnants (SNRs) or preexplosion data so far (e.g., González Hernández et al. 2012).One of the well-observed SNRs in this context is Tycho's SNR (SN 1572).A promising candidate for a surviving star, namely Tycho G, was proposed by Ruiz-Lapuente et al. (2004).Followup observations however suggest that it is unlikely associated with SN 1572 (Kerzendorf et al. 2009;Xue & Schaefer 2015).Other nearby stars have also been claimed as alternative candidates, but they are still not decisive (Ihara et al. 2007;Kerzendorf et al. 2013Kerzendorf et al. , 2018)).
Recent theories and observations suggest that the companion stars might be too faint to detect (Hachisu et al. 2012;Li et al. 2011), which leads us to search for another method to distinguish the two scenarios.The surrounding environment provides key information; as found by Zhou et al. (2016), expanding molecular bubble exists around Tycho's SNR, which is likely to be evidence for a SD progenitor.Indirect evidence for the SD scenario also comes from a proper motion study of Tycho's SNR with Chandra (Tanaka et al. 2021).They discovered a drastic deceleration of the blast waves in the southwest during the last ∼ 15 yr observations, suggesting that a dense wall is surrounding a low-density cavity and this environmental structure was formed by a mass-accreting white dwarf before the explosion (see also, Kobashi et al. 2023).
If the cavity wall exists around Tycho's SNR and is/was hit by the blast waves in all directions, another piece of supporting evidence may be found in the information along the line of sight.It is generally known that when an SNR expands into a dense material, a "reflection (reflected) shock" begins to move backward behind a contact discontinuity, as is shown by a simple fluid dynamical calculation (e.g., Hester et al. 1994).Several X-ray observations of proper motions suggest that backward-moving ejecta or filaments found in Cassiopeia A (Sato et al. 2018) and RCW 86 (Suzuki et al. 2022) are likely reflection shocks produced by ambient molecular clouds or cavity walls.On the other hand, a line-of-sight velocity structure is, in general, difficult to detect due to a lack of the energy resolution of current X-ray detectors.
In the case of Tycho's SNR, line-of-sight velocities were partially measured (Sato & Hughes 2017;Millard et al. 2022) and global expansion structure was reported (Furuzawa et al. 2009;Hayato et al. 2010).A more comprehensive investigation using Chandra data is recently reported by Godinaud et al. (2023).In this paper, in order to assess the presence of the cavity wall, we aim to reveal a global velocity structure of the expanding ejecta of Tycho's SNR using the large dataset of XMM-Newton.Throughout the paper, the distance of Tycho's SNR is assumed to be 3 kpc (cf.Williams et al. 2016).Errors of parameters are defined as 1σ confidence intervals.

OBSERVATIONS
For the following analysis, we used XMM-Newton data of Tycho's SNR taken for the instrument calibration during 2005-2009 (Table 1).Note that we did not use Chandra data, while a better angular resolution might be more suitable for our study.This is because we found that ACIS data are suffering from an effect of "sacrificial charge" caused by a bright source (Grant  et al. 2003;Kasuga 2022), which fills electron traps along readout paths of the CCD and changes charge transfer inefficiency; this effect results in an overestimation of photon energies.Among the instruments aboard XMM-Newton, we used MOS data for our analysis because of the better energy resolution than pn.The event files were processed using the pipeline tool emchain in version 19.1.0 of the XMM-Newton Science Analysis System (SAS; Gabriel et al. 2004).In all the observations, the whole remnant is within a single CCD (namely CCD1) equipped at the center of the chip array of the MOS.We used blank-sky data created from background regions of several point-source observations (Table 1); we eliminated each central source region and combined all the remaining CCD1 data.

ANALYSIS AND RESULTS
Figure 1 shows the background-subtracted MOS spectrum of Tycho's SNR.We detected prominent emission lines of intermediate-mass elements (IMEs; Si, S, Ar, and Ca), and iron-group elements (IGEs; especially Cr and Fe).Since this remnant is one of the ejectadominated SNRs (Miceli et al. 2015), these heavy elements including O, Ne, and Mg (hereafter, ONeMg) are all from shock-heated ejecta material.In Figure 1, we also display a non-X-ray background (NXB) spectrum, whose count rate is lower than ∼ 1% of the SNR flux almost in the entire energy band and its uncertainty is negligible for our analysis.The following spectral fit was done on the platform of the Xray SPECtral fitting package (XSPEC) 12.11.1 (Arnaud 1996) with the Atomic DataBase (AtomDB) 3.0.9(Foster 2015).A fitting method we used was W -statistics  (Wachter et al. 1979), which is a generalized method of C -statistics (Cash 1979) with background subtraction (see also Kaastra 2017).
In order to investigate the spatial distribution of physical properties of Tycho's SNR, we divided the entire region into square grids with a spacing of 15 ′′ (totally 34×34 = 1156 segments) as presented in Figure 2: in total 845 segments (grid regions) covering the remnant are available for the spectral analysis.For each grid region, we performed a wide-band spectroscopy between 0.5 keV and 8 keV.Spectral and detector response files were generated with mos spectra in SAS.Although time variations are reported by previous studies with Chandra (cf. Katsuda et al. 2010;Tanaka et al. 2021), a typical expected angular variation is ∼ 0.14 ′′ yr −1 , which is much smaller than the grid size and thus negligible for our analysis.We thus simply combined all the seven observations during 2005-2009 for better statistics for every single grid.We adopted marfrmf, and addrmf, which are available in the software package of the High Energy Astrophysics Software (HEASoft) 6.28 (Blackburn 1995), for combining ancillary response (arf) and response matrix files (rmf).
On the basis of the previous X-ray spectroscopic analysis of Tycho's SNR (e.g., Tamagawa et al. 2009), we adopted a multi-component Non-Equilibrium Ionization (NEI) model plus a power-law continuum; the later component is representing the synchrotron emission (Fink et al. 1994).Two Gaussians are also included at ∼ 0.7 keV and ∼ 1.2 keV to compensate for the known uncertainty of the plasma code (see discussion in Okon et al. 2020).The absorption column density N H using the Tübingen-Boulder model (Wilms et al. 2000) was fixed at 7.5×10 21 cm −2 , derived from an average of previous measurements (e.g., Yamaguchi et al. 2017;Okuno et al. 2020).
Given that the ejecta consists of pure metals stratified into layers in the shell of the remnant (Hayato et al. 2010), we applied several NEI components for a thermal emission corresponding to different element groups, namely ONeMg, IME, and IGE.For the ONeMg component, O is fixed to 10 5 whereas Ne and Mg are free; we confirmed that these abundances do not affect the following analysis and results since the ONeMg component is almost negligible in the middle-energy band that we focus on.For the IME component, Si is fixed to 10 5 whereas S, Ar, and Ca are free.The IGE component is further divided into two components (IGE1 and IGE2) with different plasma temperatures.For these components, Fe is fixed to 10 5 and linked the abundances between IGE1 and IGE2.We tied the electron temperature kT e among the ONeMg, IME and IGE1 components and varied only the ionization timescale τ for simplicity under an assumption that in each grid the temperature gradient is negligible along the line of sight.The IGE2 component has a different set of kT e and τ from the IGE1 component since the ionization state of Fe derived from the line centroid of Fe Kα seems different from that estimated from the Fe-L complex.The thermal components required for our analysis are then a combination of the low-kT e ONeMg, IME, and IGE1 plus the high-kT e IGE2 (hereafter, multi-NEI component).We however found that the above model (with a power law) still cannot reproduce spectra well, particularly for regions near the center.This seems reasonable because the ejecta of Tycho's SNR has a velocity structure as reported by several studies (e.g., Hayato et al. 2010;Sato & Hughes 2017).We thus applied the two convolution models (vmshift and gsmooth in XSPEC) to all the components for representing a Doppler shift and line broadening effect of each element group (hereafter, we call this "modified multi-NEI model").Note that we set the index parameter in gsmooth to be 1 so that the line width is proportional to the photon energy.

Parameter distributions
In Figure 3, we show spectra and best-fit results of example regions.All the spectra are well reproduced by the modified multi-NEI model.Parameters of the IGE1 and IGE2 components are especially well constrained even in regions near the outer rim where the nonthermal emission is dominant (for instance, region (3, 20) displayed in Figure 3).The results enable us to investigate the spatial distributions of the parameters for all the thermal components, as displayed in Figures 4 and  5.

Electron Temperature and Ionization Timescale
The electron temperature distributions shown in Figure 4 indicate that kT e of the ejecta is nearly uniform (∼ 2 keV) whereas that of the IGE2 (i.e., hot component of Fe) is significantly higher (> 7 keV) in the southeastern knot (Fe knot; Yamaguchi et al. 2017) than those in the other regions.These results imply that the origin of the IGE2 is different from the low-kT e components, as is pointed out by Yamaguchi et al. (2017).The overall spatial trend of each component is consistent with those shown by previous studies (e.g., Miceli et al. 2015;Matsuda et al. 2022).The ionization timescale τ for IME shows a relatively uniform distribution, whereas we found that τ for the ONeMg is somewhat higher in the northwest (> 5 × 10 10 cm −3 s) than in the other regions (∼ 2 × 10 10 cm −3 s).Since the higher-τ regions roughly coincide with a location of the RS determined by Warren et al. (2005), it might be interpreted as evidence that these elements in the northwest were heated by the revers shock (RS) earlier than in the other directions.We found that τ for the IGE1/2 has large uncertainties near the rim.This is because the synchrotron emission is dominant in these regions, which makes it difficult to estimate τ from line intensity ratios.Note that these uncertainties are not significant for the following discussion because our main purpose is to determine the kinematic structure of the IME component.

Three-Dimensional Velocity Diagnostics
Hereafter we primarily focus on the ejecta distribution of the IMEs shown in Figure 5, since systematic uncertainties of the Doppler shift and the line width are much smaller in comparison with the other components.The Doppler velocity map of the IMEs shows a clumpy distribution, in which blueshift/redshift regions are somewhat clustering with each other.The trend is consistent with the mean photon energy maps of Si Heα shown by Sato & Hughes (2017) and Millard et al. (2022).A recent comprehensive analysis with Chandra also supports our result (Godinaud et al. 2023).The maximum velocity of the IME ejecta is roughly ∼ ±2, 000 km s −1 , which is consistent with the previous measurements for the Si ejecta.The result is understandable because in the IME component the most prominent line is Si Heα.
The Doppler velocity map for the IGEs is weakly correlated with that for the IME: the northwestern and southeastern parts (around the Fe knot) are blueshifted, whereas several clumps in the south are redshifted.The correlation implies that the global ejecta structure for the IMEs spatially overlaps with that for the IGEs, especially Fe.On the other hand, the ONeMg component shows no clear correlation with the heavier elements.This may confirm that these light elements spatially distribute above the core elements such as Fe and Si as suggested also by the ionization timescale distribution of the ONeMg component (see Section 3.1.1).
From the line width maps, which are the main indicators of the ejecta kinematics, we clearly see that the IMEs show an azimuthally symmetric velocity structure.A similar trend is apparent in the line width map of the IGEs whereas that of the ONeMg shows a different asymmetric distribution.Since the lighter elements are distributed in outer regions, this component is likely affected by an inhomogeneous ambient density structure.The heavier elements, especially the IMEs, are expanding symmetrically from a single point (probably from the center of the remnant), as previously reported (Furuzawa et al. 2009;Hayato et al. 2010;Sato & Hughes 2017).The result will be helpful for constraining the 3D expanding structure of Tycho's SNR with a line-of-sight velocity of the ejecta.
We constructed azimuthally averaged radial profiles of the obtained line widths for each octant, as displayed in Figure 6.The resultant profiles are similar: line widths (σ) are highest near the center and gradually decrease toward the rim of the remnant, which is reasonable given that Tycho's SNR is approximately spherically symmetric even along the ling of sight.Only the profile for θ = 90 • -135 • deviates from this trend, which can be partially attributed to thermal broadening in the south-  eastern knot reported by Williams et al. (2020).We note that thermal broadening is also expected from the contact discontinuity to the reverse shock front.According to Badenes et al. (2006), this effect is significant only in a narrow region (less than 10% of the line-of-sight volume) and is negligible beyond 200 arcsec from the center.
From Figure 6, we found that σ falls to ∼ 0 eV far inside the shell.This trend is clearly exemplified in the profiles for θ = 180 • -225 • and 315 • -360 • .It may be somewhat counterintuitive because if we assume a normal radial velocity variation of expanding material, σ should reach 0 eV at the outer edge of the profile.As reported by Godinaud et al. (2023), proper motion velocities (V xy ) near the rim of the remnant exceed ∼ 3, 000-4,000 km s −1 .On the other hand, the line width σ we measured represents an average of the superposition of the velocity components in the line-of-sight direction.These results may be consistently interpreted by assuming that the ejecta has a double-shell velocity structure: an inner thick high-velocity component overlaid with a thin lower-velocity shell.This assumption is probably reasonable because Chandra HETG observations (Millard et al. 2022) suggest lower line-of-sight velocities near the shell.Given that the trend is common across all the octants, we infer that our result may reflect a global velocity structure of Tycho's SNR.  6 implies that some deceleration may have occurred near the edge of Tycho's SNR.This scenario is possible if we consider a recent interaction between the blast waves and ambient dense gas.Since previous studies (e.g., Zhou et al. 2016;Tanaka et al. 2021) suggest the presence of such dense material around the remnant, it is reasonable to consider that the forward shock recently hit the cavity wall essentially in all the directions, which can affect the global radial velocity structure.In order to clarify the global structure of Tycho's SNR, we measured the line-of-sight expansion velocity v ∥ of the IME in each segment by refitting the spectra with a Doppler-broadening model, with reference to Hayato et al. (2010).We focused on the data between 1.6 keV and 6.0 keV because the IME component is generally dominant in the middle-energy band.For simplicity, we assume that the observed line broadening can be explained by a two-NEI (plus a power law) model with different velocities v r = v ∥ and v b = −v ∥ , representing the red/blue-shifted (i.e., moving backward and forward) components, respectively.The parameters of the two NEI components, except for the normalization, were linked to each other and fixed at the best-fit values of the multi-NEI model.We also included a common offset parameterized as v 0 for representing a bulk motion of Tycho's SNR and/or a systematic uncertainty in the instruments.By averaging the line shift at the outermost edges of the remnant, we obtained v 0 = −420 km sec −1 .A local red/blue shift is in this case explained by an intensity ratio of v r and v b components.We confirmed that all the spectra were well fitted with the above model.An example fit is shown in Figure 7.
Figure 8 displays the v ∥ profiles for different directions based on the best-fit results with the Dopplerbroadening model.Since v ∥ is measured independently from v 0 and directly reflects the expansion structure, the uncertainty in the latter does not critically affect the following discussion.Incidentally we found that a similar analysis with Chandra gives a much broader line shape and a stronger blue shift in almost all the regions.A possible cause of this discrepancy is the effect of "sacrificial charge" mentioned in Section 2. While it is out of the scope of this paper to argue this overestimation and calibration issue concerning Chandra, we exemplify spectral comparison between two instruments in Apendix A that will be useful for future studies.
As indicated in Figure 8, the profile shapes of v ∥ are similar to those of the line widths (Figure 6).We then constructed a model profile for v ∥ under the following assumption.Since the reverse shock of this remnant has not reached the center yet, the X-ray emitting region is within r RS ≤ r < r ejecta , where r RS and r ejecta are the positions of the reverse shock and the outer boundary of the ejecta, respectively.Assuming a uniform and lowdensity environment around Tycho's SNR, the expan- sion velocity of the ejecta (hereafter, v IME ) within this region can be considered as roughly being constant with radius (e.g., Blondin & Ellison 2001).As explained in Appendix 5 and Figure 9, we can describe v ∥ (x ′ ), which is an average velocity along the line-of-sight at a position x = x ′ , as follows: where, (2) We fitted the model v ∥ (x ′ ) described by equation (1) to the v ∥ profiles obtained by the Doppler model applied to the XMM-Newton data as shown in Figure 8. From the best-fit results, we found a significant discrepancy between the model and the data except for the southeast direction (θ = 90 • -135 • , 135 • -180 • ) where the Feknot is located (Yamaguchi et al. 2017).In Figure 8, we also show the locations of the contact discontinuity r CD , which is estimated as the average radius of the outer boundary of the ejecta in each sector.If the ejecta is expanding in a uniform and low-density environment without deceleration, r ejecta should coincide with the contact discontinuity (r CD ), i.e., r ejecta = r CD .This assumption, however, fails to reproduce the observed v ∥ profiles, leading us to favor a non-uniform expanding structure of the ejecta.

Velocity Structure of Tycho's SNR
From the best-fit v ∥ profiles (Figure 8), we obtained v IME , r RS , and r ejecta as the fitting parameters, for each direction as summarized in Table 2.The averaged ejecta velocity v IME ∼ 3900 km s −1 is between the predicted for a distance of 3 kpc).We also found that v IME at θ = 180 • -270 • (southwest region) is larger than those at the other directions.It is roughly aligned with previous studies suggesting encountering dense gas toward the northeast direction in past (Lee et al. 2004;Zhou et al. 2016;Williams et al. 2016;Godinaud et al. 2023), and more recent deceleration of blast waves in the southwest (Tanaka et al. 2021).The average position of the reverse shock r RS = 175 ′′ ± 3 ′′ is also fairly consistent with the previous estimations, i.e., 183 ′′ (Warren et al. 2005) or 158 ′′ (Yamaguchi et al. 2014).These results indicate that the velocity structure and the ejecta morphology are naturally explained by the standard SNR evolution model (e.g., Truelove & McKee 1999) except for the outer layers.
In order to understand the cause of the discrepancy between the v ∥ (x) model and the v ∥ profiles near the rim as seen in Figure 8, we converted the measured v ∥ to v IME by solving an inverse problem of equation ( 1).Here, we define v IME as an "expected" averaged velocity of the ejecta along the line of sight.Using equation (1), if the ejecta is not decelerated and has a constant velocity all over the emitting region, the value v IME is expected to be constant throughout the remnant interior (i.e., v IME (x ′ ) = constant, irrespective of x ′ ).If the ejecta is somehow decelerated, v IME should become lower and not follow the constant trend.Therefore v IME is not equivalent to the velocity itself but should be considered as an indicator of the deceleration, i.e., "deviation" from expected value assuming only a high-velocity expansion.
With this property in mind, we present profiles of calculated v IME in Figure 10.Interestingly, v IME keeps constant within the radius of ∼ 200 ′′ above which it significantly drops from ∼ 4000 km s −1 to ∼ 1000 km s −1 in almost every direction.Although the parameter v IME does not directly mean a decelerated velocity as noted above, it is obvious that the velocity structure evidently contradict the expectation from the uniform expansion.We thus conclude that the ejecta was decelerated in recent past by dense gas in the vicinity of Tycho's SNR.
When an expanding SNR shock collides with a dense material, numerical simulations suggest that a "reflection (reflected) shock" traverses back into the ejecta (cf.Dwarkadas 2005;Orlando et al. 2022).Several lines of observational evidence are also reported from corecollapse SNRs interacting with molecular clouds (e.g., Inoue et al. 2012;Sato et al. 2018).In the case of Tycho's SNR, as the previous works suggest (e.g., Tanaka et al. 2021;Kobashi et al. 2023), it is more reasonable to expect a cavity wall formed by the progenitor system before the SN explosion.Assuming that the cavity had a radius of 4.5 pc with a density of 0.3 cm −3 (see also Zhou et al. 2016), we present a hydrodynamical calculation of shock propagation as shown in Figure 11.In our setup, a blast wave generated by a normal Type Ia SN is colliding with a wall with a number density of 100 cm −3 , which is consistent with the estimate given by Tanaka et al. (2021).The result demonstrates that the blast wave hits the wall and is decelerated rapidly, which gives a discontinuous low-velocity structure near the edge.We confirmed that the velocity of the decelerated region is typically ∼ 1000 km s −1 during ∼ 20 yr after the collision.
As our simulation predicts, the generated reflection shock will move backward and converge onto the reverse shock within a few tens of years.It provides important constraints on when (and where) Tycho's SNR reached the cavity wall.Notably this timescale is fairly in agreement with the previous proper-motion measurement (Tanaka et al. 2021), in which the authors observed a shock deceleration during last ∼ 15 yr in the southwest.We also infer that the cavity was nearly spherically symmetric around the SN explosion.We further found that the reflection shock forms a positive velocity gradient to the outside (R = 4.5-4.7 pc in Figure 11); the trend is suggestively similar to the v IME profiles shown in Figure 10 (for example, θ = 180 • -225 • ).We thus conclude that our result provides further evidence of the cavity formed by a wind from the white dwarf during mass accretion in the SD progenitor system (Hachisu et al. 1996).Although there remains large uncertainties in the measurement with the CCDs, future observations with higher energy resolution will provide a more detailed velocity structure, which hints at the environment of Tycho's SNR and its progenitor system.Our method with high angular resolution microcalorimeters such as Athena (Barret et al. 2018) and Lynx (Gaskin et al. 2019) will provide a new way to distinguish between the SD and DD progenitor scenarios for Type Ia SNRs.

CONCLUSIONS
We revealed a global three-dimensional velocity structure of Tycho's SNR using XMM-Newton data, and investigated the presence of a wind bubble which was suggested by a recent proper-motion observation (Tanaka et al. 2021).We performed a spatially resolved spectroscopy of the remnant (845 grid cells in total) and found that emission lines of the ejecta such as Si Heα can be explained by a blue/red-shifted thermal model in almost all regions.While the obtained velocity structure approximately indicates a spherically symmetric expansion, the radial profile drops down far inside the edge in almost all directions, which implies a significant deceleration of blast waves around the remnant.
We converted the obtained line-of-sight velocities v ∥ to expansion velocities v IME and found possible evidence that ejecta just behind the CD was slowed down to v IME ∼ 1000 km s −1 , whereas global expansion velocity is v IME = 3870 ± 40 km s −1 .We also confirmed that average radii of RS and ejecta are r RS = 175 ′′ ± 3 ′′ and r ejecta = 210 ′′ ± 1 ′′ , respectively, which are in good agreement with other estimations.Our result strongly suggests a recent collision of Tycho's SNR with dense surrounding materials, supporting the presence of the cavity wall (wind bubble).
We performed a hydrodynamical calculation of shock propagation into a dense material.The result indicates that a reflection shock is formed behind the forward shock and causes a significant deceleration, which can consistently explain the observed velocity structure.We thus conclude that our result supports the presence of the cavity and a dense wall, and therefore the SD scenario is preferable for the origin of Tycho's SNR.Since our study is based on the CCD data, a further analysis with a better energy resolution with microcalorimeters will reveal more detailed velocity structure of this remnant.Also note that the method we present here is applicable to other Type Ia SNRs if we use future observatories like Athena and Lynx.In addition to ways to test explosion models from SNR morphology (Ferrand et al. 2021) or abundances (Yamaguchi et al. 2015), the velocity structure will provide useful information on the presence of stellar wind bubbles and hence will be helpful to identify the SD or DD progenitor scenario.T.K. appreciates all supports for his short staying at Kyoto by Takeshi Go Tsuru, Toshihiro Fujii, and their group members at Kyoto University, and Hein Mallee and Takako Okamoto at the Research Institute for Humanity and Nature.The authors are very grateful to Paul Plucinsky, Hirokazu Odaka, Shigehiro Nagataki, and Gilles Ferrand for giving us precious advice on this study.The authors also thank Tsuyoshi Inoue for the helpful discussion on the shock physics.This work is mainly supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Nos.JP21J11443 (T.K.), JP22H01265 (H.U.), JP20H00174 (K.M.), JP19H01936 (T.T.), JP21H04493 (T.T.), and JP19K03908 (A.B.).

A. SACRIFICIAL CHARGE ISSUE ENCOUNTERED IN ANALYSIS OF CHANDRA DATA
In this appendix, we briefly present an effect of "sacrificial charge" that we encountered in our analysis of Chandra ACIS data.This issue will become severe when a bright diffuse source is observed with a CCD detector.Since a signal charge is lost during readout with a certain probability due to traps of the lattice of silicon, an enhancement of charge transfer inefficiency (CTI) cannot be generally ignored for in-orbit photon-counting CCDs (Townsley et al. 2000).In order to recover energy losses, the CTI correction method is therefore applied in a standard reduction and analysis of Chandra ACIS dataset (Plucinsky et al. 2002).However, if a very high-count-rate X-ray object is observed with the ACIS CCDs, charge deposited along the readout path fills the traps during readout and thus the CTI is unintentionally and excessively corrected.This is called sacrificial charge, which gives a wrong energy peak and line width (Grant et al. 2003).The sacrificial charge issue would become more considerable for high-angular resolution instruments.
Figure 12 shows a comparison between spectra of Tycho's SNR, for which the same region is extracted with XMM-Newton and Chandra.We found a clear difference between the two spectra and confirmed that the ACIS spectrum shows a broader and more blue-shifted feature, even if taking into account the difference in energy resolution between each other.Most of the ACIS spectra are suffering from this problem; estimated expansion velocities exceed 3000 km s −1 even around the edge of the remnant.These results are clearly attributed to the effect of the sacrificial charge explained above.We therefore did not use the Chandra data for our analysis.

B. MODELING OF LINE-OF-SIGHT VELOCITY STRUCTURE
As shown in Figure 9, we assume that Tycho's SNR is spherically symmetric (Hayato et al. 2010) and the X-ray emitting ejecta (IME) is filled between the RS and the CD with an expansion velocity v IME .In a two-dimensional plane, which is a cut surface of the angular direction of the remnant, we define the x and z axes as the angular and the line-of-sight directions, respectively.In any position x = x ′ and z = z ′ (r ′ ≡ √ x ′2 + z ′2 ), the line-of-sight velocity is described as v(x ′ ; z ′ ) = v IME cos φ, where φ is the angle between the direction of v IME and the z axis: tanφ = x ′ /z ′ .
Here, an integrated line-of-sight velocity f (x ′ ; R) at x = x ′ is calculated as where R is the radius of the circumference of the ejecta on the two-dimensional plane.By defining Θ ≡ arcsin (x ′ /R), we can modify this equation as  (B5) The observed line-of-sight velocity F 1 (x ′ ; R) should be the density-averaged value of f (x ′ ; R).If the ejecta density is uniform, F 1 (x ′ ; R) is calculated as On the other hand, if the ejecta is filled between R = R 1 and R = R 2 (R 1 < R 2 ), the line-of-sight velocity F 2 (x ′ ; R) can be described as

Figure 2 .
Figure 2. Chandra ACIS image (734 ks in total) of Tycho's SNR (1.80-1.92keV; Si Kα line band) overlaid with the grids we used for our analysis.The definitions of the region number (X, Y ) and the azimuthal angle θ (from the north in counterclockwise) are also displayed.The square regions enclosed by the thick lines are those from which we obtained example spectra shown in Figure 3.

Figure 3 .
Figure 3. X-ray spectra of the example regions shown in Figure 2. The black points are binned observation data and the green solid line is the best-fit model.The blue and red solid lines represent the ONeMg and IME components, respectively.The solid and dashed orange lines indicate the IGE1 and IGE2 components, respectively.The magenta dotted line is the synchrotron emission component.The black dotted lines represent the Gaussians at 0.7 keV and 1.2 keV.The residuals are displayed in each lower panel.

Figure 4 .
Figure 4. Best-fit parameter maps of the electron temperatures kTe and the ionization timescale τ for each component: low-kTe ONeMg, IME, and IGE1 plus high-kTe IGE2 (see text).
Measurement of the ejecta velocityAs claimed in Section 3.1.2,the result shown inFigure

Figure 5 .
Figure 5. Doppler shift (left) and line width (right) maps of Tycho's SNR.Results estimated from the ONeMg, IME, and IGE components are shown in the top, middle, and bottom panels, respectively.Numbers around the maps represent the xy coordinate of the region (X, Y ).The size of each dot represents a typical size of the 1σ error bars; smaller dots are more uncertain.Reference values of the errors are shown in the bottom-left corner of each panel, where the five dots from the smallest to largest indicate ×3, ×2, ×1, ×1/2, and ×1/3 of the reference values.The regions where no colored dot is displayed indicate that the obtained parameter has a large uncertainty and the error exceeds ×3 of the reference values.

Figure 6 .Figure 7 .
Figure6.Radial profiles of the line width from the center to the outside for each octant.The x-axis corresponds to the projection position x defined in the text and Figure9region(22, 15)

Figure 8 .Figure 9 .
Figure 8. Radial profiles of v ∥ from the center to the outside for each octant.Best-fit results of equation (1) is shown in green.The vertical dotted lines indicate positions of the contact discontinuity rCD.

Figure 11 .
Figure 11.Results of a hydrodynamical simulation.Top and middle panels display the density and velocity profiles with radius, respectively.The bottom panel shows the zoomin view of the velocity profile near the outer edge, where the reflection shock is generated.

Figure 12 .
Figure 12.Spectra of a region in Tycho's SNR obtained with the XMM-Newton MOS (red) and the Chandra ACIS (black).The top left and right panels show the same spectra but x-axis is in the logarithmic and linear scales, respectively.The bottom panels show the ACIS (left) and MOS (right) spectra displayed with the black data points.The best-fit model obtained with the MOS data is overlaid onto each spectrum; the black, red, and green lines indicate the power-law, thermal, and total model shapes, respectively.