© ESO 2020

1 Introduction

Core-collapse (CC) supernovae (SNe), the final fate of massive stars, are known to play a major role in the dynamical and chemical evolution of galaxies (by injecting mass and energy) and in driving the chemical enrichment of the diffuse gas. Nevertheless, many aspects of the processes governing the SN engine and the final stages in the evolution of their progenitor stars are not fully understood. Severe limitations in these studies are the rarity of SN events in our galaxy (on average about one every 50 yr; Diehl et al. 2006), the large distances of extragalactic SNe which make them unresolved point-like sources, and their unpredictability and transient nature. These issues make it extremely challenging to extract key information on the progenitor stars and on the explosion processes associated with SNe from observations.

Otherwise, nearby supernova remnants (SNRs), the objects left over from SN explosions, are extended sources for which the structure and chemical composition of the stellar debris ejected in the SN outburst (i.e., the ejecta) can be studied in detail. More specifically, young and nearby SNRs encode valuable information in their expanding debris and evolution that can reveal details about the inherent asymmetries of the explosion itself and the products of nucleosynthesis occurring during the explosion. Moreover, SNRs probe the circumstellar medium (CSM) surrounding the SNe, which can be shaped by the progenitor stars through their powerful stellar winds. However, with increasing age, all these pieces of information, which are useful for establishing a SN–SNR connection, are diluted by the increasing interaction of the remnant with the surrounding environment. Therefore, young and nearby remnants are most likely the best candidates to investigate the intimate link that exists between the morphological properties of a SNR and the complex phases in the SN explosion, although in the Milky Way less than a dozen SNRs are <1000 yr old (Ferrand & Safi-Harb 2012; Green 2014).

Hydrodynamic and magnetohydrodynamic (HD–MHD) models can be powerful tools to bridge the gap between CC–SNe and their remnants. However, up until recent years, the key strategy was to describe either the SN evolution or the expansion of the remnant due to the very different time and space scales involved in these two phases of evolution and the inherent three-dimensional (3D) nature of the phenomenon. As a result, SN models described the early SN evolution, leaving out an accurate description of their subsequent interaction with the ambient environment; SNR models described only the expansion of the remnant and its interaction with the ambient medium, leaving out a self-consistent and accurate description of the distribution of mass and energy of the ejecta soon after the SN explosion. Only recently, seminal studies (Patnaude et al. 2015; Utrobin et al. 2015; Orlando et al. 2015, 2016, 2019, 2020; Wongwathanarat et al. 2017; Ferrand et al. 2019) have described the whole evolution from the SN to the SNR and have shown that HD–MHD models can be very effective in studying the SN–SNR connection. Thus, these studies bring together two distinct physical scenarios (SN and SNR) on the logical and temporal level, creating continuity between the two.

To date, the mixing of chemically homogeneous layers of ejecta after the SN event has been investigated in detail only in the aftermath of SN explosions (e.g., Kifonidis et al. 2006; Joggerst et al. 2009, 2010), not considering the transition from the phase of SNto that of SNR. Very recently, Ono et al. (2020) studied the matter mixing in aspherical CC-SNe through accurate 3D simulations, aiming to link the ejecta structure soon after the shock breakout at the stellar surface (i.e.,when the shock produced by the explosion leaves the stellar surface) with the nature of the progenitor star.

Here, we aim to analyze how matter mixing occurs at later times during the remnant expansion and interaction with the CSM to investigate how the various chemically homogeneous layers at the time of the explosion map into the resulting abundance pattern of the SNR, and how large-scale anisotropies typically observed in CC–SNe (possibly originated by instabilities developed during the CC; e.g., Nagataki et al. 1997, 1998; Nagataki 2000; Kifonidis et al. 2006; Takiwaki et al. 2009; Gawryszczak et al. 2010; Wongwathanarat et al. 2017; Bear & Soker 2018) affect the evolution of these layers. To this end, we developed a 3D MHD model describing the evolution of the remnant from the onset of the SN to the development of its remnant at the age of 5000 yr, focusing on the case of a progenitor red supergiant (RSG). We analyzed the effects of large-scale anisotropies that may develop soon after the SN on the final remnant morphology. These anisotropies may be responsible for the knotty and clumpy ejecta structure observed at different wavelengths in many CC SNRs (e.g., the Vela SNR, Aschenbach et al. 1995; Miceli et al. 2008; García et al. 2017; G292.0+1.8, Park et al. 2004; Puppis A, Katsuda et al. 2008; and Cassiopeia A, Fesen et al. 2006; Milisavljevic & Fesen 2013; N132D, Law et al. 2020) which is commonly interpreted as a product of HD instabilities developed during the complex phases of the SN explosion (Gawryszczak et al. 2010). However, the role of these anisotropies in the evolution of the chemical stratification of the remnant is still unclear.

The paper is organized as follows: in Sect. 2 we describe the MHD model, in Sect. 3 we discuss the results of the simulations, and in Sect. 4 we draw our conclusions.

2 The numerical setup

The model describes the evolution of a SNR from soon after the CC of the progenitor star to the full-fledged remnant at the age of 5000 yr, following the expansion of the remnant through its pre-SN environment. As in Orlando et al. (2016), the initial distribution of ejecta is derived from a 1D SN model which describes the SN evolution from the CC to the shock breakout at the stellar surface, covering about 24 h of evolution. The output of the SN simulation provided the initial radial profile of the ejecta distribution, including their chemical composition. We then mapped this output into the 3D domain and started the 3D SNR simulations.

2.1 The initial conditions

As initial conditions for our 3D SNR simulations we adopted one of the SN simulations presented in Ono et al. (2020). These authors performed both 1D and 3D HD simulations of SN explosions, with the aim of investigating the matter mixing during the expansion of the blast wave through the stellar interior. Their 1D simulations assume spherical symmetry, whereas the 3D simulations consider aspherical explosions triggered by an initial asymmetric injection of kinetic and thermal energy around the composition interface of the iron core and silicon-rich layer. The simulations were initialized by considering different pre-SN stellar models, including progenitor blue and red supergiants, and were performed with the adaptive-mesh refinement HD–MHD code FLASH (Fryxell et al. 2000).

The SN model takes into account the effects of gravity (both self-gravity and gravitational effects of the central proto-neutron star), the fallback of material onto the central compact object, the explosive nucleosynthesis through a nuclear reaction network including 19 species (n, p, ¹H, ³ He, ⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne, ²⁴Mg, ²⁸Si, ³²S, ³⁶Ar, ⁴⁰Ca, ⁴⁴Ti, ⁴⁸Cr, ⁵²Fe, ⁵⁴Fe and ⁵⁶Ni), the feedback of nuclear energy generation, and the energy deposition due to radioactive decays of isotopes synthesized in the explosion (see Ono et al. 2020 for more details). For our purposes, we considered the case of a 1D SN simulation from a non-rotating, solar-metallicity progenitor RSG with a main-sequence mass of 19.8 M_⊙ which reduces to 15.9 M_⊙ at collapse (Sukhbold et al. 2016), and with injected an energy of 2.5 × 10⁵¹ erg (see Sect. 3.1 in Ono et al. 2020).

The aim of our SNR simulations is not a detailed description of the effects of matter mixing before the shock breakout (indeed carefully modeled and analyzed through the 3D SN simulations by Ono et al. 2020). The aim of this work is to highlight the contributionto matter mixing during the expansion of the blast wave through the CSM and to investigate the stability of the 1D “onion-shell”structure of the progenitor star during the evolution of the SNR. To this end, we explored the case of a spherically symmetric SN explosion and the case of a SN characterized by large-scale anisotropy. In this latter case, once the output of the 1D SN simulation is mapped to the 3D domain, we parametrize the initial large-scale anisotropy of ejecta that may have been developed soon after the SN explosion. The aim is to investigate how the matter mixing and the final remnant morphology depend on the physical (density and velocity) and geometric (position and size) properties of the initial anisotropy. The physical properties of the initial condition, that is, the output from the CC–SN simulation, consist of an initial remnant at t ~ 6 × 10⁴ s (≈ 17 h) after the CC with a radius of ~ 1.85 × 10¹⁵ cm, an ejecta mass of 14.4 M_⊙, and a kinetic energy of 1.35 × 10⁵¹ erg. At this epoch, hours after the breakout of the shock at the stellar surface, the ejecta expand almost homologously through the tenuous wind of the progenitor star (e.g., Gawryszczak et al. 2010; Wongwathanarat et al. 2017; Ono et al. 2020). The initial radial density profile of the ejecta presents a less dense inner region surrounded by a denser shell located in the region of radius 6.2 × 10¹⁴ cm, corresponding to about 30% of the radius of the ejecta (see Figs. 1 and 2; see also Ono et al. 2020 for a description of the evolution of the SN and the formation of the large-scale structure in the density radial profile). This density profile arises from the SN simulation and depends on the progenitor star considered. Here we considered a generic case of a RSG; a different progenitor, such as for instance a blue supergiant, could produce a much smoother density profile with structures that may significantly differ from the dense shell discussed above (see e.g., Ono et al. 2020). The velocity profile (see Fig. 1) increases linearly with the radius of the ejecta from zero at the origin to a maximum of 5.3 × 10⁸ cm s⁻¹ at the external shell of the initial remnant.

Fig. 1 Initial conditions for the 3D simulations. Top panel: radial profiles of initial ejecta density (black line) and pressure (red line). Bottom panel: radial profiles of initial ejecta velocity (black line) and Mach number (red line).

2.2 Modeling the evolution of the SNR

The numerical setup for the SNR simulations is analogous to that described in Orlando et al. (2019, 2020). As the remnant expansionisdescribed during the free expansion and Sedov expansion phases, the radiative losses (see Orlando et al. 2005, also treated in Miceli et al. 2006; Lee et al. 2015) are neglected. We adopted a 3D Cartesian coordinate system (x, y, z) and we carried out the simulations by numerically solving the full time-dependent ideal MHD equations in conservative form: $$\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho v)=0,\\ \frac{\partial \rho v}{\partial t}+\nabla \cdot {\left[\rho vv-BB+I{P}_{\text{*}}\right]}^T=0,\\ \frac{\partial \rho E}{\partial t}+\nabla \cdot [v(\rho E+P_{\text{*}})-v\cdot BB]=0,\\ \frac{\partial B}{\partial t}+\nabla \cdot (vB-Bv)=0,\end{array}$$

where $$P_{\text{*}}=P+\frac{B^2}{2},E=\varepsilon +\frac{1}{2}{v}^2+\frac{1}{2}\frac{B^2}{\rho }$$(5)

are the total pressure and the total gas energy (internal energy, ɛ, kinetic energy, and magnetic energy), respectively, B is the magnetic field, ρ is the density, v is the velocity, I is the identity matrix, and P is the pressure and t is the time. The set of equations is closed by the equation of state (EOS) for ideal gases, P = (γ − 1)ρɛ, where γ = 5∕3 is the adiabatic index.

The 3D simulations of the expanding SNR were performed using PLUTO (see Mignone et al. 2007), a numerical code developed for the solution of hyper-sonic flows that has been widely and successfully used in several studies. The MHD equations are solved using the MHD module available in PLUTO, which was configured to compute inter-cell fluxes with the Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver, while third order in time is achieved using a Runge-Kutta scheme. To preserve the condition of ∇ ⋅ B = 0, a mixed hyperbolic and parabolic divergence cleaning technique is used (Mignone et al. 2007).

The distribution of elements is traced through a set of advection equations: $$\frac{\partial \rho X_{\text{i}}}{\partial t}+\nabla \cdot (\rho X_{\text{i}}v)=0,$$(6)

which are solved for each of the elements considered in addition to the MHD equations, where X_i is the mass fraction of the element of index i, and i runs over the 19 species listed in Sect. 2.1. Unlike the 1D SN simulation used as initial condition, the 3D SNR simulations do not include the radioactive decays of the explosive nucleosynthesis products synthesized in the SN explosion, such as ⁵⁶Ni or ⁴⁴Ti. The main effect of radioactive decay would be a slight inflation of the instability-driven structures developed in regions near the decaying elements (Gabler et al. 2017). As the decay time of these species is much lower than the evolution time analyzed here, we assumed at the beginning of our simulations that all the ⁵⁶Ni and ⁴⁴Ti are already decayed in their decaying products (⁵⁶Fe and ⁴⁴Ca respectively; see Nadyozhin 1994; Ahmad et al. 2006).

Our computational domain describes only one octant of the whole remnant to reduce the computational cost. At the initial condition, the explosion is located at the origin (0,0,0) of the coordinate system. To follow the wide range of space and time scales covered during the expansion of the remnant from the onset of the explosion to the first 5000 yr, we adopted the same approach as Orlando et al. (2019). Initially, the domain extends between 0 and 2.19 × 10¹⁵ cm in all directions and is covered by a uniform grid with 256³ zones. The domain is then gradually extended as the forward shock propagates outward and the physical values are remapped in the new domains. More specifically, each time the forward shock gets close to 90% of the domain size, the domain is extended by a factor of 1.2 in all directions, maintaining a uniform grid of 256³ grid points. As discussed by Ono et al. (2013) this approach is reliable because it does not introduce errors larger than 0.1% after 40 remappings. In such a way, the finest spatial resolution is ~ 10¹³ cm at the beginning of our 3D simulations, and because of the expansion of the system the resolution is ~ 10¹⁷ cm at the end of simulation time. The simulations therefore span from an initial domain of (2.19 × 10¹⁵ cm)³ up to (~5.8 × 10¹⁹ cm)³. Equatorial symmetric boundary conditions¹ with respect to x = 0, y = 0 and z = 0 planes are considered, while fixed pre-shock values conditions are assumed at all other boundaries. We checked the effects of spatial resolution on the remnant evolution in Appendix A, and verified that our main conclusions are not significantly affected by the resolution adopted (see Appendix A for further details).

The ejecta distribution is expected to be characterized by small-scale clumping of material and larger-scale anisotropies (e.g., Orlando et al. 2016). Therefore, following Orlando et al. (2012), we introduced small-scale perturbations in the initial ejecta distribution, mimicking the small-scale clumping of ejecta. These inhomogeneities were modeled as per-cell random-density perturbations derived from a power-law probability distribution. All the clumps have an initial size of ≈ 10¹³ cm and the density perturbation of each clump is calculated as the ratio of the mass density of the resulting clump to the local (unperturbed) average density in the region occupied by the clump; the maximum density perturbation allowed is a factor of 4. In Sect. 2.3 we discuss how post-explosion large-scale anisotropies in the ejecta distribution are modeled.

We describe the CSM around the progenitor star as a spherically symmetric steady wind with a gas density proportional to r⁻² (where r is the radial distance from the center of explosion): $${\rho }_{\text{wind}}={\rho }_{\text{ref}}{\left(\frac{r_{\text{ref}}}{r}\right)}^2 .$$(7)

We assumed that the wind is isobaric with pressure P_wind ~ 10⁻¹⁰ dyn cm⁻² to prevent the wind from evolving towards the boundaries of the domain, causing “jumps” in the profile of the wind whenever the domain is extended. This assumption does not affect the results because the forward shock propagates much faster than the stellar wind during the whole evolution.

The CSM is characterized by three parameters: the reference density, the radial distance of the reference density, and the temperature (n_ref = 1 cm⁻³, r_ref = 2.5 pc, and T_ref = 100 K, respectively). The values considered are within the range of typical values for a progenitor RSG, which is subject to high mass-loss rates (10⁻⁵−10⁻⁴ M_⊙ yr⁻¹) and very slow winds (20–100 km s⁻¹) (Dwarkadas 2005; Crowther 2001). Assuming Ṁ = 5 × 10⁻⁵ M_⊙ yr⁻¹ and v_w = 20 km s⁻¹, we derived the reference mass density as ${\rho }_{\text{ref}}=\dot{M}/(4\pi v_w{r}_{\text{ref}}^2)$ which corresponds to a particle number density n_ref ≈ 1 cm⁻³. In such a tenuous environment, the radiative losses would become relevant after t_rad > 10⁴ yr (see Blondin et al. 1998; Orlando et al. 2005; Miceli et al. 2006; Lee et al. 2015), a timescale much greater than the simulated time (t = 5000 yr).

The simulations include the effect of an ambient magnetic field. This field is not expected to influence the overall dynamics of the system which is characterized by high values of the plasma β (i.e., the plasma pressure dominates over the magnetic one). We note that MHD codes are generally more diffusive than HD codes because they use approximate Riemann solvers. However, we decided to include the effects on an ambient magnetic field in our simulations because the field can play a relevant role in limiting the growth of HD instabilities that would develop at the border of ejecta clumps and anisotropies and that are responsible for their fragmentation (Fragile et al. 2005; Shin et al. 2008; Orlando et al. 2008, 2012, 2019). Following Orlando et al. (2019), the field configuration adopted is the “Parker spiral”, which is the field resulting from the rotation of the progenitor star and theexpanding stellar wind (Parker 1958). The Parker spiral can be described in spherical coordinates (r, θ, ϕ) as: $$\begin{array}{ll}{B}_{\text{r}}\hfill & =\frac{A_{\text{1}}}{r^2}\hfill \\ B_{\varphi }\hfill & =-\frac{A_{\text{2}}}{r}\sin \theta .\hfill \end{array}$$

The previous parameters are set to A₁ = 3 × 10²⁸ G cm² and A₂ = 8 × 10¹⁰ G cm in order to produce a magnetic field of the order of 10 G at the surface of the progenitor star (e.g., Tessore et al. 2017 and references therein).

Fig. 2 Density (left-hand quadrants) and temperature (right-hand quadrants) sections (x, 0, z) showing examples of the initial conditions. The values are color coded according to the scale shown for each quadrant. Upper quadrants show the case of a spherically symmetric explosion, and lower quadrants show a case with a dense, isobaric spherical anisotropy (namely run Si-R5-D750-V5 in Table 1).

2.3 Post-explosion large-scale anisotropy

We considered the case of a post-explosion large-scale anisotropy located in the inner initial ejecta distribution. Our aim is to study whether or not and how the remnant structure and morphology keep an imprint or “memory” of post-explosion anisotropies and whether or not these anisotropies have an effect on the evolution and mixing of the different metal-rich layers of the progenitor star. The anisotropy included in our model does not arise from the 1D SN simulation, but is set in the initial conditions of our SNR model. Explaining the physical origin of these initial anisotropies is well beyond the scope of this work; indeed we aimed to study the effects of a large-scale post-explosion anisotropy on the evolution of a SNR and to investigate how the different initial parameters of the anisotropy might determine the final shape of the remnant.

The anisotropy is described as an overdense sphere in pressure equilibrium with the surrounding ejecta whose center lies on the z-axis (therefore only a quarter of the sphere is modeled). We investigated possible numerical effects which may develop by assuming the overdense sphere propagating along the z-axis by exploring a case of a clump with its center located at ~45° in the plane (x, 0, z; see Appendix B for further details). We find that simulations assuming the sphere propagating either along the z-axis or at ~ 45° in the plane (x, 0, z) produce very similar results on the timescale considered (5000 yr) and we concluded that the numerical effects (if any) due to the propagation of the bullet along the z-axis have a negligible impact on the evolution of the ejecta bullet.

The geometrical properties of the anisotropy are param- etrized by its distance from the center (D), and its radius (r). We considered two values of D, associated with different layers of the exploding star, namely: the interface between ⁵⁶Fe/²⁸Si-rich regions at D = 0.24 R_SNR (where R_SNR is the initial radius of remnant) and the ⁵⁶Fe-rich region at D = 0.15 R_SNR. We explored two values for the clump radius, namely r = 0.05 R_SNR and r = 0.04 R_SNR, both consistent with the characteristic size of clumps generated by Rayleigh-Taylor (RT) instability (seeded by flow structures resulting from neutrino-driven convection) in SN explosion simulations (e.g., Kifonidis et al. 2006; Gawryszczak et al. 2010). We explored anisotropies with density between 500 and 750 times that of the surrounding ejecta at distance D (density contrast χ_n), and with radial velocity between two and seven times that of the surrounding ejecta (velocity contrast χ_v), making surethat the maximum speed of the clump does not exceed the velocity of the fastest ejecta in the outermost layers. Wongwathanarat et al. (2017) showed that fast and high-density large-scale anisotropies/plumes with high concentration of ⁵⁶Ni and ⁴⁴Ti can emerge from neutrino-driven SN explosions. Some other authors have predicted that high-velocity ejecta components might be produced in the presence of magnetic fields that are anchored to a rapidly rotating proto-neutron star (e.g., Takiwaki et al. 2009). The anisotropies parametrized in this paper assume the same chemical composition as the surrounding ejecta region, i.e., iron alone or iron and silicon when D = 0.15 and D = 0.24, respectively.A summary of the cases explored is given in Table 1.

Previous simulations of ejecta bullets in a SNR performed by Miceli et al. (2013) and Tsebrenko & Soker (2015) showed that clumps with initial values of χ_n ≤ 100 are able to reach the SNR shock front without being dispersed by the interaction with the reverse shock. At odds with those simulations which adopt power-law density profiles for the ejecta distribution and start tens of years after the SN event, our simulations include a self-consistent description of the initial ejecta density profile (derived from 1D SN simulations) which shows a high-density thick shell in the innermost regions of the ejecta (see left boxes of Fig. 2). As a result, we find that higher density contrasts (χ_n > 500) are necessary for the clumps to overcome this overdense structure and to eventually protrude the forward shock. The maximum initial mass of the anisotropies considered is ~0.25 M_⊙ (<2% of the initial ejecta mass) with a maximum kinetic energy of ~4 × 10⁴⁹ erg (<3% of the total initial kinetic energy of the ejecta).We note that these values are very similar to those estimated for the Si-rich jet-like feature observed in Cassiopeia A (Orlando et al. 2016) and to those inferred for the shrapnel G in the Vela SNR (García et al. 2017).

3 Results

3.1 The case of a spherically symmetric explosion

We first considered the case of a spherically symmetric SN explosion with small-scale (isotropically distributed) clumps of ejecta, without including any large-scale anisotropies (hereafter, the spherical model) to study the spherically symmetric expansion of the remnant and the evolution of its chemically homogeneous layers.

Figure 3 (see also online Movie 1 and 2) shows the 2D cross-sections through the (x, 0, z) plane of density at different evolutionary stages for the spherical model (left panel). After 100 yr of evolution, it is possible to see HD instabilities (Rayleigh-Taylor, Kelvin-Helmholtz shear instability; e.g., Gull 1973; Chevalier 1976; Fryxell et al. 1991; Chevalier et al. 1992) developing as finger-like structures at the interface between ejecta and shocked CSM (as already shown by numerical simulations in Dwarkadas & Chevalier 1998; Wang & Chevalier 2001), thus driving the mixing between the ejecta and the ambient surrounding material. During this evolutionary phase, the over-dense internal shell in the ejecta expands and is clearly visible in the density distribution. At the age of 1000 yr, the reverse shock reaches the over-dense internal shell. This allows the over-dense shell to enter in the mixing region causing the destruction of the shell itself by the HD instabilities at t = 5000 yr (see online Movie 1). At the end of the simulation, the SNR reaches a radius of ≈ 4 × 10¹⁹ cm (≈ 13 pc).

Among all the elements considered in the simulations, radioactive nuclei (such as ⁵⁶Ni and ⁴⁴Ti) synthesized during the explosion and their decay products are very important because they are produced in the innermost regions of the exploding star, which makes them excellent probes of the internal conditions that lead to the determination of the explosion mechanism and of the shock-wave dynamics during the earliest phases of SN outbursts. For this reason, the study of these elements is key to connecting explosion models to observations and to deducing important constraints on the underlying processes. Radioactive ⁵⁶Ni and ⁴⁴Ti transform into the stable isotopes ⁵⁶Fe and in ⁴⁴Ca, respectively, with a half-life of <100 yr (Nadyozhin 1994; Ahmad et al. 2006). As mentioned in Sect. 2, since we are interested in investigating the matter mixing and evolution of the remnant at epochs longer than 100 yr, and recalling that our model does not include the radioactive decay of the elements (see Sect. 2.2), we assume that all the ⁵⁶Ni and ⁴⁴Ti are already decayed to ⁵⁶Fe and ⁴⁴Ca, respectively,at the beginning of our simulations.

Figure 4 displays the ejecta mass distribution of selected elements versus the radial velocity at different epochs during the remnant expansion. The elements appear stratified: ⁵⁶Fe is concentrated in the slower and innermost part of the ejecta, approximately at 30% of the radius of the initialremnant, inside the dense shell seen in the density distribution in Fig. 2, along with the other heavy elements (as ⁴⁴Ca). ²⁸Si envelopes the inner ⁵⁶Fe, whereas ¹⁶O and ¹²C are prominent at the bottom of the ⁴He shell. The outer half of the radial ejecta distribution is dominated by ⁴He and ¹H.

By looking at the mass distribution of the selected elements at various times, we observe how the overall composition of the ejecta changes during expansion of the remnant (see Fig. 4). The overall structure and composition of the unshocked ejecta do not change significantly over the timescale considered, although some small differences in the innermost layers are present due to some numerical diffusivity of the code. This indicates that the unshocked ejecta expand almost homologously until they start to interact with the reverse shock. This is shown in Fig. 4, where the profiles of the elements with velocity <1500 km s⁻¹ remain almost unchanged at all epochs. In the unshocked ejecta, we observe a large fraction of ⁵⁶Fe and ⁴⁴Ca, together with the low-velocity tail of the ²⁸Si distribution. Instead, the profiles of the elements at a greater distance from the center of the explosion, that is, the elements with high-velocity components, are strongly influenced by their interaction with the reverse shock and the CSM. These elements are strongly decelerated in the intershock region and consequently mixed by HD instabilities. This produces a homogenization of the different profiles in these regions: the overall shapes of the high-velocity tails of the mass distributions of shocked elements are more similar to each other and characterized by steep slopes, suggesting a significant mixing between layers of different chemical composition due to HD instabilities developing at the contact discontinuity (see also Wongwathanarat et al. 2017 for the effects of HD instabilities on the mass distribution versus the radial velocity).

Nevertheless, it is important to note that the chemical distribution generated during the immediate aftermath of the SN (i.e., the onion-skin chemical layering) is roughly preserved after 5000 yr of evolution. Thus, for a spherically symmetric SN explosion, it is not possible to reproduce the inversion among the innermost chemical layers (e.g., Wang & Chevalier 2002; Wang 2011; Orlando et al. 2016), as observed for example in Cassiopeia A (e.g., Hughes et al. 2000), if the ejecta are characterized only by small-scale clumping.

Fig. 3 Density distributions in the (x, 0, z) plane at different simulation times for a spherically symmetric explosion (left panel) and for model Fe-R5-D750-V6 of Table 1 (right panel). From top left clockwise: t = 10, t = 100, t = 1000 and t = 5000 yr from the explosion. The units in the color bars are g cm−3 logarithmically scaled. The color-coded density scale is shown close to each quadrant. We highlight the different scales used in each quadrant for both the density and the distance along the axis. Right panel: the contours enclose the computational cells consisting of the original anisotropy material by more than 50% (solid line) and 10% (dotted line). See also online Movie 1 and 2 for an animated version.

Fig. 4 Mass distributions of selected elements as a function of radial velocity for the spherically symmetric explosion at the labeled years from the explosion. Upper left panel: initial conditions from Ono et al. (2020). Only the dominant fractions are considered. Mi is the total ejecta mass of element i, ΔMi is the mass of the ith element in the velocity range between v and v +δv. The size of the velocity bins δv is 100 km s−1. The black dashed vertical line shows the approximate reverse shock position in each panel.

3.2 Effect of large-scale anisotropy

As seen in Sect. 3.1, the ideal case of a spherically symmetric explosion cannot reproduce the complex morphology of an evolved SNR, which is in good agreement with the fact that most of the CC-SNe are believed to explode asymmetrically (Maeda et al. 2008; Wongwathanarat et al. 2013; Janka et al. 2016; Janka 2017; O’Connor & Couch 2018; Burrows et al. 2019). The physical characteristics of the SN explosion generate asymmetric expansion and strong inhomogeneities that have a complex influence on the evolution of the SNR. Here we consider a possible example of a dynamical event that can cause strong mixing and overturning of the chemical layers in the ejecta. We studied the effects of a post-explosion anisotropy generated soon after the shock breakout. Multiwavelength observations show the presence of clump structures in the ejecta of CC SNRs (Aschenbach et al. 1995; Miceli et al. 2008; García et al. 2017; Park et al. 2004; Katsuda et al. 2008; Fesen et al. 2006). We model the anisotropy as a spherical clump of radius r, with its center located on the z axis, at a distance D from the origin, as explained in Sect. 2.3.

As shown in Table 1, we explored the case of a clump located either at D = 0.15 or at D = 0.24, corresponding to the initial ⁵⁶Fe-rich internal layer and between the ⁵⁶Fe- and ²⁸Si-rich layers, respectively (see the top left panel in Fig. 4).

The right panel of Fig. 3 shows a 2D cross section through the (x, 0, z) plane of density distribution of the model Fe-R5-D750-V6 at four different evolution times (but qualitatively the overall evolution is similar for all the cases explored). The top-left quadrant shows the system 10 yr after the beginning of the simulation, when the clump has already overcome the dense shell initially located at about 30% of the radius of the ejecta and the topmost part of the clump has reached the reverse shock. At this stage, the clump has an elongated structure along the z-axis, from the central region of the remnant to the reverse shock. Only the head of the clump is still significantly overdense with respect to the surrounding ejecta. The top-right quadrant shows the system at 100 yr. The clump has not changed its morphology, but its topmost part interacts with the intershock region. Rayleigh-Taylor instabilities start to develop at the interface between ejecta and shocked CSM. After 1000 yr of evolution (bottom right quadrant) the part of the clump in the intershock region is partially eroded by the HD instabilities, while the passage of the clump has left a low-density strip region in the inner part of the ejecta. This phase lasts ≈ 2000 yr; the wake behind the ejecta clump is then filled in eventually as the pressure gradient drives lateral expansion of the ejecta to fill the void. For model Fe-R5-D750-V7 and Si-R5-D600-V4, we find that at ~ 2500 yr, the clump reaches a velocity of ~3000 km s⁻¹, similarly to what was found very recently for the runaway knot in N132D (3650 km s⁻¹, Law et al. 2020). The bottom-left quadrant shows the clump at 5000 yr and it is possible to see the characteristic supersonic bow shock of the clump protruding beyond the SNR forward shock. At this evolutionary stage, the low-density region has been filled by the surrounding material. An exception to this trend is represented by models Fe-R5-D750-V3 (see Fig. 5) and Fe-R5-D500-V3. In these two cases, the clump penetrates into the high-density shell without overcoming it, remaining stuck inside it. In this case, the clump expands together with the surrounding ejecta and eventually becomes completely fragmented by the interaction with the reverse shock.

Figure 6 (see also online Movie 3, 4, 5 and 6) shows the spatial distribution of ⁵⁶Fe (red), ²⁸Si (green) and ¹⁶O (blue) in the ejecta of the remnant at t = 5000 yr for different model setups (see Table 1). Figure 6A shows four models in which the large-scale clump is initially located at D = 0.24, that is, between the ²⁸Si and the ⁵⁶Fe shells, with r = 0.05. In each of these cases, the clump can produce a remarkable jet-like structure of ²⁸Si-rich ejecta. This structure protrudes beyond the SNR forward shock in all cases, except for model Si-R5-D500-V3 (Fig. 6A-1). We point out that the ²⁸Si-rich collimated structure has overcome the ¹⁶O-rich layer of ejecta at this evolutionary stage. Moreover, we find that after ≈ 350 yr, for models Si-R5-D500-V4 (Fig. 6A-2) and Si-R5-D600-V4 (Fig. 6A-4), this Si-rich jet-like structure has already overcome the ¹⁶O layer, which is similar to a feature observed in the northeast and southwest regions of Cassiopeia A (e.g., DeLaney et al. 2010; Milisavljevic & Fesen 2013).

Our simulations show that a velocity contrast χ_v > 2 is necessary to allow the clump to overcome the high-density shell. The length of the jet-like feature depends on the initial density contrast of the clump, as can be seen in quadrant (1) of Fig. 6A (χ_n = 500, χ_v = 3), where the position of the bow-shock of the clumphas its lowest value (4.12 × 10¹⁹ cm; ≈ 13 pc) among the cases considered here. The shape of the jet-like structures appears somehow “curly”, with the exception of Si-R5-D600-V3 (quadrant 3 of Fig. 6A), which has a more regular shape and is more collimated along the z-axis. The clump has been strongly eroded at this time.

The final distribution of ⁵⁶Fe, ²⁸Si, and ¹⁶O in the remnant in the case of D = 0.15 (clump entirely made of ⁵⁶Fe) and r = 0.05 is shown in Figs. 6B, C. In both cases, the clump produces a ⁵⁶Fe-rich region along the z-axis, extending from the center of the remnant to the intershock region where it mixes with the ¹⁶O layer. We find that the clump can protrude the forward shock only if its velocity contrast is χ_v ≥ 6, which is a factor of two higher than that obtained for D = 0.24. This is because, for D = 0.15, the clump is located in an inner and slower region (due to the linearly increasing velocity profile of the ejecta, described in Sect. 2), and therefore a higher velocity contrast is needed with respect to the case with D = 0.24. Moreover, for D = 0.24, the clump is at the inner border of the high-density shell, and so the interaction with the shell occurs at the early stages of evolution when the clump-to-shell mean density ratio (in the following η) is η ≥ 14. On the other hand, for D = 0.15, the clump evolves in a lower pressure region, and therefore the clump undergoes a significant expansion before interacting with the dense shell and η ≤ 0.25 at the interaction. This effect causes the clump to struggle to overcome the dense shell, leading to a wider jet-like structure with respect to cases with D = 0.24. We also find that the initial clump size (r) does not affect the mixing of ⁵⁶Fe in the outer layers (see Fig. 6D) within the range of values explored.

In the cases explored, we observed that the clump pushes the heavy elements (especially ⁵⁶Fe and ²⁸Si) into the outer and high-speed regions of the ejecta, as has been observed in other core-collapse SNe such as SN1987A (Haas et al. 1990; Utrobin et al. 1995), Cassiopeia A (Grefenstette et al. 2014, 2017; Siegert et al. 2015), Vela (García et al. 2017), and, very recently, N132D (Law et al. 2020). According to the exploration of the parameter space we have performed, the average chemical stratification of the ejecta is roughly preserved in regions not affected by the large-scale anisotropy (e.g., x-axis in Fig. 6). In these regions, the evolution of the ejecta is similar to that described in Sect. 3.1 and the forward shock reaches approximately the same distances in all the cases explored (the differences shown in quadrants (2) and (4) of Figs. 6B and C are due to slightly different simulation times, of the order ≈ 10 yr).

Shock-heated plasma emits X-ray thermal emission, especially through line emission, giving us a valuable tool to measure the elemental composition and spatial distribution inside the remnant. From the simulations, we derived the mass of shocked ⁵⁶Fe and ²⁸Si normalized to the value of the total mass of the respective element for the different models as a function of time (see Figs. 7–8). For both D = 0.15 and D = 0.24, the final amount of shocked ⁵⁶Fe is at least two times larger than that obtained for the case of a spherically symmetric explosion.

For D = 0.15 (see Fig. 7), during the first 300 yr, the amount of shocked ⁵⁶Fe depends both on the size and on the initial density contrast of the clump, increasing with size and contrast. However, after 3000 yr, the reverse shock interacts with the ⁵⁶Fe shell even without any initial clump, and therefore the amount of shocked ⁵⁶Fe observed at the end of the simulations is not entirely due to the initial clump. As for the shocked ²⁸Si, we find that the final amount of shocked mass does not show any strong dependence on the initial clump parameters or on its initial dimension or position. The main cause of this is the interaction of the reverse shock with the Si-rich shell, which shocks a relevant part of the ²⁸Si-rich material regardless of the presence of the initial clump, thus making the contribution of the clump negligible. For D = 0.24 (see Fig. 8) on the other hand, the amount of shocked ⁵⁶Fe after 5000 yr is not strongly affected by the initial clump parameters. This is expected because the jet-like structure is mostly made of ²⁸Si for models with D = 0.24, as shown in Fig. 6. On the other hand, the velocity contrast χ_v plays an important role in determining the beginning of the interaction between the reverse shock and the ⁵⁶Fe-rich ejecta. We observed that, for χ_v = 4, the reverse shock heats the ⁵⁶Fe-rich material ~90 yr earlier than the case with χ_v = 3 (see upper panel of Fig. 8). However, as seen for the ⁵⁶Fe, we find that the initial velocity contrast χ_v affects the timing of the interaction between the reverse shock and the ²⁸Si-rich ejecta, but only for D = 0.24.

To deepen our understanding of the role played by the clump, we also computed the amount of shocked ²⁸Si (upper panel of Fig. 9) and ⁵⁶Fe (lower panel of Fig. 9) as a function of time in three different regions of the ejecta, namely: (i) cells within 45° of the z-axis, that is, those enclosing the ejecta clump and mostly affected by its evolution; (ii) cells within 15° of the equatorial plane, which is those least affected by the clump; and (iii) all the other cells. The shocked mass is normalized to the total mass of the respective element in the corresponding region. We compared the spherically symmetric explosion (green lines) and model Fe-R5-D750-V7 (red lines, see Table 1). For model Fe-R5-D750-V7, in the region within 45° of the z-axis the mass of the ²⁸Si (upper panel) has been shocked by a larger fraction compared to the other regions, but only up to t = 4000 yr. At odds with the ⁵⁶Fe (lower panel), after t = 4000 yr the reverse shock reaches the ²⁸Si layer in the regions away from the jet-like structure making these regions the main contributors to the shocked mass of ²⁸Si. Looking at the evolution of the shocked mass of ⁵⁶Fe (lower panel) for model Fe-R5-D750-V7, the effect of the clump stands out in cells within 45° of the z-axis, while the other regions follow the same evolution as the spherically symmetric explosion. However, the presence of the clump influences the regions within 45° of the equatorial plane shocking 40% more mass of the region compared to the spherically symmetric explosion.

Fig. 5 Same as the right panel of Fig. 3 but for model Fe-R5-D750-V3 (see Table 1).

Fig. 6 Color-coded images of the logarithm of the ejecta mass fraction (>10−4) distributionsat the end of the simulation (t ≈ 5000 yr) for different models (the model as presented in Table 1 is reported near each quadrant) in the (x, 0, z) plane of 56Fe (red), 28Si (green), and 16O (blue). Black contours enclose the computational cells consisting of the original clump material by more than 50% (solid line) and 10% (dotted line). The white dashed line represents the projected position of the forward shock; see also online Movie 3, 4, 5, and 6 for an animated version.

Fig. 7 Mass of shocked 28Si (upper panel) and 56Fe (lower panel) vs. time for models assuming a clump initially located at D = 0.15 and characterized by different parameters (see Table 1). The shocked mass is normalized to the total mass of the relative element.

3.3 Synthesis of X-ray emission

In this paper we show that the presence of an anisotropy can significantly influence the shape of the remnant leading to the formation oflarge-scale jet-like structures and protrusions of the remnant outline. The matter-mixing strongly depends on the physical andchemical characteristics of the initial anisotropy and it is now interesting to have a probe of such mixing by looking at the X-rayspectra. When the shock wave interacts with the matter, the latter is heated up to temperatures of several million degrees, leading to X-ray emission. Thus, X-ray spectral analysis is a powerful tool to distinguish between the different spectral signatures of the regions affected by an ejecta bullet and the rest of the evolving remnant between the forward and the reverse shock.

By adopting the approach described in Greco et al. 2020 (see also Miceli et al. 2019), we self-consistently synthesized X-ray spectra from our simulations in a format virtually identical to that of Chandra ACIS-S observations. The synthesis includes the deviations from equilibrium of ionization and from temperature-equilibration between electrons and protons, and takes into account the chemical composition of ejecta and ISM resulting from the simulations in each computational cell (see e.g., Miceli et al. 2019; Orlando et al. 2019 for more details). For the spectral synthesis, we assume a distance of 1 kpc and a column density N_H= 2 × 10²¹ cm⁻² with an exposure time of 100 ks. As an example, we analyzed run Fe-R5-D750-V7 (see Table 1) where a clear jet-like feature is evident at an age of t = 5000 yr. We selected three regions, identified by the red, black, and green boxes in the upper panel of Fig. 10, corresponding to the ejecta bullet protruding the remnant outline (in the following the clump), the wake (namely the low-density strip-region behind the clump, already mentioned in Sect. 3.2), and an area of the shell not affected by the anisotropy, respectively. For a distance from the observer of 1 kpc, each of these regions corresponds to a box with angular side l = 2′. Synthetic spectra extracted from the three regions are shown in the lower panel of Fig. 10.

The main spectral difference between the spectrum extracted from the region enclosing the clump and the other spectra is the bright “false-continuum” emission which is due to the blending of Fe XIV-XXIV emission lines at energies around 1 keV. This result showsthat the clump spectrum is dominated by the emission of Fe-rich ejecta. Emission lines of other elements such as Mg (~ 1.3 keV) and Si (~1.8 keV) are fainter in the clump spectrum, because the protrusion of the Fe-rich clump drags away lighter elements. Enhanced Ne (at ~ 0.92 keV), Mg, and Si emission lines are visible in the spectrum of the shell region (green spectrum in Fig. 10). The wake spectrum (black spectrum in Fig. 10) shows bright Mg and Si emission lines, very low Fe emission (even lower than that in the shell region), and an overall lower emission measure. The emission is fainter because the density is low, as we can see in the integrated density distribution map displayed in the upper panel of Fig. 10. Our analysis shows that the effects of matter mixing in the ejecta produced by the evolution of a post-explosion anisotropy lead to spectralfeatures in the X-ray band which are detectable through a spatially resolved analysis of X-ray spectra of SNRs.

Fig. 8 Same as Fig. 7 but for models with a clump initially located at D = 0.24.

Fig. 9 Mass of shocked 28Si (upper panel) and 56Fe (lower panel) vs. time for the spherically symmetric explosion (green lines) and for model Fe-R5-D750-V7 (red lines, see Table 1). Solid and dashed lines show the shocked mass from cells within 45° of the z-axis and from 15° of the equatorial plane respectively, while the dotted lines show the mass from the others cells. The shocked mass is normalized to the total mass of the relative element in the region considered.

Fig. 10 Upper panel: density distribution map integrated along the line of sight (in this case the y-axis) of the run Fe-R5-D750-V7, showing the regions chosen for the spectral synthesis: red for the anisotropy protruding the remnant outline (the clump), black for the wake of the anisotropy, and green for a region of the shell not affected by the anisotropy (the shell). Lower panel: X-ray Chandra ACIS-S synthetic spectra in the 0.5-2 keV energy band for the clump (red), the wake (black), and the shell (green).

4 Summary and conclusions

We investigated the evolution of the ejecta in core-collapse SNRs from the onset of the SN to the full-fledged remnant. The aim was to analyze how matter mixing occurs during the remnant expansion and interaction with the CSM, how the various chemically homogeneous layers at the time of the explosion map into the resulting abundance pattern observed when the remnant is fully developed, and how the presence of large-scale anisotropies affects the evolution of these layers. To this end, we developed a 3D MHD model describing the evolution of a SNR starting soon after the SN explosion and following the interaction of the SN ejecta with the CSM for 5000 yr. As initial conditions for the SNR simulations, we adopted the output of a spherical SN explosion – model by Ono et al. (2020) (see also Ono et al. 2013) – soon after the breakout of the shock wave at the stellar surface (about 17 h after the core-collapse). More specifically, we considered the case of a 19.8 M_⊙ progenitor RSG which exploded as a SN withan injected energy of 2.5 × 10⁵¹ erg. The initial ejecta structure was modeled by including small-scale clumping of material and larger-scale anisotropy. The simulations are multi-species to trace the life cycle of elements during the whole evolution, from the nuclear reaction network of the SN to the enrichment of the CSM through mixing of chemically homogeneous layers of ejecta.

To ascertain the effect of the spherically symmetric expansion of the remnant on the matter-mixing, we first explored a case with small-scale clumpy structures of ejecta without considering any large-scale anisotropy. The ejecta soon after the SN explosion show the original onion-like structure, with elements arranged in concentric shells. Although the onion-like structure is maintained after 5000 yr of evolution, which is due to the homologous expansion of the unshocked layers of the ejecta, we observed a significant matter mixing in the inter-shock region of the remnant, as also expected from previous studies (e.g., Gull 1973; Chevalier 1976; Fryxell et al. 1991; Chevalier et al. 1992; see Fig. 4). This mixing is caused by the development of HD instabilities generated by the interaction of the ejecta with the CSM. However, this mixing is not able to cause a spatial inversion between chemical layers, leaving the initial chemical configuration of the ejecta roughly preserved.

We then performed a set of simulations to study the evolution of post-explosion large-scale anisotropy (the clump) in the ejecta and its effects on the matter-mixing, investigating the role of its initial parameters (position, dimension, density, and velocity contrasts). We primarily explored two cases: (1) a clump initially located in the ⁵⁶Fe/²⁸Si region (D = 0.24) and (2) a clump in the ⁵⁶Fe region (D = 0.15). In both cases, the clump is located behind a thick and dense shell of lower-Z ejecta. For this thick and dense shell, a higher density contrast is required for the clump to reach the forward shock compared to previous works (Miceli et al. 2013; Tsebrenko & Soker 2015). The dense-shell structure could strongly depend on the progenitor model. This is because different progenitors could produce different density profiles after the explosion, and therefore we should obtain different initial conditions and the clump would evolve in a different environment. The progenitor discussed in this paper (a 19.8 M_⊙ progenitor RSG) is a generic case of a RSG among the different progenitors simulated in Ono et al. (2020). The magnetic field is not relevant for the overall evolution of the remnant (the value for plasma-β being much larger than 1) but plays some role in reducing the fragmentation of the clump caused by the HD instabilities that would develop at its border (Fragile et al. 2005; Shin et al. 2008; Orlando et al. 2008, 2012).

From the first ≈10 yr of evolution, the clump rapidly assumes an elongated structure and its outermost region becomes denser. It takes ~100 yr for the clump to overcome the high-density shell of lower-Z ejecta and to start to interact with the reverse shock and the intershock region. After ~1000 yr, the HD instabilities begin to erode the clump by gradually fragmenting it, supported by the effects of compression and heating caused by the interaction of the clump with the reverse shock. We find that the initial position of the clump strongly affects the morphology of the clump at t = 5000 yr. A clump initially located at D = 0.24, that is, attached to the high-density shell (see right quadrant of Fig. 2), entirely overcomes the high-density shell, and at the end of the simulation only a small fraction of the clump is left behind the reverse shock (e.g., top left panel of Fig. 6 where contours enclosing cells with more than 10% of the original clump material are shown). Nevertheless, a clump initially located at D = 0.15 expands significantly before interacting with the shell. This involves a minor density contrast of the clump at the interaction with the shell, and consequently a major difficulty to overcome the shell. This translates into a higher density contrast of the clump required to overcome the high-density shell of lower Z ejecta and a larger fraction of the clump left behind the reverse shock at the end of simulation time (see top right and bottom left panels of Fig. 6).

Furthermore, we find that, for D = 0.24, a clump is able to protrude the forward shock with an initial velocity contrast χ_v = 3, while for D = 0.15, χ_v > 6 is required. For a smaller clump (r = 0.04), only in model Fe-R4-D750-V7 (χ_n = 750, χ_v = 7) is the clump able to protrude the forward shock.

The chemical stratification in the ejecta is strongly affected by the clump, especially in the region along the z-axis. The initial onion-like structure is not preserved, as the clump propagating through the remnant forms a stream of ²⁸Si(⁵⁶Fe) ejecta, for D = 0.24 (D = 0.15), from the inner regions up to the intershock region. The clump makes its way out to the external layer of the remnant by piercing the chemically homogeneous shells and pushing the layers on the side of the stream. These streams cause a spatial inversion of the chemical layers, bringing the ²⁸Si/⁵⁶Fe to the external ¹⁶O shell (see Figs. 6 and 10). Where the effects of the clump’s passage are not relevant, that is, in the (x, y, 0) plane, the evolution of the chemical layer is analogous to that described by the spherically symmetric model (see Fig. 9). The main effect of this inversion is on the amount of shocked mass of both ⁵⁶Fe and ²⁸Si as shown in Figs. 8, 7, and 9. We find that a clump initially located in the ⁵⁶Fe/²⁸Si region does not increase the final amount of shocked ²⁸Si, whereas it doubles the amount of shocked ⁵⁶Fe. Although the density contrast (χ_n) does not produce large effects, the velocity contrast (χ_v) has an important role in determining the age of the shocking time of the elements: a higher contrast causes an earlier interaction. As a result, this will produce an earlier appearance of X-ray emission from shocked ejecta. In particular, we find that for a spherically symmetric expansion, the ²⁸Si (⁵⁶ Fe) begins to be shocked after ~2000 yr (~3000 yr), while the presence of the clump causes an early interaction of both elements with the reverse shock ~ 100 yr (at most) in advance. For a clump initially located in the Fe shell (i.e., D = 0.15), we find that the trend over time of the shocked ²⁸Si mass is not affected by the presence of the initial large-scale anisotropy, while on the contrary the final amount of shocked ⁵⁶Fe increases with increasing size, velocity, and density contrast of the clump.

Our simulations allowed us to study the effect of large-scale anisotropy on the evolution and matter-mixing of the ejecta of a SNR and provided a valuable tool to reproduce their observable properties. By synthesizing X-ray spectra from the model we find that it is possible to obtain more detailed diagnostics of ejecta inhomogeneities in CC-SNe by carefully comparing X-ray observations of SNRs with numerical simulations. Thus, our findings can be a useful guide in the interpretation of observations, such as for example the Si-rich shrapnel protruding beyond the front of the primary blast shock wave in the Vela SNR (García et al. 2017) or the fast runaway knot with a significantly high Si abundance recently detected in N132D (Law et al. 2020).

Movies

Movie 1 Access here

Movie 2 Access here

Movie 3 Access here

Movie 4 Access here

Movie 5 Access here

Movie 6 Access here

Acknowledgements

We thank an anonymous referee for the careful reading of the manuscript and for constructive and helpful criticism. We acknowledge that the results of this research have been achieved using the PRACE Research Infrastructure resource Marconi based in Italy at CINECA (PRACE Award N.2016153460). The PLUTO code is developed at the Turin Astronomical Observatory (Italy) in collaboration with the Department of General Physics of Turin University (Italy) and the SCAI Department of CINECA (Italy). The SN simulation has been performed with the FLASH code, developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago (USA). The numerical computations of the SN simulation were carried out complementarily on XC40 (YITP, Kyoto University), Cray XC50 (Center for Computational Astrophysics, National Astronomical Observatoryof Japan), HOKUSAI (RIKEN). S.O., M.M., F.B. and E.G. acknowledge financial contribution from the agreement ASI-INAF n.2017- 14-H.O, and partial financial support by the main-stream INAF grant “Particle acceleration in galactic sources in the CTA era”. This work is supported by JSPS Grants-in-Aid for Scientific Research “KAKENHI” Grant Numbers JP26800141 and JP19H00693. S.N., M.O. and G.F. wish to acknowledge the support from the Program of Interdisciplinary Theoretical & Mathematical Sciences (iTHEMS) at RIKEN (Japan). S.N. also acknowledges the support from Pioneering Program of RIKEN for Evolution of Matter in the Universe (r-EMU).

Appendix A Effect of spatial resolution

To ensure that the resolution used to perform the simulations (256³ grid zones) was sufficient to catch the basic properties of the system evolution, we repeated simulation Fe-R4-D750-V6 but using a uniform grid of 512³ grid zones (hereafter simulation Fe-R4-D750-V6-512pt) and we did the same for the spherically symmetric explosion simulation. The left panel of Fig. A.1 shows density distributions in the (x, 0, z) plane at the end of simulation time for both models at different resolutions. The right panel of Fig. A.1 shows the same as left but for color coded images of the mass fraction distributions of ⁵⁶Fe, ²⁸Si, and ¹⁶O. As expected, the spatial resolution mainly affects the structure of the RT instability: the RT fingers appear more extended and branched in the case of higher spatial resolution (see right quadrants in the panels of Fig. A.1). As a consequence, the spatial resolution changes the clump structure during its evolution. We also computed the distance of the center of mass of Fe-rich ejecta from the origin of the explosion for models Fe-R4-D750-V6 and Fe-R4-D750-V6-512pt, obtaining a value of 2.9 × 10¹⁹ cm and 3.2 × 10¹⁹ cm respectively, thus showing a slight discrepancy of the order of 10% in the distribution of ⁵⁶Fe in the two cases. We checked that simulations Fe-R4-D750-V6 and Fe-R4-D750-V6-512pt yield similar results in terms of the amount of shocked material. We were mainly interested in the mixing between different chemical layers, we therefore computed the amount of ⁵⁶Fe and ²⁸Si within the cells where the mixing is more relevant for both spherical (left panel of Fig. A.2) and nonspherical (right panel of Fig. A.2) models. The profiles shown in Fig. A.2 are in good agreement between the two models. This proves that, even if the small-scale density structures developing during the propagation of the ejecta bullet depend on the spatial resolution, a resolution of 256³ uniform grid zones is able to capture the basic properties of the system evolution and of the matter-mixing among different chemical layers.

Fig. A.1 Left panel: density distributions in the (x, 0, z) plane at the end of simulation time for a SN explosion with a large-scale anisotropy (upper quadrants) (runs Fe-R4-D750-V6 and Fe-R4-D750-V6-512pt in Table 1) and for a spherically symmetric explosion (lower quadrants), at low resolution (2563 grid points; left quadrants) and high-resolution (5123 grid points; right quadrants). Right panel: same but for color-coded images of the logarithm of the mass fraction distributions of Fe (red), Si (green), O (blue). In the upper quadrants of each panel, the black contours enclose the computational cells consisting of the original anisotropy material by more than 10% (dotted line)and 50% (solid line). The white dashed line represents the approximate position of the forward shock.

Fig. A.2 Left panel: amount of 56Fe and 28Si mass vs. time in cells with a mass fraction between 10 and 50% for a spherically symmetric explosion at different resolution. Right panel: amount of 56Fe and 28Si mass vs. time in cells with a mass fraction between 10 and 50% and with anisotropy material between 10 and 90% for models Fe-R4-D750-V6 (2563 grid points) and Fe-R4-D750-V6-512pt (5123 grid points) (see Table 1).

Appendix B Boundary effects

In our numerical setup, the large-scale anisotropy, or clump, has been described as a sphere with its center located on the z-axis of the domain (see Sect. 2.3). In Sect. 3.2 we show that the presence of this clump can cause a jet-like structure in the SNR that generates a strong mixing between different chemical layers. To ensure that the arise of this structure is not due to numerical artifacts caused by the proximity of the axis to the clump, we repeated simulation Fe-R4-D750-V6 placing the center of the clump at ≈45° in the (x,0,z) plane (hereafter run Fe-R4-D750-V6-ROT). In Fig. B.1 we compare both models after t ≈ 100 yr of evolution and at the end of the simulation, namely t ≈ 5000 yr after the explosion. The evolution and final density distribution does not change significantly rotating the initial position of the clump out from the z-axis. More specifically, we note that the evolution of the ejecta bullet is indistinguishable in the two cases during the first phase of the simulations (the first 30 yr), namely before the wake of the bullet reaches the (reflective) boundaries in run Fe-R4-D750-V6-ROT. In this latter case, the interaction of the wake with the boundaries produces reflection of material which propagates back into the domain. After ~100 yr (see upper panel in Fig. B.1), when the clump overcomes the high-density shell of ejecta and starts to interact with the inter-shock region, we observe a lesser presence of HD instabilities in the model Fe-R4-D750-V6-ROT (left panel) and also some artifacts due to the reflection of material on the boundary during the expansion. For this reason, the mixing of anisotropy material in the two models is slightly different after ~ 100 yr of evolution (see Fig. B.2). The perturbation of the inner portion of the remnant by the material reflected at the boundaries in run Fe-R4-D750-V6-ROT is the reason why we preferred to consider the bullet propagating along the z axis rather than at 45° angle from all axes. Nevertheless, even with the perturbations in run Fe-R4-D750-V6-ROT, the evolution of the ejecta bullet in both simulations is very similar, with no significant effects on the final remnant structure (see Sect. 3.2). In light of this, we can confirm that the jet-like structures observed in the simulations discussed in the paper do not arise from numerical artifacts.

Fig. B.1 Density distribution maps in the (x, 0, z) plane at t ≈ 100 (upper panels) and t ≈ 5000 yr from the explosion (lower panels). Right: model Fe-R4-D750-V6 (see Table 1). Left: same model with the initial position of the clump located at 45° in the (x, 0, z) plane, namely Fe-R4-D750-V6-ROT. Left panels: the (x, 0, z) plane from the simulation is rotated to the z-axis for comparison, and thus the lower blue triangle is not part of the computational domain. We note the different scales in the upper and lower panels. The contours enclose the computational cells consisting of the original anisotropy material by more than 50% (continuous line) and 10% (dotted line).

Fig. B.2 Mass vs. time (in logarithmic scale) for the two models compared in Fig. B.1. Total mass in cells with anisotropy material by more than 10% (continuous line) and in the range 10–90% (dotted line). In every case the mass is normalized to the value of the mass in the first step of the simulation.

References

All Tables

All Figures

	Fig. 1 Initial conditions for the 3D simulations. Top panel: radial profiles of initial ejecta density (black line) and pressure (red line). Bottom panel: radial profiles of initial ejecta velocity (black line) and Mach number (red line).
In the text

	Fig. 2 Density (left-hand quadrants) and temperature (right-hand quadrants) sections (x, 0, z) showing examples of the initial conditions. The values are color coded according to the scale shown for each quadrant. Upper quadrants show the case of a spherically symmetric explosion, and lower quadrants show a case with a dense, isobaric spherical anisotropy (namely run Si-R5-D750-V5 in Table 1).
In the text

Fig. 3 Density distributions in the (x, 0, z) plane at different simulation times for a spherically symmetric explosion (left panel) and for model Fe-R5-D750-V6 of Table 1 (right panel). From top left clockwise: t = 10, t = 100, t = 1000 and t = 5000 yr from the explosion. The units in the color bars are g cm−3 logarithmically scaled. The color-coded density scale is shown close to each quadrant. We highlight the different scales used in each quadrant for both the density and the distance along the axis. Right panel: the contours enclose the computational cells consisting of the original anisotropy material by more than 50% (solid line) and 10% (dotted line). See also online Movie 1 and 2 for an animated version.

In the text

Fig. 4 Mass distributions of selected elements as a function of radial velocity for the spherically symmetric explosion at the labeled years from the explosion. Upper left panel: initial conditions from Ono et al. (2020). Only the dominant fractions are considered. Mi is the total ejecta mass of element i, ΔMi is the mass of the ith element in the velocity range between v and v +δv. The size of the velocity bins δv is 100 km s−1. The black dashed vertical line shows the approximate reverse shock position in each panel.

In the text

	Fig. 5 Same as the right panel of Fig. 3 but for model Fe-R5-D750-V3 (see Table 1).
In the text

Fig. 6 Color-coded images of the logarithm of the ejecta mass fraction (>10−4) distributionsat the end of the simulation (t ≈ 5000 yr) for different models (the model as presented in Table 1 is reported near each quadrant) in the (x, 0, z) plane of 56Fe (red), 28Si (green), and 16O (blue). Black contours enclose the computational cells consisting of the original clump material by more than 50% (solid line) and 10% (dotted line). The white dashed line represents the projected position of the forward shock; see also online Movie 3, 4, 5, and 6 for an animated version.

In the text

	Fig. 7 Mass of shocked 28Si (upper panel) and 56Fe (lower panel) vs. time for models assuming a clump initially located at D = 0.15 and characterized by different parameters (see Table 1). The shocked mass is normalized to the total mass of the relative element.
In the text

	Fig. 8 Same as Fig. 7 but for models with a clump initially located at D = 0.24.
In the text

Fig. 9 Mass of shocked 28Si (upper panel) and 56Fe (lower panel) vs. time for the spherically symmetric explosion (green lines) and for model Fe-R5-D750-V7 (red lines, see Table 1). Solid and dashed lines show the shocked mass from cells within 45° of the z-axis and from 15° of the equatorial plane respectively, while the dotted lines show the mass from the others cells. The shocked mass is normalized to the total mass of the relative element in the region considered.

In the text

Fig. 10 Upper panel: density distribution map integrated along the line of sight (in this case the y-axis) of the run Fe-R5-D750-V7, showing the regions chosen for the spectral synthesis: red for the anisotropy protruding the remnant outline (the clump), black for the wake of the anisotropy, and green for a region of the shell not affected by the anisotropy (the shell). Lower panel: X-ray Chandra ACIS-S synthetic spectra in the 0.5-2 keV energy band for the clump (red), the wake (black), and the shell (green).

In the text

Fig. A.1 Left panel: density distributions in the (x, 0, z) plane at the end of simulation time for a SN explosion with a large-scale anisotropy (upper quadrants) (runs Fe-R4-D750-V6 and Fe-R4-D750-V6-512pt in Table 1) and for a spherically symmetric explosion (lower quadrants), at low resolution (2563 grid points; left quadrants) and high-resolution (5123 grid points; right quadrants). Right panel: same but for color-coded images of the logarithm of the mass fraction distributions of Fe (red), Si (green), O (blue). In the upper quadrants of each panel, the black contours enclose the computational cells consisting of the original anisotropy material by more than 10% (dotted line)and 50% (solid line). The white dashed line represents the approximate position of the forward shock.

In the text

Fig. A.2 Left panel: amount of 56Fe and 28Si mass vs. time in cells with a mass fraction between 10 and 50% for a spherically symmetric explosion at different resolution. Right panel: amount of 56Fe and 28Si mass vs. time in cells with a mass fraction between 10 and 50% and with anisotropy material between 10 and 90% for models Fe-R4-D750-V6 (2563 grid points) and Fe-R4-D750-V6-512pt (5123 grid points) (see Table 1).

In the text

Fig. B.1 Density distribution maps in the (x, 0, z) plane at t ≈ 100 (upper panels) and t ≈ 5000 yr from the explosion (lower panels). Right: model Fe-R4-D750-V6 (see Table 1). Left: same model with the initial position of the clump located at 45° in the (x, 0, z) plane, namely Fe-R4-D750-V6-ROT. Left panels: the (x, 0, z) plane from the simulation is rotated to the z-axis for comparison, and thus the lower blue triangle is not part of the computational domain. We note the different scales in the upper and lower panels. The contours enclose the computational cells consisting of the original anisotropy material by more than 50% (continuous line) and 10% (dotted line).

In the text

	Fig. B.2 Mass vs. time (in logarithmic scale) for the two models compared in Fig. B.1. Total mass in cells with anisotropy material by more than 10% (continuous line) and in the range 10–90% (dotted line). In every case the mass is normalized to the value of the mass in the first step of the simulation.
In the text