Three-dimensional morphological asymmetries in the ejecta of Cassiopeia A using a component separation method in X-rays

Recent simulations have shown that asymmetries in the ejecta distribution of supernova remnants can still reﬂect asymmetries from the initial supernova explosion. Thus, their study provides a great means to test and constrain model predictions in relation to the distributions of heavy elements or the neutron star kicks, both of which are key to better understanding the explosion mechanisms in core-collapse supernovae. The use of a novel blind source separation method applied to the megasecond X-ray observations of the well-known Cassiopeia A supernova remnant has revealed maps of the distribution of the ejecta endowed with an unprecedented level of detail and clearly separated from continuum emission. Our method also provides a three-dimensional view of the ejecta by disentangling the red- and blue-shifted spectral components and associated images of the Si, S, Ar, Ca and Fe, providing insights into the morphology of the ejecta distribution in Cassiopeia A. These mappings allow us to thoroughly investigate the asymmetries in the heavy elements distribution and probe simulation predictions about the neutron star kicks and the relative asymmetries between the different elements. We ﬁnd in our study that most of the ejecta X-ray ﬂux stems from the red-shifted component, suggesting an asymmetry in the explosion. In addition, the red-shifted ejecta can physically be described as a broad, relatively symmetric plume, whereas the blue-shifted ejecta is more similar to a dense knot. The neutron star also moves directly opposite to the red-shifted parts of the ejecta similar to what is seen with 44 Ti. Regarding the morphological asymmetries, it appears that heavier elements have more asymmetrical distributions, which conﬁrms predictions made by simulations. This study is a showcase of the capacities of new analysis methods to revisit archival observations to fully exploit their scientiﬁc content.


Introduction
Cassiopeia A (hereafter, Cas A) is among the most studied astronomical objects in X-rays and is arguably the best-studied supernova remnant (SNR). Investigation of the distribution of metals on sub-parsec scales is possible because it is the youngest core-collapse (CC) SNR in the Milky Way (about 340 yr old; Thorstensen et al. 2001), its X-ray emission is dominated by the ejecta metals (Hwang & Laming 2012), and it is relatively close (3.4 kpc, see Reed et al. 1995;Alarie et al. 2014). Cas A benefits from extensive observations (about 3 Ms in total by Chandra), making it an ideal laboratory to probe simulation predictions regarding the distribution of ejecta metals.
In the last few years, three-dimensional simulations of CC supernovae (SNe) have begun to produce testable predictions of SNe explosion and compact object properties in models using the neutrino-driven mechanism (see reviews by Janka et al. 2016;Müller 2016). In particular, explosion-generated ejecta asymmetries (Wongwathanarat et al. 2013;Summa et al. 2018;Janka 2017) and neutron star (NS) kick velocities (DeLaney & Satterfield 2013) appear to be key elements in CC SN simulations that Cas A's data can constrain. Although it is challenging to disentangle the asymmetries produced by the surrounding medium from those inherent to the explosion, Orlando et al. (2016) have explored the evolution of the asymmetries in Cas A using simulations beginning from the immediate aftermath of the SN and including the three-dimensional interactions of the remnant with the interstellar medium. Similar simulations presenting the evolution of a Type Ia SNR over a period spanning from one year after the explosion to several centuries afterward have been made by Ferrand et al. (2019), showing that asymmetries present in the original SN can still be observed after centuries. The same may go for the CC SNR Cas A, and a better knowledge of its three-dimensional morphology could lead to a better understanding of the explosion mechanisms by providing a way to test the simulations.
An accurate mapping of the different elements' distributions, the quantification of their relative asymmetries, and their relation to the NS motion would, for example, allow us to probe the simulation predictions that heavier elements are ejected more asymmetrically and more directly opposed to the NS motion than lighter elements (Wongwathanarat et al. 2013;Janka 2017;Gessner & Janka 2018;Müller et al. 2019). On this topic, this paper can be viewed as a follow-up to Holland-Ashford et al. (2020), a study that aimed to quantitatively compare the relative asymmetries of different elements within Cas A, but which was hindered due to difficulties in separating and limiting contamination in the elements' distribution. Moreover, in that analysis, the separation of the blue-and red-shifted parts in these distributions was not possible.
Here, we intend to fix these issues by using a new method to retrieve accurate maps for each element's distribution, allowing us to further investigate their individual and relative physical properties. This method is based on the General Morphological Components Analysis (GMCA, see Bobin et al. 2015), a blind source separation (BSS) algorithm that was introduced for X-ray observations by Picquenot et al. (2019). It can disentangle both spectrally and spatially mixed components from an X-ray data cube of the form (x, y, E) with a precision unprecedented in this field. The new images thus obtained suffer less contamination by other components, including the synchrotron emission. It also offers the opportunity to separate the blue-and red-shifted parts of the elements' distribution, thereby facilitating a threedimensional mapping of the X-ray emitting metals and a comparison of their relative asymmetries. Specifically, the GMCA is able to disentangle detailed maps of a red-and a blue-shifted parts in the distributions of Si, S, Ca, Ar, and Fe, thus providing new and crucial information about the three-dimensional morphology of Cas A. This is a step forward as previous studies intending to map the distribution of the individual elements and study their asymmetries in Cas A in X-rays (Hwang & Laming 2012;Katsuda et al. 2018;Holland-Ashford et al. 2020) were not able to separate red-and blue-shifted components.
This paper is structured as follows. In Sect. 2, we will describe the nature of the data we use (Sect. 2.1), our extraction method (Sect. 2.2), our way to quantify the asymmetries (Sect. 2.3), and our method to retrieve error bars (Sect. 2.4). In Sect. 2, we will present the images resulting from the application of our extraction method (Sects. 3.1 and 3.2), and we will discuss the interpretation of the retrieved images as blue-or redshifted by looking at their associated spectra (Sect. 3.3), and will present the results of a spectral analysis on these same spectra (Sect. 3.4). Lastly, we will discuss in Sect. 4 the physical information we can infer from our results. Section 4.1 will be dedicated to the interpretation of the spatial asymmetries of each line emission, while Sects. 4.2 and 4.4 will focus respectively on the mean direction of each line's emission and on the NS velocity. A comparison with the NuSTAR data of 44 Ti will finally be presented in Sect. 4.5.

Nature of the data
Spectro-imaging instruments, such as those aboard the current generation of X-ray satellites XMM-Newton and Chandra, provide data comprising spatial and spectral information: The detectors record the position (x, y) and energy E event by event, thereby producing a data cube with two spatial dimensions and one spectral dimension. For our study, we used Chandra observations of the Cas A SNR, which was observed with the ACIS-S instrument in 2004 for a total of 980 ks (Hwang et al. 2004). We used only the 2004 data set to avoid the need to correct for proper motion across epochs. The event lists from all observations were merged in a single data cube. The spatial (of 2 ) and spectral binning (of 14.6 eV) were adapted so as to obtain a sufficient number of counts in each cube element. No background subtraction or vignetting correction has been applied to the data.

Image extraction
The main concept of GMCA is to take into account the morphological particularities of each component in the wavelet domain to disentangle them, without any prior instrumental or physical information. Apart from the (x, y, E) data cube, the only input needed is the number n of components to retrieve, which is userdefined. The outputs are then a set of n images associated with n spectra. Each couple image-spectrum represents a component: The algorithm makes the assumption that every component can be described as the product of an image with a spectrum. Thus, the retrieved components are approximations of the actual components with the same spectrum on each point of the image. Nevertheless, Picquenot et al. (2019) showed that when tested on Cas A-like toy models, the GMCA was able to extract morphologically and spectrally accurate results. The tested spectral toy models included power-laws, thermal plasmas, and Gaussian lines. In particular, in one of these toy models, the method was able to separate three components: two nearby partially overlapping Gaussian emission lines and power-law emission. The energy centroids of both Gaussians were accurately retrieved, despite their closeness. Such a disentangling of mixed components with similar neighboring spectra cannot be obtained through line-interpolation, and fitting of a two-Gaussian model region by region is often time consuming, producing images contaminated by other components with unstable fitting results.
In the same paper, the first applications on real data of Cas A were promising, in particular concerning asymmetries in the elements' distribution. For Si, S, Ar, Fe and Ca, the GMCA was able to retrieve two maps associated with spectra slightly blueor red-shifted from their theoretical position. The existence of blue-or red-shifted parts in these elements' distribution was previously known, and the Fe maps from Picquenot et al. (2019) were consistent with prior works but endowed with more details (see Willingale et al. 2002;. Thus, they constitute a great basis for an extensive study of the asymmetries in the elements' distribution in Cas A. In this paper, we will use a more recent version of the GMCA, the pGMCA, that was developed to take into account data of a Poissonian nature (Bobin et al. 2020). In the precedent version of the algorithm, the noise was supposed to be Gaussian. Even with that biased assumption, the results were proven to be reliable. However, a proper treatment of the noise is still relevant: It increases the consistency of the spectral morphologies of the retrieved components and makes the algorithm able to disentangle components with a fainter contrast.
The mathematical formalism is similar to that of the GMCA, presented in Picquenot et al. (2019). The fundamental difference is that instead of a linear representation, the pGMCA uses the notion of a Poisson-likelihood of a given sum of components to be the origin of a certain observation. The problem solved by the algorithm is thus essentially the same kind, the main difference being a change in the nature of the norm that needs to be minimized. A more precise description of this new method is available in Bobin et al. (2020).
The use of the pGMCA is also highly similar to that of the GMCA. One notable difference is that the pGMCA is more sensitive to the initial conditions, so it needs a first guess for convergence purposes. The analysis therefore consists of two steps: a first guess obtained with the GMCA and a refinement step using the Poissonian version pGMCA.
The aforementioned workflow was applied to the Cas A Chandra observations by creating data cubes for each energy band shown in large enough to have the leverage to allow the synchrotron continuum to be correctly retrieved and to be narrow enough to avoid contamination by other line emissions. The pGMCA being a fastrunning algorithm, the final energy bands were chosen after tests to find the best candidates for both criteria. For each band, the initial number of components n was 3: the synchrotron emission and the blue-and red-shifted parts of the line emission. We then tested using 4 and 5 components to ensure extra components were not merged into our components of interest. We also tested with 2 components to verify our assumption on the presence of blue-and red-shifted parts was not imposing the apparition of a spurious component. For each emission line, we then chose n as the best candidate to retrieve the most seemingly meaningful components without spurious images. For each analysis, the algorithm was able to retrieve a component that we identify as the synchrotron emission (a power-law spectrum and filamentary spatial distribution, not shown here) and multiple additional thermal components with strong line features. We were able to identify two associated images with shifted spectra from the theoretical emission line energy for all these line features except O, Fe L, and Mg.

Quantification of asymmetries
We use the power-ratio method (PRM) to quantitatively analyze and compare the asymmetries of the images extracted by pGMCA. This method was developed by Buote & Tsai (1995) and previously employed for use on SNRs (Lopez et al. 2009a(Lopez et al. ,b, 2011. It consists of calculating multipole moments in a circular aperture positioned on the centroid of the image, with a radius that encloses the whole SNR. Powers of the multipole expansion P m are then obtained by integrating the mth term over the circle. To normalize the powers with respect to flux, they are divided by P 0 , thus forming the power ratios P m /P 0 . For a more detailed description of the method, see Lopez et al. (2009b).
The P 2 /P 0 and P 3 /P 0 terms convey complementary information about the asymmetries in an image. The first term is the quadrupole power-ratio and quantifies the ellipticity/elongation of an extended source, while the second term is the octupole power-ratio and is a measure of mirror asymmetry. Hence, both are to be compared simultaneously to ascertain the asymmetries in different images.
Here, as we want to compare asymmetries in the blue-and red-shifted part of the elements' distribution, the method is slightly modified. In a first step, we calculate the P 2 /P 0 and P 3 /P 0 ratios of each element's total distribution by using the sum of the blue-and red-shifted maps as an image. Its centroid is then an approximation of the center-of-emission of the considered element. Then, we calculate the power ratios of the blue-and red-shifted images separately using the same center-of-emission. Ultimately, we normalize the power ratios thus obtained by the power ratios of the total element's distribution: where i = 2 or 3 and P i /P 0 (red or blue image) is calculated using the centroid of the total image. That way, we can compare the relative asymmetries of the blue-and red-shifted parts of different elements, without the comparison being biased by the original asymmetries of the whole distribution.

Error bars
As explained in Picquenot et al. (2019), error bars can be obtained by applying this method on every image retrieved by the GMCA applied on a block bootstrap resampling. However, as was shown in that paper, this method introduces a bias in the results of the GMCA. We show in Appendix B that the block bootstrap method modifies the Poissonian nature of the data, thus impacting the results of the algorithm. Since the pGMCA is more dependent than GMCA on the initial conditions, the bias in the outputs is even greater with this newer version of the algorithm (see Fig. B.3). For that reason, we developed a new resampling method we named "constrained bootstrap", presented in Appendix B.4. Thus, we applied pGMCA on a hundred resampled data cubes obtained thanks to the constrained bootstrap and plotted the different spectra we retrieved around the ones obtained on real data. As stated in Appendix B.4, the spread between the resamplings has no physical significance but helps in evaluating the robustness of the algorithm around a given set of original conditions. The blue-shifted part of the Ca line emission, a very weak component, was not retrieved for every resampling. In this case, we created more resamplings in order to obtain a hundred correctly retrieved components. The faintest components are the ones with the largest relative error bars, as can be seen in Figs. 3 and 4, highlighting the difficulty for the algorithm to retrieve them in a consistent way on a hundred slightly different resamplings.
To obtain the error bars for the PRM plot of the asymmetries, we applied the PRM to the hundred images retrieved by the pGMCA on the resamplings. Then, in each direction we plotted error bars representing the interval between the 10th and the 90th percentile and crossing at the median. We also plotted the PRM applied on real data. Although our new constrained bootstrap method ensures the Poissonian nature of the data to be preserved in the resampled data sets, we see that the results of the pGMCA on real data are sometimes not in the 10th-90th percentile zone, thus suggesting there may still be some biases. It happens mostly with the weakest components, showing once more the difficulty for the pGMCA to retrieve them consistently out of different data sets presenting slightly different initial conditions. However, even when the results on real data are not exactly in the 10th-90th percentile zone, the adequation between the results on real and resampled data sets is still good, and the relative positioning for each line is the same, whether we consider the results on the original data or on the resampled data sets.

Images retrieved by pGMCA
By applying the pGMCA algorithm on the energy bands surrounding the eight emission lines shown in Fig. 1, we were able to retrieve maps of their spatial distribution associated with spectra, successfully disentangling them from the synchrotron emission or other unwanted components. The O, Mg, and Fe L lines were only retrieved as single features, each associated with a spectrum, whereas Si, S, Ar, Ca, and Fe-K were retrieved as two different images associated with spectra that we interpret as being the same emission lines slightly red-or blue-shifted. Figure 2 shows the total images for all eight line emissions, obtained by summing the blue-and red-shifted parts when necessary. It also indicates the centroid of each image that is adopted in the PRM. Figure 3 shows the red-and blue-shifted parts of five line emissions, together with their associated spectra, while Fig. 4 presents the images of O, Mg, and Fe L together with their respective spectra. Figure 2 is similar to Fig. 10 from Picquenot et al. (2019), but the images here are more accurate and less contaminated by other components thanks to a proper treatment of the Poisson noise, and the associated spectra not shown in our first paper are presented here in Figs. 3 and 4.

Discussion on the retrieved images
The fact that our algorithm fails to separate a blue-shifted from a red-shifted part in the O, Mg, and Fe L images is not surprising. At 1 keV, we infer that a radial speed of 4000 km s −1 would lead to a ∆E of about 13 eV, which is below the spectral bin size of our data. We see in Fig Their spatial distributions appears similar in Fig. 2, and the division into a red-and a blue-shifted part (as found by the pGMCA) allows us to investigate their three-dimensional morphology. We also notice that the maps of Si and Ar are similar to that of the ArII in infrared (DeLaney et al. 2010).
As the reverse shock has not fully propagated to the interior of Cas A (Gotthelf et al. 2001;), our images may not reflect the full distribution of the ejecta. However, Hwang & Laming (2012) estimates that most of the ejecta mass has already been shocked: we can thus conclude that our images capture the bulk of the ejecta and our element images are likely similar to the true ejecta distributions. In addition, our Fe red-shifted image matches well with the 44 Ti image, produced by radioactive decay instead of reverse-shock heating (see Sect. 4.5).
Hence, we can quantify the asymmetries in the ejecta distribution by using the PRM method described in Sect. 2.3 on our images. Figure 5 presents the quadrupole power-ratios P 2 /P 0 versus the octupole power-ratios P 3 /P 0 of the total images from Fig. 2. Figure 6 shows the quadrupole power-ratios versus the octupole power-ratios of the red-and blue-shifted images presented in Fig. 6 normalized with the quadrupole and octupole power-ratios of the total images (Fig. 2) as defined in Eq. (2.3).

Discussion on the retrieved spectra
As stated before, it is the spectra retrieved together with the aforementioned images that allow us to identify them as "blue-" or "red-shifted" components. Here we will expand on our reasons for supporting these assertions.
The spectra in Fig. 3 are superimposed with the theoretical positions of the main emission lines in the energy range. In the case of Si, the retrieved features are shifted to lower or higher energy with respect to the rest-energy positions of the Si XIII and Si XIV lines. Appendix A shows that this shifting is not primarily due to an ionization effect as the ratio Si XIII/Si XIV is roughly equal in both cases. The same goes for S, where two lines corresponding to S XV and S XVI are shifted together while keeping a similar ratio.
A word on the Ca blue-shifted emission: This component is very weak and in a region where there is a lot of spatial overlap, making it difficult for the algorithm to retrieve. For that reason, the retrieved spectrum has a poorer quality than the others, and it was imperfectly found on some of our constrained bootstrap resamplings. Consequently, we were compelled to run the algorithm on more than a hundred resamplings and to select the accurate ones to obtain a significant envelop around the spectrum obtained on the original data.

Spectral analysis
Using the spectral components retrieved for each data subset shown in Fig. 3, we carried out a spectral fitting assuming a

Red-shift
Blue-shift Si S Ar Ca Fe K Fig. 3. Red-and blue-shifted parts of the Si, S, Ar, Ca, and Fe line emission spatial distribution and their associated spectrum as found by pGMCA. The spectra in red correspond to the application of the algorithm on real data, while the dotted gray spectra correspond to the application on a hundred constrained bootstrap resamplings illustrating statistical uncertainties. The x-axis is in keV and the y-axis in counts. The dotted vertical lines represent the energy of the brightest emission lines for a non-equilibrium ionization plasma at a temperature of 1.5 keV and ionization timescale of log(τ) = 11.3 cm −3 s produced using the AtomDB (Foster et al. 2012). These parameters are the mean value of the distribution shown in Fig. 2 of Hwang & Laming (2012).
residual continuum plus line emission in XSPEC (power-law + Gauss model). In this analysis, the errors for each spectral data point are derived from the constraint bootstrap method presented in Appendix B. This constrained bootstrap eliminates a bias introduced by classical bootstrap methods and that is critical to pGMCA, but underestimates the true statistical error. Therefore, no statistical errors on the line centroids are listed in Table 1 as, in addition, systematic errors associated with ACIS energy calibration are likely to be the dominant source of uncertainty.
The resulting line centroid and equivalent velocity shifts are shown in Table 1. To transform the shift in energy into a velocity shift, a rest energy is needed. The ACIS CCD spectral resolution does not resolve the line complex and cannot easily disentangle velocity and ionization effects. However, given the range of ionization state observed in Cas A (with ionization ages of τ ∼ 10 11 −10 12 cm −3 s, see and their associated spectra as found by pGMCA. The spectra in red correspond to the application of the algorithm on real data, while the dotted gray spectra correspond to the application on a hundred constrained bootstrap resamplings. The x-axis is in keV and the y-axis in counts. . Quadrupole power-ratios P 2 /P 0 versus the octupole power-ratios P 3 /P 0 of the total images of the different line emissions shown in Fig. 2. The dots represent the values measured for the pGMCA images obtained from the real data, and the crosses the 10th and 90th percentiles obtained with pGMCA on a hundred constrained bootstrap resamplings, with the center of the cross being the median. rest energy was chosen as the brightest line for a non-equilibrium ionization plasma with a temperature of 1.5 keV temperature and log(τ) = 11.3 cm −3 s, the mean values from Fig. 2 of Hwang & Laming (2012).
For the specific case of the Si XIII line, a very large asymmetry in the red/blue-shifted velocities is observed. This could be due to possible energy calibration issues near the Si line as shown by  in a comparison of ACIS and HETG line centroid, resulting in a systematic blue-shift effect in ACIS data. The Si XIII * line in Table 1 Fig. 6. Quadrupole power-ratios P 2 /P 0 versus the octupole power-ratios P 3 /P 0 of the red-and blue-shifted images of the different line emissions shown in Fig. 3, normalized with the quadrupole and octupole powerratios of the total images. The dots and error bars are obtained in the same way as in Fig. 5. line energy to illustrate systematic uncertainties associated with calibration issues. For the Fe-K complex of lines, we rely on the analysis of DeLaney et al. (2010) who derived an average rest line energy of 6.6605 keV (1.8615 Å) by fitting a spherical expansion model to their three-dimensional ejecta model. We can note that with this spectral analysis, what we measure here is the radial velocity that is flux weighted over the entire image of the associated component. Therefore, we are not probing the velocity at small angular scale but the bulk velocity of the entire component.
With the caveats listed above, we notice an asymmetry in the velocities where ejecta seem to have a higher velocity toward us (blue-shifted) than away from us, even in the case of Si XIII after calibration corrections. A comparison of those results with previous studies and possibles biases are discussed in Sect. 4.3. The large uncertainties associated with the energy calibration and the choice of rest energy has little impact on the delta between the red-and blue shifted centroids and hence on the ∆V. We note that all elements show a consistent ∆V of ∼6000 km s −1 .

Physical interpretation
4.1. Quantification of ejecta asymmetries Figure 5 shows that the distribution of heavier elements is generally more elliptical and more mirror asymmetric than that of A82, page 6 of 13 A. Picquenot et al.: Morphological asymmetries in the ejecta of Cassiopeia A in X-rays  (2012), Mg emission does not follow the exact same trend as the other elements: it has roughly an order of magnitude lower elliptical asymmetry (P 2 /P 0 ) than the other elements. In contrast to Holland-Ashford et al. (2020) and Hwang & Laming (2012), our Mg image (as shown in Fig. 4) presents a morphology highly different from that of the Fe L; we believe that the pGMCA was able to retrieve the Mg spatial distribution with little continuum or Fe contamination. Figure 6 presents the relative ellipticity/elongation and mirror asymmetries of the blue-and red-shifted ejecta emission compared to the total ejecta images (Fig. 2). A value of "1" indicates that the velocity-shifted ejecta has equivalent levels of asymmetry as the full bandpass emission. In the cases where we can clearly disentangle the red-and blue-shifted emission (i.e.Si, S, Ar, Ca, and Fe-K, described in previous paragraphs), we see that the red-shifted ejecta emission is less asymmetric than the blue-shifted emission. This holds true both for elliptical asymmetry P 2 /P 0 and mirror asymmetry P 3 /P 0 . Thus, we could physically describe the red-shifted ejecta distribution as a broad, relatively symmetric plume, whereas the blue-shifted ejecta is concentrated into dense knots. This interpretation matches with the observation that most of the X-ray emission is from the red-shifted ejecta, as we can also see in the flux ratios shown in Table 2 and in the images of Fig. 3, suggesting that there was more mass ejected away from the observer, NS, and blue-shifted ejecta knot. We note that there is not a direct correlation between ejecta mass and X-ray emission due to the position of the reverse shock, the plasma temperature and ionization timescale, but the indication that most of the X-ray emission is red-shifted is consistent with our knowledge of the 44 Ti distribution (see Sect. 4.5 for a more detailed discussion).
Furthermore, in all cases, the red-shifted ejecta emission is more circularly symmetric than the total images, and the blueshifted ejecta is more elliptical and elongated than the total images. Moreover, the red-shifted ejecta is more mirror symmetric than the blue-shifted ejecta, though both the red-shifted and blue-shifted Si are more mirror asymmetric than the total image. The latter result may suggest that the red-shifted and blue-shifted Si images' asymmetries sum together such that the total Si image appears more mirror symmetric than the actual distribution of the Si.  Fig. 13 of Grefenstette et al. (2017). Only the direction is relevant as the norm of this specific vector is arbitrary. Figure 7 shows the centroids of the blue-and red-shifted parts of each emission line relative to the center-of-explosion of Cas A, revealing the bulk three-dimensional distribution of each component. We note that this figure was only made using the centroids of the red-and blue-shifted images retrieved by pGMCA, without using the PRM method. We can see the red-shifted ejecta is mainly moving in a similar direction (toward the northwest), while all the blue-shifted ejecta is moving toward the east. As discussed in Sect. 4.5, this result is consistent with previous works on Cas A investigating the 44 Ti distribution with NuSTAR data (Grefenstette et al. 2017).

Three-dimensional distribution of heavy elements
The blue-shifted ejecta is clearly moving in a different direction than the red-shifted ejecta, but not directly opposite to it. The angles between the blue-and red-shifted components are all between 90 • and 140 • . This finding provides evidence against a jet and counter-jet explosion mechanism being responsible for the explosion and resulting in the expansion of ejecta in Cas A (e.g., Fesen 2001;Hines et al. 2004;Schure et al. 2008). We can also note a trend where heavier elements exhibit increasingly larger opening angles than lighter elements, from Si showing a 90 • angle to Ca and Fe that show opening angles of about 130 • −140 • .

Velocities of red-and blue-shifted structures
By fitting the line centroids, we obtained the velocities discussed in Sect. 3.4. As stated before, the effects of ionization on possible "imposter velocities" are discussed in Appendix A. Our derived velocities showed higher values for the blue-shifted component than for the red-shifted one for all elements. Those results are in disagreement with spectroscopic studies and in agreement with some others. On the one hand, the X-ray studies of individual regions Willingale et al. (2002) (Fig. 8 to the spectrum extracted in each small-scale region. In regions where both red-and blue-shifted ejecta co-exist (see Fig. 3), the Gaussian fit will provide a flux weighted average velocity value of the two components as they are not resolved with ACIS. As the red-shifted component is brighter in average, a systematic bias that would reduce the blue velocities could exist. This could be the case in the southeastern region where most of the blue-shifted emission is observed and where a significant level of red-shifted emission is also seen. Besides this, calibration issues may also play an important role. Although the GMCA method was successful in retrieving the centroid energy of nearby emission lines using a simple toy model (Picquenot et al. 2019, Fig. 7), we do not rule out that the higher velocity of the blue-shifted component is an artifact of the method. Further tests of the method with the help of synthetic X-ray observations using numerical simulations could shed light on this issue.

Neutron star kick direction
The NS in Cas A is located southeast of the explosion site, moving at a velocity of ∼340 km s −1 southeast in the plane of the sky (Thorstensen et al. 2001;Fesen et al. 2006). In Hwang & Laming (2012) it was stated that, contrary to expectations, the Fe structure was not observed to recoil in the opposite direction to the NS. Here, thanks to our ability to disentangle red-and blue-shifted structures, we find that the red-shifted ejecta is moving nearly opposite the NS : the angles between the red-shifted structures and the NS tangential motion range between 154 • (Fe) and 180 • (Si). Table 2 also shows that the bulk emission is from red-shifted ejecta (consistent with Milisavljevic & Fesen 2013). This correlation is consistent with theoretical predictions that NSs are kicked opposite to the direction of bulk ejecta motion, in adequation with conservation of momentum with the ejecta (Wongwathanarat et al. 2013;Müller 2016;Bruenn et al. 2016;Janka 2017). Specifically, observations have provided evidence for the "gravitational tugboat mechanism" of generating NS kicks asymmetries proposed by Wongwathanarat et al. (2013); Janka (2017), where the NS is gravitationally accelerated by the slower moving ejecta clumps, opposite to the bulk ejecta motion. It is impossible to calculate the NS line-of-sight motion by examining the NS alone as its spectra contains no lines to be Doppler-shifted. However, limits on its three-dimensional motion can be placed by assuming it moves opposite the bulk of ejecta and examining the bulk three-dimensional motion of ejecta. Grefenstette et al. (2017) studied Ti emission in Cas A and found that the bulk Ti emission was tilted 58 • into the plane of the sky away from the observer, implying that the NS is moving 58 • out of the plane of the sky toward the observer. This finding is supported by three-dimensional simulations of a Type IIb progenitor by Wongwathanarat et al. (2017) and Jerkstrand et al. (2020), which suggested that the NS is moving out of the plane of the sky with an angle of ∼30 • .
The results of this paper support the hypothesis that, if the NS is moving away from the bulk of ejecta motion, the NS is moving toward us. Furthermore, we could tentatively conclude that the NS was accelerated toward the more slowly moving blueshifted ejecta, which would further support the gravitational tugboat mechanism. The strong levels of asymmetry exhibited by the blue-shifted emission combined with the lower flux would imply that the blue-shifted ejecta is split into relatively small ejecta clumps, one of which would possibly be the source of the neutron star's gravitational acceleration. However, the velocities determined in Table 1 contradict this hypothesis as the blue-shifted clumps seem to move faster.

Comparison with 44 Ti
44 Ti is a product of Si burning and is thought to be synthesized in close proximity with iron. The 44 Ti spatial distribution has been studied via its radioactive decay with the NuSTAR telescope and revealed that most of the material is red-shifted and does not seem to follow the Fe-K X-ray emission (Grefenstette et al. 2014(Grefenstette et al. , 2017. In our study, we have found that 70% of the Fe-K X-ray emission (see Table 2) is red-shifted and that the mean direction of the Fe-K red-shifted emission shown in Fig. 7 is compatible with that of the 44 Ti as determined in Fig. 13 of Grefenstette et al. (2017). Yet, we can see the mean 44 Ti direction is not perfectly aligned with the mean red-shifted Fe-K direction. This may be caused by the fact that the Fe-K emission is tracing only the reverse shock-heated material and may not reflect the true distribution of Fe, whereas 44 Ti emission is from radioactive decay and thus reflects the true distribution of Ti.
In Fig. 8, we overlay the ten regions where Grefenstette et al. (2017) detected 44 Ti with our red-shifted component image. The regions 19 and 20 (which dominate our Fe-K red-shifted component image) have respective 44 Ti velocities of 2300 ± 1400 and 3200 ± 500 km s −1 , values that are compatible with our measured value of ∼2800 km s −1 shown in Table 1.
Concerning our Fe-K blue-shifted component map, its X-ray emission is fainter and located mostly in the southeast of the source (see Fig. 3). This southeastern X-ray emission is spatially coincident with region 46 in the Fig. 2 NuSTAR map of Grefenstette et al. (2017), not plotted in our Fig. 8 as the 44 Ti emission was found to be below the detection threshold.
We note that blue-shifted 44 Ti emission is harder to detect for NuSTAR than a red-shifted one as it is intrinsically fainter. In addition, any blue-shifted emission of the 78.32 keV 44 Ti line places it outside the NuSTAR bandpass, precluding detection of one of the two radioactive decay lines in this case.

Conclusions
By using a new methodology and applying it to Cas A Chandra X-ray data, we were able to revisit the mapping of the heavy elements and separate them into a red-and a blue-shifted parts, allowing us to investigate the three-dimensional morphology of the SNR. These new maps and the associated spectra could then be used to quantify the asymmetries of each component, their mean direction and their velocity. The main findings of the paper are consistent with the general results found in the previous studies cited in Part 4, and are summarized below: -Morphological asymmetries: an extensive study of the asymmetries shows the distribution of heavier elements is generally more elliptical and mirror asymmetric in Cas A, which is consistent with simulation predictions. For the elements we were able to separate into a red-and a blue-shifted parts (Si, S, Ar, Ca, Fe), it appears that the red-shifted ejecta is less asymmetric than the blue-shifted one. The red-shifted ejecta can then be described as a broad, relatively symmetric plume, while the blue-shifted ejecta can be seen as concentrated into dense knots. Most of the emission from each element is red-shifted, implying there was more mass ejected away from the observer, which agrees with past studies. -Three-dimensional distribution: the mean directions of the red-and blue-shifted parts of each element are clearly not diametrically opposed, disfavouring the idea of a jet/counterjet explosion mechanism. -NS velocity: we find that the NS is moving most opposite to the direction of the red-shifted ejecta that forms the bulk of the ejecta emission. This supports the idea of a "Gravitational Tugboat Mechanism" of generating NS kicks through a process consistent with conservation of momentum between NS and ejecta. This result implies that the NS is moving toward us, which is consistent with the findings of past studies. However, we find the blue-shifted clumps to be faster than the red-shifted ones, which is not consistent with the gravitational tug-boat mechanism's prediction that the NS is moving opposite to the faster ejecta. -Comparison with 44 Ti: our finding that the bulk of ejecta is red-shifted and moving NW is consistent with the 44 Ti distribution from NuSTAR observations. Its direction is similar to that of the red-shifted Fe-K emission, but a slight difference could be explained by the fact that the Fe-K only traces the reverse shock-heated ejecta and not the full distribution of the Fe ejecta. The component separation method presented here enabled a three-dimensional view of the Cas A ejecta despite the low energy resolution of the Chandra CCDs by separating entangled components all at once, without the need of a detailed spectral analysis on hundreds of regions. In the future, X-ray microcalorimeters will enable kinematic measurements of X-ray emitting ejecta in many more SNRs. In its short operations, the Hitomi mission demonstrated these powerful capabilities. In particular, in a brief 3.7-ks observation, it revealed that the SNR N132D had highly red-shifted Fe emission with a velocity of ∼800 km s −1 without any blue-shifted component, suggesting the Fe-rich ejecta was ejected asymmetrically (Hitomi Collaboration 2018). The upcoming replacement X-ray Imaging and Spectroscopy Mission XRISM will offer 5-7 eV energy resolution with 30 pixels over a 3 field of view (Tashiro et al. 2018). In the longer term, Athena and Lynx will combine this superb spectral resolution with high angular resolution, fostering a detailed, three-dimensional view of SNRs that will revolutionize our understanding of explosions (Lopez et al. 2019;Williams et al. 2019). While the new instruments will provide a giant leap forward in terms of data quality, development of new analysis methods are needed in order to maximize the scientific return of next generation telescopes.

B.1. Introduction
The BSS method we used in this paper, the pGMCA, is one of the numerous advanced data analysis methods that have recently been introduced for a use in astrophysics, among which we can also find other BSS methods, classification, PSF deconvolution, denoising or dimensionality reduction. We can formalize the application of these data analysis methods by writing Θ = A(X), where X is the original data, A is the nonlinear analysis operator used to process the signal and Θ is the estimator for which we want to find errors (in this paper for example, X is the original X-ray data from Cas A, A is the pGMCA algorithm and Θ represents the retrieved spectra and images). Most of these methods being nonlinear, there is no easy way to retrieve error bars or a confidence interval associated with the estimator Θ. Estimating errors accurately in a nonlinear problem is still an open question that goes far beyond the scope of astrophysical applications as there is no general method to get error bars from a nonlinear data-driven method such as the pGMCA. This is a hot topic whose study would be essential for an appropriate use of complex data analysis methods in retrieving physical parameters, and for allowing the user to estimate the accuracy of the results.

B.2. Existing methods to retrieve error bars on Poissonian data sets
Our aim, when searching for error bars associated with a certain estimator Θ on an analyzed data set, is to obtain the variance of Θ = A(X) , where the original data X is composed of N elements. When working on a simulation, an obvious way to proceed in order to estimate the variance of Θ is to apply the considered data analysis method A on a certain number of Monte-Carlo (MC) realizations X i and look at the standard deviation of the results Θ i = A(X i ). The variance of the Θ i provides a good estimation of the errors. Yet, this cannot be done with real data as only one observation is available: the observed one. Thus, a resampling method such as the jackknife, the bootstrap (see Efron 1979) or its derivatives, able to simulate several realizations out of a single one, is necessary. Ideally, the aim is to obtain through this resampling method a number of "fake" MC realizations centered on the original data: new data sets variating spatially and reproducing the spread of MC drawings with a mean equal or close to the mean of the original data. The mechanisms at stake in jacknife or bootstrap resamplings are similar. Jacknife and bootstrap resampling methods produce n resampled setsX i by rearranging the elements of X, and allow us to consider the variance of Θ i = A(X i ) for i in 1, n as an approximation for the variance of Θ. As jacknife and bootstrap methods are close to each other, and the bootstrap and some of its derivatives are more adapted to handle correlated data sets, we will in this appendix focus on a particular method, representative of other resampling methods and theoretically suited for astrophysical applications: the block bootstrap, which is a simple bootstrap applied on randomly formed groups of events rather than on the individual events.
In the case of a Poisson process, the discrete nature of the elements composing the data set can easily be resampled with a block bootstrap method. The N discrete elements composing a Poissonian data set X will be called "events". In X-rays for example, the events are the photons detected by the spectro-imaging instrument. The bootstrap consists in a random sampling with replacement from the current set events X. The resampling obtained through bootstrapping is a setX boot of N events taken randomly with replacement amid the initial ones (see Fig. B.1). This method can be repeated in order to simulate as many real-izationsX boot i as needed to estimate standard errors or confidence intervals. In order to save calculation time, we can choose to resample blocks of data of a fixed size instead of single events: This method is named block bootstrap. The block bootstrap is also supposed to conserve correlations more accurately, making it more appropriate for a use on astrophysical signals. The data can be of any dimension but for clarity, we will only show in this Appendix bi-dimensional data sets, that is images.

B.3. Biases in classical bootstrap applied on Poissonian data sets
The properties of the data resampled strongly depends on the nature of the original data. Biases may appear in the resampled data sets, proving a block bootstrap can fail to reproduce consistent data that could be successfully used to evaluate the accuracy of certain estimators. In particular, Poissonian data sets, including our X-rays data of Cas A, are not consistently resampled by current resampling methods such as the block bootstrap. A Poissonian data set X can be defined as a Poisson realization of an underlying theoretical . On the right, the black histogram correspond to the original data X = P(X * ). The red histograms are those of the data setsX boot i obtained through resampling of the original data and the blue ones are the histograms of a Poisson realization of the original data P(X) = P(P(X * )). It appears that the resampled data sets have histograms highly similar to that of the original data with additional Poisson noise. model X * , which can be written: where P(.) is an operator giving as an output a Poisson realization of a set.
A look on the histogram of a data set resampled from a Poissonian signal shows the block bootstrap fails to reproduce accurately the characteristics of the original data. Figure B.2, top, compares the histogram of the real data X, a simple image of a square with Poisson noise, with the histograms of the resampled data setsX boot i , and highlights the fact that the latter are more similar to the histogram of a Poisson realization of the original data P(X) = P(P(X * )) than to the actual histogram of the original data X = P(X * ), where X * is the underlying model of a square before adding Poisson noise. This is consistent with the fact that the block bootstrap is a random sampling with replacement, which introduces uncertainties of the same nature as a Poisson drawing. Figure B.2, bottom, shows the comparison between the histogram of the toy model Cas A image and the histograms of the data sets resampled with a block bootstrap. We can see the resampling is, in this case too, adding Poissonian noise and gives histograms resembling P(P(X * )) rather than P(X * ). The same goes with our real data cube of Cas A: Fig. B.3 shows an obvious instance of this bias being transferred to the results of the pGMCA, thus proving the block bootstrap cannot be used as such to retrieve error bars for this algorithm.

B.4. A new constrained bootstrap method
Bootstrap resamplings consisting in random drawings with replacement, it is natural that they fail to reproduce some characteristics of the data, among which the histogram that gets closer There is an obvious bias in the results, the resampled data spectra being consistently underestimated.

Original image
Defining a new histogram Imposing the new histogram on the original image to the histogram of a Poisson realization of the original data than to the histogram of the actual data. The block bootstrap method is therefore unable to simulate a MC centered on the original data: the alteration of the histogram strongly impacts the nature of the data, hence the differences in the morphologies observed by looking at the wavelet coefficients. It is then necessary to find a new method in which we could force the histogram of the resampled data sets to be similar to that of the original data.
A natural way to do so would be to impose the histogram we want the resampled data to have before actually resampling the data. To allow this constraint to be made on the pixel distribution, we can no longer consider our events to be the individual elements of X or a block assembling a random sample of them. We should directly work on the pixels and their values, the pixels here being the basic bricks constituting our data. Just as the block bootstrap, our new method can work with data of any dimension. In the case of images, the "basic bricks" correspond to actual pixels values. In the case of X-rays data cubes, they are tiny cubes of the size of a pixel along the spatial dimensions, and the size of an energy bin along the spectral dimension. The same goes for any dimension of our original data. The method can also be adapted for uni-dimensional data sets. The key of our new method is then to work on the histogram of the data presenting the pixels' values rather than on the data itself, event by event. We can either change the value of a pixel or exchange the value of a pixel for that of another one. The first operation simultaneously adds and subtracts 1 in the corresponding columns of the global histogram while the second operation does not produce any change in it. A good mixture of these two operations would then allow us to obtain the histogram we want to impose in our resampled data sets, and following a Poisson probability law to select the pixels to exchange would introduce some spatial variations, in order to reproduce what a MC would do.
Our new constrained bootstrap method is thus composed of two steps, that are described below and illustrated in Fig. B.4: Firstly, obtaining the probability density function of the random variable underlying the observed data histogram using the Kernel Density Estimation (KDE), and randomly generating n histograms from this density function with a spread around the data mimicking that of a MC, with a constraint enforcing a Poissonian distribution of the total sums of pixel values of the n histograms.
Secondly, producing resampled data sets associated with the new histograms by changing the values of wisely chosen pixels in the original image.
During these steps, the pixels equal to zero remain equal to zero, and the nonzero pixels keep a strictly positive value. This constraint enforces the number of nonzero pixels to be constant and avoids the creation of random emergence of nonzero pixels in the empty area of the original data. While this is not completely realistic we prefer constraining the resampled data sets in this way than getting spurious features. We could explore ways to release this constraint in the future. Figure B.5 highlights the similarities between the original histogram and those obtained through MC realizations and our new constrained bootstrap resamplings, while Fig. B.6 and the spectra in Figs. 3 and 4 show that even after being processed by the sensitive pGMCA algorithm, this resampling method shows little to no biases. Hence, our new constrained bootstrap method brings a first and successful attempt at solving the problem of biases in bootstrapping Poissonian data sets.
The comparison of our resampled data sets to a group of MC realizations of the same simulation of Cas A appears to be promising for the variance induced by our method. However, when applying a complex estimator such as the pGMCA on both the MC realizations and our resampled data sets, it appears that the variances obtained through our method fail to accurately reproduce those of the MC realizations. For that reason, the error bars retrieved by our constrained bootstrap method do not have a physical signification. Nevertheless, they constitute an interesting way to assess the robustness of our method around a certain line emission. The different resamplings explore initial conditions slightly different from the original data, thus evaluating the dependence of our results on the initial conditions. Figures 3 and 4 indeed show that for some line emissions, the dispersion between the results on different resampled data sets is far greater than for others.
This new constrained bootstrap method is a first and promising attempt at retrieving error bars for nonlinear estimators on Poissonian data sets, a problem that is often not trivial. In nonlinear processes, errors frequently cannot be propagated correctly, so the calculation of sensitive parameters and the estimation of errors after an extensive use of an advanced data analysis could benefit from this method. We will work in the future on a way to constraint the variance of the results to be more closely related to that of a set of MC realizations in order to ensure the physical signification of the obtained error bars.