From the supernova to the supernova remnant: the three-dimensional imprint of a thermonuclear explosion

TL;DR

This study links 3D supernova explosion models to supernova remnant evolution, showing that initial asymmetries significantly influence remnant morphology observed centuries later.

Contribution

It introduces a method to connect supernova asymmetries with supernova remnant structures using combined 3D simulations.

Findings

Initial SN asymmetries affect SNR morphology after hundreds of years.

Previous symmetric models failed to reproduce observed large-scale structures.

Asymmetries from SN are crucial for realistic SNR modeling.

Abstract

Recent progress in the three-dimensional modeling of supernovae (SN) has shown the importance of asymmetries for the explosion. This calls for a reconsideration of the modeling of the subsequent phase, the supernova remnant (SNR), which has commonly relied on simplified ejecta models. In this paper we bridge SN and SNR studies by using the output of a SN simulation as the input of a SNR simulation carried on until 500~yr. We consider the case of a thermonuclear explosion of a carbon-oxygen white dwarf star as a model for a Type Ia SN; specifically we use the N100 delayed detonation model of Seitenzahl et al 2013. In order to analyze the morphology of the SNR, we locate the three discontinuities that delineate the shell of shocked matter: the forward shock, the contact discontinuity, and the reverse shock, and we decompose their radial variations as a function of angular scale and time. Assuming a uniform ambient medium, we find that the impact of the SN on the SNR may still be visible after hundreds of years. Previous 3D simulations aiming at reproducing Tycho's SNR, that started out from spherically symmetric initial conditions, failed to reproduce structures at the largest angular scales observed in X-rays. Our new simulations strongly suggest that the missing ingredient was the initial asymmetries from the SN itself. With this work we establish a way of assessing the viability of SN models based on the resulting morphology of the SNR.

Tables

Table 1. Table 1: Data for the radioactive decay of Ni 28 56 subscript superscript Ni 56 28 {}_{28}^{56}\mathrm{Ni} and Co 27 56 subscript superscript Co 56 27 {}_{27}^{56}\mathrm{Co} (from Jeffery, 1999 ) .

	${}_{28}{}^{56}\mathrm{Ni}$	${}_{27}{}^{56}\mathrm{Co}$
$t_{1/2}$ = half-life (days)	6.077(12)	77.27(3)
$t_e$ = $e$-folding time (days)	8.767(12)	111.48(4)
$Q$ total (MeV)	2.135(11)	4.566(2)
$Q$ photons+kinetic (MeV)	1.729(17)	3.74(4)
$C$ (erg.s-1.g-1)	$3.94\left(4\right)\times {10}^{10}$	$7.27\left(7\right)\times {10}^9$

Equations

$$r_{\mathrm{ch}}$$

$$u_{\mathrm{ch}}$$

$$t_{\mathrm{ch}}$$

$$\tilde{\mathbf{r}}$$

$$\mathrm{d}\tilde{t}$$

$$\tilde{\rho }$$

$$\tilde{\mathbf{u}}=\frac{\mathrm{d}\tilde{\mathbf{r}}}{\mathrm{d}\tilde{t}}=a^{\beta }\left(\mathbf{u}-H\mathbf{r}\right)$$

$$H=\frac{1}{a}\frac{\mathrm{d}a}{\mathrm{d}t}=\frac{\mathrm{d}\ln a}{\mathrm{d}t}$$

$$\tilde{P}$$

$$\tilde{E}$$

$$\widetilde{ϵ}$$

$$\frac{\partial \tilde{\rho }}{\partial \tilde{t}}+\tilde{\nabla }.\left(\tilde{\rho }\tilde{\mathbf{u}}\right)=(\alpha -\nu )\tilde{H}\tilde{\rho }$$

$$\frac{\partial \tilde{\rho }\tilde{\mathbf{u}}}{\partial \tilde{t}}+\tilde{\nabla }.\left(\tilde{\rho }\tilde{\mathbf{u}}\tilde{\mathbf{u}}\right)=(\alpha -\nu +\beta -1)\tilde{H}\tilde{\rho }\tilde{\mathbf{u}}-\tilde{\mathbf{F}}$$

$$\frac{\partial \tilde{E}}{\partial \tilde{t}}+\tilde{\nabla }.\left(\left(\tilde{E}+\tilde{P}\right)\tilde{\mathbf{u}}\right)=\left(\alpha -\nu +2(\beta -1)\right)\tilde{H}\tilde{E}-\tilde{\mathbf{F}}.\tilde{\mathbf{u}}+\left(\frac{2}{\gamma -1}-\nu \right)\tilde{H}\tilde{P}$$

$$\tilde{H}(t)=a^{2}H(t),$$

$$\tilde{\mathbf{F}}\left(\tilde{\mathbf{r}},\tilde{t}\right)=\tilde{f}\left(\tilde{t}\right)\tilde{\rho }\left(\tilde{\mathbf{r}},\tilde{t}\right)\tilde{\mathbf{r}}$$

$$\tilde{f}=-\left(\frac{\mathrm{d}\tilde{H}}{\mathrm{d}\tilde{t}}-\beta {\tilde{H}}^2\right)=-\left(a^{2\beta +1}\frac{{\mathrm{d}}^{2}a}{\mathrm{d}{t}^2}\right)$$

$$a(t)=a_{\star }{\left(\frac{t}{t_{\star }}\right)}^{\lambda }$$

$$\frac{\mathrm{d}a}{\mathrm{d}t}$$

$$\frac{{\mathrm{d}}^{2}a}{\mathrm{d}{t}^2}$$

$$\frac{\tilde{t}-{\tilde{t}}_0}{{\tilde{t}}_{\star }}=\{\begin{array}{ll}\frac{1}{1-2\lambda }\left({\left(\frac{t}{t_{\star }}\right)}^{1-2\lambda }-{\left(\frac{t_0}{t_{\star }}\right)}^{1-2\lambda }\right)& \lambda \ne \frac{1}{2}\\ \ln \left(\frac{t}{t_{\star }}\right)-\ln \left(\frac{t_0}{t_{\star }}\right)=\ln \left(\frac{t}{t_0}\right)& \lambda =\frac{1}{2}\end{array}$$

$$\lambda (t)=\frac{da/a}{dt/t}=\frac{d\ln (a)}{d\ln (t)}$$

$$H(t)=\frac{1}{a}\frac{da}{dt}=\frac{\lambda (t)}{t}$$

$$\frac{1}{a}\frac{d^{2}a}{d{t}^2}=\frac{\lambda (t)(\lambda (t)-1)+\mu (t)}{t^2}$$

$$\mu (t)=\frac{d\lambda (t)}{d\ln (t)},$$

$$R\left(\theta ,\varphi \right)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{+\ell }{a}_{\ell }^m{Y}_{\ell }^m\left(\theta ,\varphi \right)$$

$$C_{\ell }=\overline{{\left|a_{\ell }^m\right|}^2}=\frac{1}{2\ell +1}\sum _{m=-\ell }^{+\ell }{\left|a_{\ell }^m\right|}^2,$$

$${\int }_{\theta =0}^{\pi }{\int }_{\varphi =0}^{2\pi }|R\left(\theta ,\varphi \right){|}^2\mathrm{d}\Omega =\sum _{l=0}^{\infty }\sum _{m=-\ell }^{+\ell }{\left|a_{\ell }^m\right|}^2=\sum _{l=0}^{\infty }\left(2\ell +1\right)C_{\ell }$$

$$C=\frac{Q}{\mathrm{amu}A{t}_e},$$

$$f(t)=f(t_0)\exp \left(-\frac{t}{t_e}\right),$$

$$C=\frac{Q}{\mathrm{amu}{A}_P(t_e-t_{e,P})},$$

$$f(t)=f_P(t_0)\left(\exp \left(-\frac{t}{t_e}\right)-\exp \left(-\frac{t}{t_{e,P}}\right)\right),$$

$${}_{28}^{56}\mathrm{Ni}\to {}_{27}^{56}\mathrm{Co}+\gamma +{\nu }_e$$

$${}_{27}^{56}\mathrm{Co}\to \{\begin{array}{ll}{}_{26}^{56}\mathrm{Fe}+\gamma +{\nu }_e& (81\%)\\ {}_{26}^{56}\mathrm{Fe}+\gamma +{\nu }_e+e^{+}& (19\%).\end{array}$$

Full text

supernovae, supernova remnants

††facilities: Simulations were performed on the iTHEMS clusters at RIKEN.††software: HEALPix (Gorski et al., 2005), SciPy (Jones et al., 2001), Matplotlib (Hunter, 2007)

1 Introduction

Supernovae (SNe) are some of the most extreme events in the Universe. As endpoints of stellar evolution (of either massive stars or white dwarfs), they are important to understand the life cycle of stars in galaxies and play a crucial role in the synthesis of chemical elements. Despite substantial effort and recent progress, the very mechanisms by which stars explode have not been fully elucidated. For core-collapse supernovae (type II or Ib/c SNe), some self-consistent simulations have been able to unbind a star, although computational resources are insufficient to assess whether such models produce the required energetics in the end. It appears that multi-dimensionality, while challenging to simulate, is essential: instabilities play a key role in reviving the stalled shock inside the star (for recent reviews, see Janka (2012); Janka et al. (2016)). For thermonuclear supernovae (type Ia SNe), it is still unknown whether they are produced by single-degenerate or double-generate progenitors – or a mixture of both (see Hillebrandt et al. (2013) for a review). Yet, the (empirical) properties of their light curves have become a cornerstone of precision cosmology (Riess et al., 1998; Perlmutter et al., 1999).

The interaction between the stellar ejecta and the ambient medium generates a supernova remnant (SNR). Strong shock waves are produced, which convert the kinetic energy of the explosion into heat, supra-thermal particles, and magneto-hydrodynamic turbulence. The supersonic motion of the ejecta triggers a forward shock (FS) running into the interstellar medium (ISM), while the slowing down of the ejecta triggers a reverse shock (RS) running back into the ejecta. A shell of shocked matter forms, which is an X-ray emitter for typically thousands of years. The SNR is thus one of the ways to study the explosion mechanism, after the SN phase per se. The multidimensional dynamics of SNRs have been studied using numerical simulations for some time (e.g. Chevalier et al., 1992; Blondin & Ellison, 2001). The SNR phase presents a hydrodynamic instability: the Rayleigh-Taylor instability (RTI) grows at the interface between the ejecta and the circumstellar medium, called the contact discontinuity (CD). Ferrand et al. (2010) and Warren & Blondin (2013) studied the evolution of a young SNR taking into account both the RTI at the CD and efficient particle acceleration at the FS. A common trait of all these previous works is the use of idealized initial conditions. Ferrand et al. (2010) relied on the Chevalier (1983) semi-analytical solutions for the shocked profiles in the early self-similar phase, while Warren & Blondin (2013) adopted an exponential ejecta profile following Dwarkadas & Chevalier (1998). Note that such initial conditions were effectively one-dimensional (1D), being functions of radius only, even though the subsequent SNR evolution was three-dimensional (3D). The main point of this new work is to use realistic initial conditions for the SNR phase.

Commonly, researchers who study the evolution of SNRs do not simulate the supernova itself, and, conversely, researchers working on the explosion mechanism typically stop their simulations after several seconds or minutes. Our work fills this gap. Our postulate is that the initial stage of the SNR can reveal information about the SN explosion itself. Since 3D explosions are found to be asymmetric to some extent, we are looking for the imprint of the supernova on the observed morphology of the remnant, in terms of the spatial distribution of the ejecta. Recent studies that bridge SN physics and SNR physics, for the case of core-collapse explosions, are the works by Orlando et al on SN 1987A (2015) and on Cas A (2016). In both simulations, a 1D SN model was mapped to 3D before computing the SNR evolution. In the latter one, asymmetries were added in the ejecta, so as to match the asymmetries of the observed SNR. This demonstrates the relevance of morphological studies to investigate the SN to SNR connection. Our approach is to start with a 3D SN model as obtained from state-of-the art simulations, and compute the hydro evolution from there. In this paper we consider only the thermonuclear case. It is known that, statistically, Type Ia SNRs tend to be more spherical and symmetric than core-collapse SNRs (Lopez et al., 2011). However their ejecta also exhibit complex structures, and we show in this work that this leads to observable features in the remnant structure. For the case of core-collapse SNe, we refer the reader to the preliminary reports by Ellinger et al. (2013) and Gabler et al. (2016).

The structure of our paper is as follows. In Section 2 we summarize the SN simulation used and how we proceed for the SNR simulation, we also explain how we analyze the 3D structure of the evolving SNR. In Section 3 we show the development over the first 500 years of a Tycho-like SNR, and the imprint of the initial SN conditions. In Section 4 we discuss the reliability of our results, and comment on implications regarding observations. In Section 5 we summarize and present the next steps.

2 Method

In this section we describe the numerical simulation, from the SN to the SNR, including tools for the analysis of the morphology.

2.1 SN model

We treat the case of a thermonuclear explosion. For this first work connecting SN and SNR simulations, we concentrate on the single degenerate scenario, that is one accreting white dwarf, and we consider the most popular explosion model, a delayed detonation (Khokhlov, 1991). One of the questions regarding thermonuclear explosions is whether the explosion propagates as a sub-sonic front (deflagration), as a supersonic front (detonation), or as a deflagration that at some point(s) transitions to a detonation (the deflagration-to-detonation model, known as DDT). A set of 2D DDT models recovers the range of luminosities, spectra, and colors of SNe Ia (Blondin et al., 2011), 3D DDT models show reasonable agreement with the spectra of normal SNe Ia although they fall short to recover the width-luminosity relation and colors (Sim et al., 2013). The newer 3D simulations were performed with updated methods (Ciaraldi-Schoolmann et al., 2013). While persisting shortcomings can potentially be attributed to the modeling approaches, it is also possible that they point to a general failure of Chandrasekhar-mass models to explain normal SNe Ia (sub-Chandrasekhar mass explosions seem to be a promising alternative, e.g. Sim et al. (2010); Shen et al. (2018)). In this work we explore whether the simulated ejecta structure of a 3D DDT explosion model translates into observable features in the remnant phase, which allows to further constrain this scenario. We note that a DDT model is also favored by Williams et al. (2017) from observations of Tycho’s SNR.

Seitenzahl et al. (2013) presented a suite of 14 three-dimensional, high-resolution hydrodynamical simulations of a DDT explosion of a Chandrasekhar-mass white dwarf. The ignition is an important yet hard to simulate process, so it was parametrized via different choices of ignition kernel locations. We pick the N100 model, which is the most promising one for reproducing the typical brightness of normal SNe Ia by producing about half a solar mass of 56Ni. In the Seitenzahl et al. (2013) simulations nucleosynthesis was computed at runtime only for the few species that significantly contribute to the energetics. A full nuclear reaction network, including 196 isotopes for 30 elements from Z = 1 to 32, was run in post-processing, using the technique of tracer particles. The results of the post-processing were then remapped from a particle representation back to a grid of resolution $200^{3}$. We use these data cubes. The SN simulation was run until a time of about $100\,\mathrm{s}$. By that time the evolution is observed to be essentially self-similar (Röpke, 2005).111We note that things would be different for the core-collapse case, where it could take a week for the ejecta to settle. The radial extent of the ejecta is $\simeq 2.9\times 10^{11}$ cm, that is $9.4\times 10^{-8}$ pc or $5$ solar radii. The total kinetic energy of the ejecta is $E_{\mathrm{SN}}=1.4\times 10^{51}$ erg, for a total mass $M_{\mathrm{ej}}=1.4M_{\odot}$, implying a maximal speed of about 28,000 km s*-1*.

The 3D structure produced by the SN simulation is substantially different from a 1D averaged profile. The 1D radial profile of the mass density, after averaging over all angles, is shown in Figure 1, compared with the exponential and power-law profiles typically used in the SNR community. One can see that it may be approximated by a 2-step exponential profile. 2D slices of the mass density field are shown in Figure 2.222Interactive 3D models are hosted online on Sketchfab at https://skfb.ly/6pKYW. These show iso-contours, of the mass density, and of the abundance of three of the main tracer species (${}^{12}\mathrm{C}$, ${}^{16}\mathrm{O}$, and ${}^{56}\mathrm{Ni}$). These iso-contours, as well as volume renderings of the full data cubes, are part of an immersive science exhibit developed at the Astrophysical Big Bang Laboratory at RIKEN, that can be experienced using room-scale virtual reality (VR) hardware (Ferrand & Warren, 2018). They reveal the complex morphology of the ejecta. Most of the mass is concentrated in the central part, with a compact nickel core, and irregular outer shells of carbon and oxygen. The small-scale structures in the centre are the result of the deflagration phase, while the sharper outer interfaces were produced by multiple detonation fronts.

2.2 SNR evolution

To simulate the SNR phase, we perform hydro simulations in a way similar to Ferrand et al. (2010) and following papers (2012; 2014). We review here the main features of the simulations, concentrating on those aspects that have been updated or are treated differently. We are using a custom version of the code RAMSES (Teyssier, 2002), an Eulerian code on a Cartesian grid. Even though spherical coordinates may feel like a natural choice to describe a SN or SNR, we are precisely interested here in deviations from spherical symmetry. Our computational domain is however spherical, with boundary conditions applied in the corners of the cube. We further comment on the grid geometry and resolution in Section 4.1. We use a grid of size $256^{3}$, which contains the entire SN data cube obtained after remapping. New compared to our previous SNR papers, we simulate the entire remnant, rather than an octant, since we no longer assume any symmetry.

Initial and boundary conditions.

We use the output from the SN simulations described in the previous section as initial conditions for the SNR simulations. We systematically perform two kinds of simulations, one using the angle-averaged radial profile shown in Figure 1, mapped to the 3D grid in all directions, and one using the actual spatial distribution shown in Figure 2. Thus the initial conditions are effectively 1D (radial only) in the former case (we will refer to this as the “1Di” case), while fully 3D in the latter case (we will refer to this as the “3Di” case). The CD, the interface between the ejecta and the ISM, is unstable and develops Rayleigh-Taylor (RT) fingers even when starting from “clean” profiles (see Figure 3).333On our Cartesian grid, the RTI is naturally seeded by numerical noise stemming from discretization errors. The location of the ejecta in the simulation box is known at any time thanks to a passive scalar $f$, set to 1 inside the ejecta at the start time and advected with the fluid.

To run the SNR simulation we need to specify the ambient medium (in the SN simulation, the explosion happens in a quasi-vacuum). For simplicity it is assumed to be uniform, with hydrogen density $n_{H}=0.1$ cm*-3*, equivalent to a total mass density $\rho_{\mathrm{ISM}}=1.4\,n_{H}m_{p}$ assuming standard ISM composition. This choice is motivated so that the values of the physical quantities are representative of Tycho’s SNR (Williams et al., 2013), which is one of the most promising targets for this study amongst young Galactic SNRs (age $\sim 450$ yr, radius $\sim 3-6$ pc depending on actual distance). The characteristic scales of the problem are controlled by the three quantities $E_{\mathrm{SN}}$, $M_{\mathrm{ej}}$, and $\rho_{\mathrm{ISM}}$ (see e.g. Truelove & McKee (1999)), we define them here in the same way as in Warren & Blondin (2013):

[TABLE]

Since no micro-physics is included, the simulations presented in this paper may be re-scaled to other sets of initial conditions using relations (1), (2), (3).

Time evolution.

Assuming self-similarity of the early evolution, we start the simulations at $t_{0}=1$ day, after simply re-scaling the SN data from 100 s to the SNR start time. Self-similarity implies that radial distance scales as time $t$, while mass density gets diluted as $t^{-3}$. We also performed simulations starting at $t_{0}=1$ year, there are no significant differences visible on the SNR structure (and thus on the maps and spectra introduced in the next section) after a few years of evolution, including the numbers and sizes of the RT fingers. This shows that self-similarity is indeed a safe assumption for the very young SNR. We present the results up to 500 yr of evolution, similar to the age of Tycho’s SNR. We emphasize that our purpose here is to study the SNR phase (from a realistic SN model), not the SN phase per se. In particular we are not concerned with the modeling of the light curve, which requires complex radiative transfer simulations to describe the first days/weeks after the explosion (see Sim et al. (2013) for the N100 model). However we checked the possible effect of heating from radioactive decay, see the discussion in Section 4.2 and more details in Appendix C. Radiative cooling is not included in the simulation, since it is negligible for this phase of the SNR dynamical evolution.

One distinctive technical aspect of the simulations by Ferrand et al. (2010), re-used and extended here, is that we are working in an expanding grid. The number of cells is fixed, but their physical size increases, so that the SNR keeps about the same size relative to the grid. Factoring out the global evolution of the remnant allows us to maintain the same relative resolution in the ejecta and shocked region. From 1 day to 500 years, the SNR extends by a factor of about 50,000. The principle of the comoving transformation with equations is given in Appendix A. It affects all variables, hydro quantities as well as space and time. The key quantity is the scale factor $a(t)$, which describes the physical grid size as a function of time. In previous works, we considered the ideal case where the scale factor can be written as a power-law $a(t)\propto t^{\lambda}$ with constant parameter $\lambda$. New in this work, we simulate the SNR as it smoothly transitions from the free-expansion phase into the ejecta-dominated phase, up to the time when the swept-up CSM becomes important (the Sedov-Taylor time for our model SNR happens to be about 500 yr). We are therefore using a variable expansion law. Then no simple analytical relation exists between comoving time and physical time, we numerically compute the comoving time as a function of physical time given the observed scale factor $a(t)$. See details in Appendix A.

2.3 Analysis of the SNR morphology

The main focus of our work is on the (3D) morphology of the SNR. Young SNRs are best observed in the X-ray domain, where one can image the emission from the hot plasma (Vink, 2012). We therefore focus our analysis on the shocked region, bounded by the reverse and forward shocks. While most SNR studies assume spherical symmetry, we want to know how the shocks develop from a more realistic ejecta distribution. In between the shocks, the ejecta are delimited by the CD, where the RTI generates turbulence. As 3D simulations of SN explosions have become available, a question arises, that to our knowledge has not been addressed yet, of how the RTI of the SNR phase develops on top of the three-dimensional structures resulting from the SN phase.

To investigate the signature of the N100 SN model, we thus choose to analyze over time the morphology of the three wave fronts: forward shock (FS), contact discontinuity (CD), reverse shock (RS). The CD is located using the ejecta fraction, while the shocks (RS and FS) are detected as pressure jumps. Since we are working in a Cartesian grid, we extract a three-dimensional mesh made of quadrilateral cells. Our simple detection scheme ensures that we extract a contiguous mesh for each wave surface. Implementation details are given in Appendix B.

The location of each wave (radial distance from the SN center) is then projected on the surface of a sphere for visualization. For the tessellation of the sphere we are using the efficient HEALPix scheme (Gorski et al., 2005). The maps in the left panels of Figures 4, 5, and 6 show the normalized wave radius $R=(r-\langle r\rangle)/\langle r\rangle$, where $\langle r\rangle$ is the average over all directions at a given time. (It is this average value, for the CD, that defines the global expansion law of the grid.) We use the Mollweide projection of the sphere, which is an equal area projection. We checked that the maps presented here are converged w.r.t. the projection of the 3D meshes.

Finally, to quantify the asymmetries observed, we expand the function $R(\theta,\phi)$ in spherical harmonics. The spectra in the right panels of Figures 4, 5, and 6 decompose the wave surface as a function of the angular scale of the fluctuations of $R$.444We chose to analyze three surfaces of particular significance. In principle one could also analyze the 3D density field as a whole, by performing a spherical Fourier-Bessel decomposition, with spherical Bessel functions $j_{\ell}(kr)$ describing the radial part and spherical harmonics $Y_{\ell}^{m}(\theta,\phi)$ describing the angular part (Leistedt et al., 2012). At angular wavenumber $\ell$, the typical angular scale probed is $\pi/\ell$. The power $C_{\ell}$ plotted is normalized in such a way that each grayed bin is the contribution of wavenumber $\ell$ to the total variance of $R$ over the sphere (see Appendix B for details).

3 Results

We are now equipped to examine the time evolution of the SNR resulting from the N100 SN model.

Density maps

In order to show the inner structure of the remnant, we present slices of the mass density field in Figure 3. A movie from 1 yr to 500 yr is available online, we show here three representative snapshots at 1 yr, 100 yr, and 500 yr.555A volume rendering of the 3D data cube, evolving in time, can also be explored in the VR demo previously mentioned, as explained in Ferrand & Warren (2018). At 1 yr (top snapshot) the ejecta distribution is still very similar to the initial distribution assumed at 1 day. The shocks are present on this plot, although barely visible even on a logarithmic scale. The FS actually appears on the first few time steps, quickly followed by the RS. At 100 yr (middle snapshot) the shell of shocked material, bounded by the RS and the FS, becomes more apparent; at 500 yr (bottom snapshot) it is all that can be seen on the map. From about 10 yr the RTI becomes visible by eye, producing the characteristic pattern of holes and fingers along the CD.

On all the snapshots we compare the two types of initial conditions: effectively 1D on the left (1Di case, depending on $r$) versus fully 3D on the right (3Di case, depending on $r$, $\theta$, $\phi$). Maps on the left are representative of previous works on SNRs (e.g. Blondin & Ellison (2001); Ferrand et al. (2010); Warren & Blondin (2013)). Maps on the right are more realistic, being obtained from an actual SN simulation. And encouragingly, they do look closer to the morphology of a young SNR like Tycho (see Warren et al. (2005) for the X-ray images, more on this in the next section). The 1Di case of spherically symmetric initial conditions is used as a reference. It shows the RTI generated purely from the SNR phase, while the 3Di case shows the RTI growing on top of ejecta with an already complex morphology. This is particularly obvious in the first hundred years of evolution, less so by the final time of 500 yr. One can see how the details of the SN explosion, that is the initial shape of the ejecta, are progressively forgotten over time. Still, the SN imprint is clearly visible at 100 yr, and not fully erased at 500 yr. We also note that, when ejecta have asymmetries from the start, it appears easier for the RT fingers to get into contact with the FS, a phenomenon that has been observed in young SNRs but has proven difficult to reproduce without introducing pre-existing structures (Jun et al., 1996) or increasing the compressibility of the fluid (Blondin & Ellison, 2001). This is in line with the findings of Orlando et al. (2012) on the role of ejecta clumping.

Projected maps and power spectra

The relative variations of the position of the CD, RS, and FS are shown in Figures 4, 5, and 6 respectively, at the same times as in Figure 3 (movies from 1 yr to 500 yr are available online). Maps on the left sides allow for visual inspection of (the 3D surface of) the waves, while spectra on the right sides allow for quantification of the asymmetries. Note that the scales used for the maps and for the spectra are the same on all plots, across the different kinds of waves (CD, RS, FS), initial conditions (1Di vs. 3Di), and times (from 1 yr to 500 yr), to facilitate comparison.

For the shocks, the 1Di case provides a measure of the numerical precision of our code, especially the FS where nothing is expected to happen. We see that the only measurable anomaly is a small $\ell=4$ mode from the cubic geometry of the hydro grid. For the more realistic 3Di case, we see how the shocks first follow the contours of the ejecta, then homogenize in radius thanks to the high pressure in the shocked material. At late times the RS shows slightly rising power at $\ell\geq 10$. At 500 yr the RS is still very close from the CD (see Figure 3), and their surfaces are well correlated: we mostly see the feet of the RT fingers. At 500 yr the FS still shows some protrusions, at scales $\ell<10$. Again comparing with the CD surface, we see that these correspond to the tips of elongated RT fingers, growing along large ridges of ejecta that were present from the beginning. This confirms that the clumpiness of the ejecta enhances their interaction with the FS. We note that, in such a young SNR, the CD as well as the shock fronts probe the outer layers of the ejecta, and thus mostly the effect of the detonation fronts that ended the SN explosion.

For the CD, we observe the same overall regularizing of the initial morphology, plus the RTI growing on top of it. The 1Di case (bottom half of the figures) shows the RTI from the SNR phase, while the 3Di case (top half of the figures) shows the RTI from the SN+SNR phases, comparison between the two allows us to separate the contribution from the SN phase. The SN power is present at small $\ell$ (large angular scales), simply going down in time, while the RTI power is present at larger $\ell$ (small angular scales), going up as well as shifting to smaller $\ell$ in time. This is the expected behavior: the growth rate of the RT instability (in the initial linear regime) is known to scale as the square root of the wavenumber of a perturbation, so it first appears at the smallest scales resolved and grows to larger scales over time (initially by finger growth, then through finger mergers). On the first time frame, we see that the instability had actually started to develop by the end of the SN simulation.

Total power as a function of time

The time evolution of the total power in the fluctuations of the wave positions is shown in Figure 7 (top panel, this corresponds to the gray area on the power spectra). Looking first at the 1Di case (dashed lines), the curve for the FS is a measure of the numerical precision, the curve for the CD is purely due to the RTI, while the curve for the RS shows a faint imprint of the RTI. Looking now at the 3Di case (solid lines), we see the trends discussed with maps and spectra. For the shocks, the power is continuously decaying. For the CD, it also decays until about 300 yr, when it starts to increase. It then follows the same temporal evolution as in the 1Di case, showing that power is now dominated by the RTI of the SNR phase.

As shown in the previous paragraphs, from the comparison of the power spectra between the 1Di and 3Di cases one can separate the ranges of $\ell$ for which the SN phase and the SNR phase dominate. In the bottom panel of Figure 7 we thus plot the integrated power for $\ell<10$ (large scales, thick lines) and $\ell>10$ (small scales, thin lines). Let us concentrate on the evolution of the CD (green lines). We see that the power at small scales in the 1Di case (thin dashed line) is the same as the total power in this case (dashed line of the top plot). And the power at small scales in the 3Di case (thin solid line) follows the same time evolution as in the 1Di case (thin dashed line). This validates our choice of $\ell=10$ to isolate the RTI. Then, the extra power at large scales for the 3Di case shows the contribution of the SN (thick solid line). This component is only decaying over time. Curves for the CD power at $\ell<10$ and at $\ell>10$ cross at $t=178$ yr = $0.4\,t_{\mathrm{ch}}$ (thick vs. thin solid green lines), which is our estimate of the transition time between an SN-dominated morphology and an SNR-dominated morphology for the N100 model.

4 Discussion

Our analysis has shown that the SN explosion can have a measurable impact on the morphology of the subsequent SNR. We now discuss some numerical and physical effects that may impact our results, and their applicability to observations.

4.1 Numerical robustness

We first assess the robustness of our findings with respect to the particular simulation setup we have used.

Effect of numerical resolution.

Our simulations and analysis are performed at a resolution appropriate for the data cubes that we obtained from the SN simulation. To check the effect of spatial resolution, we varied the resolution of the Cartesian grid by a factor of two lower and higher (along each dimension). Using the effectively 1D initial conditions (1Di model), as expected RT fingers appear at increasing $\ell$ as the spatial resolution increases (they will always appear at the smallest scale resolved), however after a few tens of years the power spectrum looks similar, and most importantly the smallest $\ell$ reached is the same as a function of time (meaning that the start time is early enough to properly catch the RTI at these scales); in particular it never reaches below $\ell=11$ (except for the small and known to be artificial $\ell=4$ mode). Using the fully 3D initial conditions (3Di model), differences between the spectra are small when doubling or halving the spatial resolution, and all the previous discussion still holds. So we are confident that the general trends observed on the maps and spectra presented here are revealing the physics of the SN to SNR connection.

We also comment on the use of adaptive mesh refinement (AMR), which is offered by the RAMSES code in the form of octree refinement. To increase the resolution in the most interesting regions, near the shocks and the CD, one can set refinement rules based on gradients of respectively the pressure and the ejecta fraction, as was done for the Ferrand et al. (2010) simulations. This technique was not used in the simulations presented in this paper, since it was found to introduce some noise in the spectra, for a limited extra gain in computing time. However the diagnostics we are discussing, in terms of overall distribution of scales as a function of time, were not significantly affected.

Effect of grid geometry.

Both the SN and the SNR simulations presented here were performed in Cartesian geometry. In such a geometry, slightly different behaviors may be observed when shock waves are propagating along the principal axes versus along a diagonal. This spurious effect is visible at the outer edge of the simulated SNR on Figure 3, for the 1Di case at the final age. However it can be seen from the angular spectra of Figure 6 that the extracted shock surface remains spherical with sufficient numerical precision.

We now compare our results with the work by Warren & Blondin (2013), who did simulations with the code VH1 in spherical geometry (using a Yin-Yang grid). They assumed spherically symmetric ejecta with an exponential profile, which can be compared with our angle-averaged initial conditions (labeled 1Di). They computed the power spectrum from a spherical harmonics expansion, after summing the ejecta density in radial columns (see their Figures 3 and 6). Although the analysis is different, it should provide similar results. They use a dimensionless time, in units of $t_{\mathrm{ch}}$ defined by Eq. (3). On their Figure 3, at times $t=0.15\;t_{\mathrm{ch}}\simeq 67$ yr and $t=0.75\;t_{\mathrm{ch}}\simeq 335$ yr, the power spectrum peaks around respectively $\ell=60$ and $\ell=45$, and drops below 1% of the peak value around respectively $\ell=20$ and $\ell=10$, which in good agreement with our own power spectra for the 1Di case (see Figure 4).

Effect of initial ejecta profile.

As seen from Figure 1 the average radial profile of our SN model is quasi-exponential. For completeness, we also tried to initialize the SNR simulation with the Chevalier (1982) profile. Since this is a 1D model, we compare it to the N100 1Di case. We start the simulations at $t_{0}=1$ yr, since the Chevalier semi-analytical solution is not valid too early on. The width of the shock region is slightly larger at a given time, and accordingly the peak of the power spectrum is at a somewhat larger scale (lower $\ell$), as can be seen in Figure 8. But the minimum $\ell$ reached by the RTI at a given age is similar, so our conclusions are unchanged, regarding the possible contributions of the SN phase versus the SNR phase.

4.2 Possible additional physics

We now consider three energy reservoirs that were not included in the simulations: radioactive decay, the magnetic field, and non-thermal particles.

Effect of heating from radioactive decay.

In this work our focus is not on the SN phase per se, but on its possible impact on the SNR phase. However, one may ask about the possible role of heating from radioactive decay, which powers the optical light curve of the SN. The most important species is ${}^{56}\mathrm{Ni}$, which decays via beta decay to iron in two steps: ${}_{28}^{56}\mathrm{Ni}\rightarrow{}_{27}^{56}\mathrm{Co}\rightarrow{}_{26}^{56}\mathrm{Fe}$. In our SN model the mass of synthesized ${}^{56}\mathrm{Ni}$ is about 40% of the mass of the exploded WD; i.e. $0.6\;M_{\odot}$. The energy released per ${}^{56}\mathrm{Ni}$ nuclear decay is a few MeV (see values in Table 1 in Appendix C), we compute the total energy deposition in our SN model to be $1.1\times 10^{50}$ erg, that is less than 8% of the SN kinetic energy. This may produce significant heating of the ejecta (see the discussion in Woosley et al. (2007)), but cannot alter the overall dynamics of the SNR. And the nickel core is located deep inside the ejecta, it is barely reached by the RS at $t=500$ yr, so we do not expect it can strongly impact the morphology of the shocked region.666The situation may be different for a core-collapse SNR. Compare against the simulations by Wongwathanarat et al. (2015), that show huge pillars of nickel growing toward the CD. To check this, we included heating from radioactive decay in our SNR simulation, under the assumption of local deposition of all the energy released (in photons and in the kinetic energy of positrons), which constitutes an extreme case. See Appendix C for more details. We observed no significant differences on the projected maps and power spectra, so that our above results were not affected. A more precise modeling would require hydrodynamic simulations coupled with radiative transfer, which is beyond the scope of this work.

Effect of the magnetic field.

Although the magnetic field is not expected to be strong enough to alter the development of young SNRs, and the key driver of the RTI is the deceleration of the ejecta, the RTI growth rate may be affected by the magnetization of the ejecta at the CD. From linear perturbation theory, a magnetic field tangential to the interface slows the growth of the instability due to the tension of the field. The instability is actually suppressed for wavelengths smaller than some critical value that is proportional to the energy density of the field (Jun et al., 1995). If the RTI from the SNR phase is suppressed, presumably it would be easier to see the contribution from the SN phase. The development of the RTI into the non-linear stage in a strong magnetic field has been studied by Stone & Gardiner (2007) in 3D, and including particle acceleration by Schure et al. (2009) in 2D. We defer the use of MHD simulations to future work.

Effect of particle acceleration.

Particle acceleration is expected to happen at the blast wave of SNRs via the mechanism of diffusive shock acceleration (DSA), and if particle acceleration is as efficient as it is believed to be, so that SNRs are significant producers of Galactic cosmic rays, then the fraction of the SN energy channeled into energetic particles may reach 10% or more, which will impact the dynamics of the SNR (e.g. Ferrand et al. (2010) and references therein). The presence of energetic particles increases the compressibility of the fluid (lower effective adiabatic index), so the general effect is to make the shocked region more compact, meaning both narrower and denser. This effect is measurable and probably observed in Tycho’s SNR (Warren et al., 2005), although it is not a huge modification.

We believe that particle acceleration alone cannot strongly affect our conclusions regarding the morphology of the CD. The simulated X-ray maps produced by Ferrand et al. (2012) do not show obvious differences in the size of the structures present inside the SNR (the effect of particle acceleration is more visible in a differential comparison of the emission at different energy bands, that probe different levels of back-reaction). In Warren & Blondin (2013) the limited effect of lowering the adiabatic index on the power spectrum of density fluctuations can be seen in their Figure 6, at their final simulation time $t=2\;t_{\mathrm{ch}}=895$ yr that is almost twice our final time. They compare the standard case $\gamma=5/3$ with the extreme cases of $\gamma=4/3$ and $\gamma=6/5$, which model the effect of the acceleration and escape of the relativistic particles (these simulations were made under the simplifying assumptions of reducing the adiabatic index everywhere in space and at all times).

Regarding the shocks, a reduced distance between the CD and the FS will make it easier for the RT fingers to reach the FS, so that the morphology of the FS may contain some echo of the CD structure. Acceleration at the shock will not impact the shock surface by itself, unless the injection efficiency varies significantly on small scales all along the shock surface (observations of Tycho’s SNR do not show evidence for this). The FS is more likely to be affected by other, external factors. As for the RS, efficient acceleration is not expected to occur in general, see the discussion in Ellison et al. (2005) (and there is no evidence for it in Tycho’s SNR).

4.3 Implications for observations

Finally we compare our results with observations, especially of the young SNR Tycho (G120.1+1.4 = SN 1572, $447$ yr), which is the closest to our simulations.777 This SNR seems to be the best candidate to apply our findings, given its age and globally simple appearance. Other young type Ia Galactic SNRs include G1.9+0.3, the youngest one known ($\simeq 150$ yr, Borkowski et al. (2013)); Kepler (G4.5+6.8 = SN 1604, $415$ yr), with a more complicated morphology presumably from the progenitor’s mass loss history (Reynolds et al., 2007); and SN 1006 (G327.6+14.6, $1013$ yr), with a characteristic bilateral morphology from the orientation of the magnetic field (Winkler et al., 2014)). Other young SN Ia SNRs in the LMC are SNR 0519-69.0 ($\simeq 450$ yr, Kosenko et al. (2010)) and SNR 0509-67.5 ($<1000$ yr, Warren & Hughes (2004)). The most relevant energy band for this study is the X-ray band, which reveals the shocked matter, and for which we have instruments with sufficient angular resolution to study the detailed morphology. For Tycho’s SNR, we refer the reader to the Chandra observations presented in Warren et al. (2005) and in Williams et al. (2017).

Simulations of this SNR by Ferrand et al. (2012) and Warren & Blondin (2013), that were made from spherically symmetric initial conditions, have reproduced the fleecy aspect of the inner ejecta, but have both failed to produce structures on scales as large as seen on the X-ray maps. Warren & Blondin (2013) noted that their power spectrum is significantly different from observations by Warren et al. (2005), who found that the peak for the angular fluctuations of the CD (as seen in projection) occurs for $\ell=6$. Our 1Di simulations confirm this: it appears that simulations made with different codes and both assuming spherical ejecta are in agreement with each other, and are in disagreement with X-ray observations of Tycho’s SNR. And by lifting the assumption of spherically symmetry, our new work provides an explanation for this difference: the power observed at low $\ell$ does not come from the RT growth but from the SN explosion. More precisely, we have observed in Section 3 that the contribution to the morphology of the CD from the SN phase is visible at low angular wavenumber $\ell<10$, and decreasing over time, while the contribution from the SNR phase is visible at larger $\ell$, increasing over time as well as shifting to lower $\ell$. Interestingly, around 500 yr (about the age of Tycho’s SNR) it so happens that the two contributions nicely match on the spectrum plot (last panel of Figure 4): the power spectrum looks like a single distribution, just extending over a wider range of $\ell$ than would be expected from RTI alone. This could easily be confused with an enhanced RTI growth. And this actually looks like what is being observed in Tycho’s SNR.

In Figure 9 we show an approximation of the thermal emission from the shock-heated ejecta, which scales as the square of the density and is observed summed up in projection along the line of sight. The more realistic SN model (used on the right) does produce a more realistic-looking mock SNR map, with a limb that shows more pronounced gaps and peaks like on Tycho’s X-ray image. The emission from the shocked plasma actually depends on the electronic density, the electronic temperature, and the ionization states of all the species present, which for young SNRs depend on the history of the shocked material since the crossing of the shock (see Ferrand et al. (2012) and references therein). In a forthcoming paper we will calculate X-ray maps and present a dedicated analysis of Tycho’s SNR data.

For the case of Tycho’s SNR, Williams et al. (2017) compared 3D simulations made from smooth initial ejecta profiles and from clumpy initial ejecta profiles (with a characteristic scale set by hand). Measuring the deceleration parameter in ejecta knots, they were not able to tell apart the two models using their X-ray observations, although visual analysis of the 2D projected images shows clear differences in the morphology of the SNR depending on the initial conditions (compare their Figure 10 to our Figure 9). Sato et al. (2019) use the genus statistics, a topological method, to quantify these differences. They find that the observed morphology of Tycho better matches the case of clumpy initial ejecta. This comforts our view that the SN explosion has a visible imprint on this SNR.

In this work we deliberately assumed a uniform ISM, so as to reveal the contribution from the SN. The power in the fluctuations of the RS and FS is continuously decaying over time, and after hundreds of years is very small, probably too difficult to detect. From observations of young SNRs it is relatively easy to delineate the FS, which is traced by the synchrotron emission of accelerated electrons (see again Warren et al. (2005) as well as Cassam-Chenaï et al. (2007) for the X-ray observations of Tycho’s SNR, and Ferrand et al. (2014) for the corresponding simulations). But by that age and on those scales the FS will most probably be impacted by inhomogeneities in the ambient medium, which are not included in our simulations.888For an example applied to Tycho’s SNR, see the recent work by Fang et al. (2018), who considered the impact of a stellar outflow from the progenitor system to explain irregularities along the outer shape of the remnant. The RS, that runs inside the ejecta, will not be as immediately affected as the FS by the ISM. However it is in general more difficult to locate from observations, and as a result comparisons will be subject to some degree of interpretation (e.g., for the case of Tycho, relying on a particular emission line, the Fe K$\alpha$ line, as suggested by Warren et al. (2005)). So we expect it will be more difficult to draw firm conclusions from studies of the SNR shocks. Density fluctuations in the shocked ejecta (as delineated by the CD) seem to contain more exploitable information. We note that, as time goes by, the RS will reveal in X-rays the ejecta located deeper and deeper inside the ejecta. At the current age of Tycho, the core of the ejecta is invisible in X-rays, and unfortunately the cold ejecta have not been detected in a type Ia SNR at other wavelengths, be it in radio (molecular lines), infrared (dust) or gamma (radioactive lines). Although such detections are challenging, they would open a new window to the SN and SNR physics, that we could take advantage of using our simulations.

5 Conclusion

In this paper we have studied the imprint of a SN explosion on the morphology of the resulting SNR, by the means of 3D hydrodynamics simulations. We have concentrated on the thermonuclear case (for type Ia SNe), and used the recent N100 model from Seitenzahl et al. (2013), a delayed detonation (DDT) of a Chandrasekhar mass white dwarf, that we put in a low-density uniform medium thought to be a good first approximation for Tycho’s SNR. Interestingly, using a realistic 3D SN model leads to larger scale and more irregular structures, which were not seen in SNR simulations made from (semi-)analytical SN models, and which better match X-ray observations of Tycho’s SNR. To quantify the effect, we have extracted the time-dependent surfaces of the forward shock, contact discontinuity, and reverse shock, and analyzed their angular scale structure (see Figures 4, 5, 6). By comparing, as initial conditions, the actual three-dimensional distribution of the ejecta obtained from the SN simulation with an angularly smoothed version of it, we have shown that morphological signatures of the explosion can be seen clearly in the first hundred years, and can still be detected after a few hundred years. For the forward and reverse shocks the angular power from the SN only decreases over time, and for the forward shock may be further erased by inhomogeneities in the ambient medium. For the CD the SN power decreases over time while the SNR power increases, due to the Rayleigh-Taylor instability. Thanks to the expansion in spherical harmonics we have devised a way to separate the contributions of the SN and SNR phases. For the CD they are dominant respectively above and below scales of about $20\deg$; in terms of total power the SNR phase takes over at around 200 yr (see Figure 7). For our Tycho-like SNR the two contributions to the angular spectrum, from the SN phase and from the SNR phase, are found to have similar power at the current observed age, which mimics a faster RTI growth (see power spectrum in Figure 4 (c)). We suggest this is what is actually observed in Tycho’s SNR.

The obvious extension of this work is to re-do our simulation and analysis with different SN models. We have access to a grid of DDT models, as well as a grid of pure deflagration models. We can also test sub-Chandrasekhar mass Type Ia supernova models, as well as other channels for Type Ia SNe (double degenerate scenario). By investigating all these different models, we will be able to assess the range of behaviors that may be observed for the case of type Ia SN(R)s, in terms of the surviving asymmetry of the ejecta. We note that the N100 model used here is a 3D explosion model having a quite high degree of spherical symmetry. Other 3D DDT models, and especially merger models, can exhibit much greater asymmetry. We also note that in the single-degenerate case there is a surviving companion, which would also interact with the ejecta (e.g. Liu et al., 2012; Pakmor et al., 2008), and could perhaps lead to additional structure in the remnant. With this program, we can determine if and how the SN signatures that we have established in this paper can be used to discriminate between different thermonuclear SN theories, which is an important step to be able to use SNRs as probes of the explosion physics.

This work was supported in part by RIKEN Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS). This work was supported in part by a RIKEN pioneering project “Extreme precisions to Explore fundamental physics with Exotic particles (E3-Project)”. The work of FKR was supported by the Klaus Tschira Foundation and by the German Research Foundation (DFG) via the Collaborative Research Center SFB 881 “The Milky Way System”. IRS was supported by Australian Research Council Grant FT160100028. We thank the anonymous referee for comments that helped clarify the text.

Appendix A Co-expansion transformation

We detail here the transformation that makes the computational grid co-expand with the SNR, without any need to add grid cells or refinement layers. The technique is commonly used for cosmological simulations, it was applied to SNRs by Fraschetti et al. (2010), following Poludnenko & Khokhlov (2007).999The supernova explosion itself was computed on an expanding grid, using a different “moving grid” technique (Röpke, 2005). Everywhere in this section tilde quantities refer to transformed quantities, i.e. after the change of variable.

A.1 The general transformation

We need to define the transformation for the independent space, time, and mass dimensions. Distance $r$, time $t$, and mass density $\rho$ are transformed as

[TABLE]

where $\alpha$ and $\beta$ are free parameters, set below. This in turns implies the velocity transformation

[TABLE]

with expansion factor (“Hubble factor”)

[TABLE]

and also

[TABLE]

for pressure $P$, total energy $E$, energy per unit mass $\epsilon$.

The Euler equations (here written in conservative form) become

[TABLE]

[TABLE]

[TABLE]

where $\nu$ is the number of spatial dimensions (for our case $\nu=3$), $\tilde{H}$ is the comoving expression of the expansion factor:

[TABLE]

and we have introduced the (comoving) inertial force

[TABLE]

where

[TABLE]

In order to keep the Euler equations as close as possible from their original conservative form, we choose $\alpha=\nu=3$ and $\beta=1$. Then, for a fluid with adiabatic index $\gamma=5/3$ (no internal degrees of freedom), the only extra terms are the ones involving the force $\mathbf{\tilde{F}}$. This stems from the fact that we are now working in a non-inertial frame.

A.2 Application to an expanding SNR

The expansion law of a SNR is commonly written as a power-law

[TABLE]

where $a_{\star}$ is the size at some reference time $t_{\star}$ (for us the physical start time $t_{0}$ of the SNR simulation), for some constant value of index $\lambda$. Self-similar phases obey such a relation: in the initial free expansion phase $\lambda=1$, the Chevalier solution for the ejecta-dominated phase with power-law profiles of the ejecta (index $n$) and of the ISM (index $s$) gives $\lambda=(n-3)/(n-s)$, in the later Sedov phase conservation of energy implies $\lambda=0.4$. In that case we have simple expressions for the derivatives

[TABLE]

and we can integrate the relation between physical and comoving time:

[TABLE]

where $\tilde{t}_{0}$ is the arbitrary starting value of the comoving time and $\tilde{t}_{\star}=t_{\star}/a^{2}_{\star}$.

When the scale factor has more complicated time dependence,101010For an analytic description of the continuous evolution in time of a (non-radiative) SNR, see Tang & Chevalier (2017). we can always define the instantaneous index

[TABLE]

and we still have

[TABLE]

and in comoving coordinates relations (A12) and (A14). We now write

[TABLE]

with

[TABLE]

which generalizes Eq. (A17).

In practice the simulation is done in two passes. In the first pass we use a fixed $\lambda\lesssim 1$ and an oversized box, and record the average position $r_{\mathrm{CD}}$ of the contact discontinuity over time. We define $a(t)=r_{\mathrm{CD}}(t)/r_{\mathrm{CD}}(t_{0})$, and integrate the physical time $t$ to compute the comoving time $\tilde{t}$ according to Eq. (A2) (from $\tilde{t}=0$ when $t=t_{0}$). We also compute the derivatives of $a(t)$ w.r.t. $\ln(t)$: $\lambda(t)$ and $\mu(t)$, needed to compute $\tilde{H}$ and $\mathbf{\tilde{F}}$, from a spline interpolator of $a(t)$ (we used scipy.interpolate.spalde). In the second pass we use the tabulated scaling law to define the transformed quantities at each time-step.

Finally we note that a side-effect of having a quasi-stationary SNR in the computational grid is the development of a numerical instability at the FS front (carbuncle instability, e.g. Quirk (1994)). Compared to our previous studies, this is far less of a concern when using a 3D SN model, since the FS is never perfectly spherical. Besides, the high-spatial frequency noise generated (at the cell scale) is irrelevant to the effects discussed here.

Appendix B Analysis of the wave structure

In this section we detail the practical aspects, that are specific to our simulation setup, of the extraction and analysis of surfaces inside the SNR.

B.1 Extracting and recording wave fronts

During runtime, at every hydro time step, we look for the position of the waves inside the simulation box. These interfaces are located by finding where some quantity crosses a threshold, along each of the three spatial dimensions. For the CD, this is the ejecta fraction $f$, we set the threshold to $f=0.5$. Varying this threshold between 0.1 and 0.9 causes only minor changes on the maps and spectra. For shocks, the most obvious marker is the pressure jump, the threshold on $\log(P)$ is set to be the arithmetic mean of the lowest (ambient) and highest (shocked) pressure in the box. Numerically the width of the shocks is only of a few cells, so the exact value of the threshold is not critically important, as long as it is applied consistently in time and across simulation runs. The distinction between the FS and the RS can readily be made knowing the ejecta fraction. Working on a Cartesian grid, we thereby collect an (unordered) set of faces that are boundaries between grid cells. We also record their orientation, obtained from the gradient of the quantity $f$ or $P$, for visualization purpose. This scheme ensures that we extract a contiguous surface for each wave. We note that the surface of the CD may present folds, meaning that its radius is a multi-valued function along a single ray from the centre.

The location of each wave (radial distance from the SN center to the face center) is then projected on the surface of a sphere, which is tesselated using the HEALPix scheme (Gorski et al., 2005). We adopt the RING ordering. The HEALPix tessellation is controlled by a resolution parameter $N_{\mathrm{side}}=2^{k}$ where $k$ is an integer, the total number of pixels being $N_{\mathrm{pixel}}=12N_{\mathrm{side}}^{2}=12\times 2^{2k}$. The resolution of the Cartesian grid on the other hand is controlled by the maximum level number $k_{\mathrm{max}}$ of the RAMSES octree, so that the number of cells on the finest level is $2^{k_{\mathrm{max}}}$ along each of the 3 dimensions. We set $k=k_{\mathrm{max}}-1$ so that the angular resolution of the HEALPix map be appropriate for the Cartesian grid (of course the individual Cartesian faces will project at a variety of angles). We tried varying the parameter $k$ for a given $k_{\mathrm{max}}$, and checked that differences between maps of varying HEALPix resolution are modest, and would not affect our diagnostics. Our baseline is $k_{\mathrm{max}}=8$ so that we have $256^{3}=16,777,216$ cells (and $6\times 256^{2}=393,216$ outer faces), and $k=7$ so that for the HEALPix maps $N_{\mathrm{side}}=128$ and $N_{\mathrm{pixel}}=196,608$.

To project the faces that make up the wave surfaces, we break them down recursively into four smaller faces, and assign their radius to the HEALPix pixels in the direction of all the leaf sub-faces. Four recursion steps are sufficient to obtain a converged projected map. This high-resolution projection is done only when outputting the wave data to disk (which is done independently of the grid outputs). Furthermore, for the CD, which might have a complicated morphology in 3D, we tried two projection schemes: collecting all the faces, thereby averaging over folds, or collecting only the outermost faces in a given direction. The two schemes result in small differences that are visible on the maps, but barely affect the spectra defined below. We present here maps made using the first scheme.

B.2 Quantifying asymmetries with spherical harmonics expansion.

At this stage, we have a representation of the morphology of each of the wave fronts (RS, CD, FS) as a function on the sphere, $r(\theta,\phi)$, the radius at which this wave front is found along direction $(\theta,\phi)$. We normalize this function as $R=(r-\langle r\rangle)/\langle r\rangle$ where $\langle r\rangle$ is the average value over all angles, so that we are working with relative fluctuations in radius. To quantify these fluctuations, we expand $R(\theta,\phi)$ in spherical harmonics:

[TABLE]

where $Y_{\ell}^{m}$ is the Laplace spherical harmonic of degree $\ell\geq 0$ and order $m\in\left[-\ell,+\ell\right]$, and the expansion coefficients $a_{\ell}^{m}$ are found by projection of $R$ on the orthogonal basis of the $Y_{\ell}^{m}$ functions. This is done in post-processing on the HEALPix map, using the Python package healpy.111111https://healpy.readthedocs.io The degree $\ell$ determines the spatial frequency of the function over the sphere, which can be estimated as $\pi/\ell$.

The angular power spectrum is defined as

[TABLE]

and Parseval’s theorem reads

[TABLE]

so that the contribution of scale $\ell$ to the total variance of $R$ over the sphere is $\left(2\ell+1\right)C_{\ell}/4\pi$. In all the figures we show the spectrum as a function of $\log\ell$ and therefore plot this latter quantity multiplied by $\ell$, so that the grayed area be the variance of the function $R$.121212We note that the convention used in all CMB plots is neither of these, but $\ell(\ell+1)C_{\ell}/2\pi$, for historical reasons specific to cosmology (the Sachs-Wolfe effect).

Appendix C Energy deposition from radioactive decay

In this section we outline the basics of energy deposition from radioactive decay. The general formalism is taken from Jeffery (1999), see also Nadyozhin (1994).

The energy deposition per unit volume and time is noted $\rho\,\epsilon$ with a time-dependent energy production rate $\epsilon=C\,f(t)$, where $C$ is the energy generation rate coefficient and $f(t)$ is the fraction of species left at time $t$. For primary species:

[TABLE]

[TABLE]

and for secondary species ($P$ designates the parent species):

[TABLE]

[TABLE]

where $Q$ is the mean energy deposit per decay in photons or e+/e- kinetic energy (the energy that can be absorbed by matter) and $t_{e}$ is the $e$-folding time ($t_{e}=t_{1/2}/\ln(2)$ where $t_{1/2}$ is the half-life).

We are here interested in the (completely dominant) decay chain

[TABLE]

So in the hydro code we add a source term for the internal energy $\rho\left(\epsilon_{\mathrm{Ni}}+\epsilon_{\mathrm{Co}}\right)$, where $\rho$ is the ejecta mass density and the initial fraction of ${}^{56}\mathrm{Ni}$ is known from the SN simulation. The relevant numerical values are compiled in Table 1. We indicate the total energy deposit as well as the part that goes into photons and kinetic energy of positrons, that may be converted into heat (excluding neutrinos, which nearly entirely escape).

In our simple tests, we assume local energy deposition, without radiative transfer. Initially the optical depth is very large, caused by Compton scattering. The transition time between optically thick and optically thin periods is estimated by Jeffery (1999) to be about 40 yr for their fiducial Type Ia SN (their Equation (28)). For comparison, 97% of the energy liberated from radioactive decay of ${}^{56}\mathrm{Ni}$ is deposited within the first year after the explosion. The total energy deposition for the N100 model is $1.1\times 10^{50}$ erg, that is less than 8% of the SN kinetic energy.