Introduction

Supernova remnants (SNRs) are formed from stellar ejecta driven out by supernova explosions. High-resolution observations have made it possible to study the three-dimensional structure in many SNRs, e.g., Tycho's SNR (Williams et al. 2017), N132D (Vogt & Dopita 2011;Law et al. 2020), Cassiopeia A (DeLaney et al. 2010;Milisavljevic & Fesen 2013, 2015), and SN 1987A (Larsson et al. 2016).

Hydrodynamic (Chevalier & Klein 1978;Velazquez et al. 1998) and magnetohydrodynamic (Jun & Norman 1995;Bucciantini et al. 2004) Rayleigh-Taylor Instabilities (RTIs) are expected to be possible factors responsible for these anisotropies. The interaction of the ejecta with the circumstellar material (CSM) drives a forward shock and a reverse shock into the CSM and the ejecta, respectively. The contact discontinuity, or the interface between the shocked ejecta and the shocked CSM, has been shown to be Rayleigh-Taylor unstable (Chevalier et al. 1992). However, the influence of seed anisotropies on the RTI structures is unclear. Large-scale anisotropies have often been attributed to those present in the SN explosion itself, e.g., Tycho's SNR (Ferrand et al. 2019) and SN 1987A (Orlando et al. 2020). For instance, the standing accretion shock instability (SASI; Blondin 2005) is thought to be one of the primary drivers of large-scale anisotropies in the ejecta structure of core-collapse supernovae (SNe) during the first seconds after the collapse (Blondin & Mezzacappa 2007;Iwakami et al. 2008). Other important drivers of anisotropies endemic to the SN include jets (González-Casanova et al. 2014) and stochastic processes occurring shortly after core collapse due to events such as radioactive decay by 56 Ni plumes and neutrino-driven convection (Wongwathanarat et al. 2015(Wongwathanarat et al. , 2017;;Orlando et al. 2016Orlando et al. , 2021;;Gabler et al. 2021). In some other cases, it has been argued that dense, asymmetric circumstellar environments are required to explain atypical features in some young SNRs (Celli et al. 2019;Sano et al. 2020aSano et al. , 2020b)).

Numerical studies of SNRs have identified power spectral analysis as a powerful tool for quantifying anisotropies in the remnant models, (e.g., Warren & Blondin 2013;Ferrand et al. 2019Ferrand et al. , 2021;;Gabler et al. 2021;Polin et al. 2022). In particular, Polin et al. (2022) use this technique to identify a dominant angular scale θ 0 in RTI structures formed in SNRs, which separates features formed by turbulence from those endemic to the SN. They find that θ 0 is a function of the density scale height of the ejecta. They also find that the dominant angular scale or the shape of the angular power spectrum of the SNR structure is not affected by fluctuations in the surrounding media in general. However, their value of θ 0 and the shape of the turbulent power spectrum at large wavenumbers do not agree with corresponding results found by Warren & Blondin (2013). In addition, Polin et al. (2022) obtain their results using 2D axisymmetric numerical calculations, which are in general different from 3D calculations, and do not necessarily give results that agree with them. For example, laboratory experiments and 3D numerical calculations generally exhibit a turbulent cascade of energy from large length scales to smaller ones, a feature that is absent in 2D calculations of turbulence. Hence, we perform a follow-up investigation of hydrodynamic instabilities in SNRs in this work using a suite of 3D numerical models. These models were developed using a new expanding-mesh hydro code called Sprout. In general, this is a complex problem to study because of the multitude of phenomena involved. SN explosions are often highly asymmetric, as seen for Cas A (Rest et al. 2011), W49B (Lopez et al. 2013;Zhou & Vink 2018), RCW 86 (Broersen et al. 2014), SN 1987A (Wang et al. 2002), and many other SNRs. Moreover, CSM environments may develop largescale asymmetries and complex density profiles as a result of the mass loss and evolutionary history of the progenitor system and have to be constrained using theoretical models (Dragulin & Hoeflich 2016) and multiwavelength observations (Chomiuk et al. 2012;Cendes et al. 2020;Sand et al. 2021). For instance, there is a growing subclass of SNe that are expected to have run into a circumstellar shell much denser than the ambient medium (Smith et al. 2008;Jencson et al. 2016;Dickinson et al. 2023), which are formed as a result of a brief period of very high mass loss by the progenitor. Our models do not attempt to simulate all these scenarios explicitly. Instead, we perform a systematic analysis on a range of idealized models to understand the physical laws governing hydrodynamic instability resulting from the expansion of SNR ejecta, and their dependence on system parameters such as the progenitor density structure. We also incorporate seed anisotropies or "clumps" in the CSM and in the SN ejecta for some of our models via a parameterized model. This allows us to study the effects of these clumps on the turbulent structures in terms of their size and mass. Our models, on account of being mostly spherical, are more directly comparable to remnants of Type Ia SNe since most of them are expected to result from nearly spherical explosions (Soker 2019, and references therein). Nevertheless, the conclusions we obtain about the nature of turbulence in SNRs are expected to be general and hold for core-collapse SNRs as well.

This paper is organized as follows. In Section 2, we describe the numerical method employed by Sprout briefly along with the model setup. Section 3 describes the power spectral analysis technique used to study the models. The results are presented in Section 4. Sections 5 and 6 describe how the SNR structure is affected by the presence of anisotropies in the ejecta and the CSM, respectively. In Section 7, we summarize our results and discuss our expectations for structure of observed SNRs.

Numerical Setup

Numerical Method

Sprout (Mandal & Duffell 2023) solves Euler's equations of ideal fluid dynamics in conservative form on a Cartesian mesh using a second-order upwind Godunov method. Secondorder spatial accuracy is obtained through piecewise linear reconstruction of the primitive variables, using a minmod slope limiter. Sprout utilizes the moving mesh methodology (Springel 2010;Duffell & MacFadyen 2011) to expand its mesh with time, preserving all aspect ratios. The mesh motion thus results in all position vectors r being updated to r + δr after a time step Δt as follows:

Thus we can resolve fluid motion for several orders of magnitude in distance, although it is only useful to resolve one or two orders of magnitude at any given instant for a reasonable duration of computation time. This also gives us the advantage that each order of magnitude in size should take roughly the same amount of time to solve for, provided the fluid velocity remains constant.

Initial Conditions

The outer ejecta in young SNRs are expected to be cold and expanding homologously into a stationary CSM, with a density profile declining steeply with radius. We choose a power-law density profile (Arnett 1988), in accordance with calculations of adiabatic explosion of polytropes (Matzner & McKee 1999;Chomiuk et al. 2016), although there are other models for SN ejecta; for example, Warren & Blondin (2013) employed an exponential ejecta profile consistent with the popular "W7" model (Nomoto et al. 1984) in their initial conditions. Even more realistic initial conditions were employed by Ferrand et al. (2019Ferrand et al. ( , 2021) ) and Orlando et al. (2021), obtained from 3D models of SN explosion. We favor a simple analytical model even though it is not ideal for directly comparing to data, because it allows us to systematically vary the relevant parameters and obtain predictions for physical laws governing their impact on turbulent phenomena in SNRs.

One advantage of using the power-law density profile is that it provides a self-similar solution for the SNR (Chevalier 1982a), allowing the turbulence to reach a statistical steady state (as will be shown in Section 4.2). After a steady state has been attained, the SNR model can be studied at any single instant of time to completely characterize the turbulence, which is applicable to all times. Our initial conditions are then given by

r P Ejecta: , CSM: , 0 10 bothcases . 2 n n s 3 6

We use the above initial conditions to obtain 3D numerical self-similar solutions for supernova remnants. The value of n is expected to be between 7 and 12 (Colgate & McKee 1969;Chevalier 1982b), and s is chosen to be either 0 (for a CSM of constant density) or 2 (for a stellar-wind environment). We vary n between 9 and 11 for both values of s. We choose an adiabatic index γ = 5/3 for all our calculations. The calculations start at some time t i (in code units), and are continued for four decades, that is, until t = 10 4 t i . Our calculations assume a fiducial value of t i = 1 day in physical units. Thus, the calculations end at a physical time corresponding to about 30 years. As will be shown in Section 4.2, we find that the turbulence reaches a steady state by this time. The solution should remain statistically self-similar after this point, until the reverse shock reaches the core and the mass swept up by the shock becomes comparable to the mass of the ejecta. This is when the SNR transitions from the ejecta-dominated phase to the Sedov-Taylor phase, which typically happens around 10 4 years after the SN explosion (Woltjer 1972).

The expansion of the computational domain ensures that the position of the contact discontinuity is roughly constant with respect to the domain. This requires

where H(t) is as described by Equation (1) and R c (t) ∝ t ( n-3)/( n-s) is the speed of the contact discontinuity.

Our runs have a fiducial resolution of 512 3 . Most of our computations model only one octant of the spherical SNR to save computation time. For one of the cases (n = 9 and s = 0), we model the full spherical remnant at the same resolution to ensure we get the same results.

Power Spectral Analysis

The anisotropies seen in the shocked region in our solutions are best quantified by defining some quantity on the (nearly) spherical ejecta surface. The relative strength in the anisotropies of all angular scales can be obtained by expanding this quantity (which we now call f (θ, f)) as a linear combination of spherical harmonics:

where a lm , or the amplitude of a given harmonic mode, can be calculated using

There are many possible choices for f (θ, f), including the distributions of contact discontinuity and forward shock radii of the remnant (Ferrand et al. 2019(Ferrand et al. , 2021)). We follow the choice made by Polin et al. (2022) and obtain spherical surface maps for radial velocity and density integrated along some radial lines of sight originating from the center of the supernovae. They are defined as follows:

The pressure weighting seeks out the relevant quantities in the shocked region since the fluid is cold everywhere else. The surface maps are then normalized through division by the angle-averaged value (of velocity or density). Given the amplitudes of the harmonics, it is then possible to calculate a power spectrum as a function of the harmonic number l:

We obtain power spectra for our normalized surface maps using the SHTOOLS package (Wieczorek & Meschede 2018).

Results

Figure 1 shows velocity and density isosurfaces from our 3D solution for n = 9 and s = 0. The outermost isosurfaces were chosen to have a value such that they lie around the contact discontinuity. Thus they provide a qualitative view of the typical size of the Rayleigh-Taylor fingers. In Figure 2, we show the radially integrated velocity and density maps (as defined in Equation ( 6)) obtained for this 3D solution. Figure 3 shows the resulting power spectra (for the n = 9, s = 0 solution). The left panel shows the velocity spectra obtained at different computational resolutions, demonstrating convergence of the spectra. The spectra are seen to peak at some intermediate spherical harmonic (hereafter l 0 ). There is another peak at l = 4. This is a numerical artifact arising from a very mild carbuncle instability in our solutions, which is also visible in the Mollweide projections in Figure 2 as small squareshaped structures. They form due to the shock not being aligned to the grid (Quirk 1994;Robinet et al. 2000). We . As these represent idealized self-similar solutions, these results apply to any time after the solution has converged to a statistically self-similar state. For a typical Type Ia SN, they apply to any time after about 30 yr and before 10 4 yr, which is when an SNR transitions to the Sedov-Taylor phase (Woltjer 1972). The values have been chosen so as to focus on the shocked region between the ejecta and the CSM, which houses the turbulent RTI structures. The visualization has been generated using the open-source software VisIt (Childs et al. 2012).

smoothen the spectra by convolving the raw data with a Gaussian kernel (shown in black in the left panel). This generates a peak at l = 4 with a width approximately equal to that of the smoothing kernel. Since this is a purely numerical feature, we develop another smoothened version of the spectra (shown in green) that ignores the bump at l = 4. This is useful for identifying features visually and will be used for subsequent plots, although all analysis was done using the raw data. The right panel shows the smoothened velocity and density power spectra. We find that both power spectra display a broken power law, similar to previous studies (Warren & Blondin 2013;Polin et al. 2022). For l smaller than l 0 , or larger angular scales, the power spectra can be approximated by a power law with a positive slope. For smaller angular scales (l > l 0 ), the spectra behave like a power law with a negative slope.

General Features of the Power Spectra

The velocity and density power spectra for all the different ejecta and CSM density profiles implemented in our models are shown in Figure 4. We find that they are fairly consistent with those from previous studies. The break wavenumber, l 0 , is found to be proportional to the ejecta density power-law index, n, the same as found by Polin et al. (2022) and Warren & Blondin (2013). The exact value is found to be

which is close to ∼3n as found by Polin et al. (2022). This harmonic corresponds to a dominant angular scale θ 0 (l 0 ∼ π/θ 0 ). Our result is consistent with their interpretation that the density scale height of the ejecta,

is roughly of the order of the arc-length s enclosed by the dominant angular scale,

Warren & Blondin (2013) also obtain a dominant angular scale in their 3D SNR models with the exponential density profile for the outer ejecta (r µ - e v v e , v = r/t, and v e is the ) (as defined by Equation ( 6)) for the n = 9, s = 0 Chevalier solution. The repeating square-like features (more prominent in the velocity map) are formed by the presence of a carbuncle instability in our solutions. They form in directions aligned to the grid and occur because the shock is not aligned with the grid everywhere. Since our solutions are idealized and self-similar, these results apply to any time after the solution has converged to a statistically self-similar state. For a typical Type Ia SN, they apply to any time between 30 and 10 4 yr. They have the shape of a broken power law, depicted by the gray dashed lines. The vertical line is at the wavenumber where the power law is broken, denoted by l 0 . Since our solutions are idealized and self-similar, these results apply to any time after the solution has converged to a statistically self-similar state. For a typical Type Ia SN, they apply to any time between 30 and 10 4 yr. characteristic velocity of the ejecta). In this case, the density scale height is time-dependent:

This leads to an effective value of n (which, from Equation (9), may be interpreted as the ratio between the ejecta radius and the density scale height), which depends on time as

where the radius of the reverse shock is considered to be the relevant radius. Warren & Blondin (2013) obtain results that are consistent with l 0 ≈ 10n eff (corresponding to a dominant angular scale smaller than we find), as pointed out by Polin et al. (2022). This suggests that although the dominant angular scale of SNRs may be dependent on the nature of the density profile of the outer ejecta, such a scale should be expected to arise out of the hydrodynamic interaction between the ejecta and the surrounding CSM. In fact, works starting with 3D initial conditions from thermonuclear SN explosion models (Ferrand et al. 2019(Ferrand et al. , 2021)), instead of 1D analytical models such as ours, also find a similar break wavenumber or harmonic for their SNR models. The power spectrum increases as a steep power law (C l ∝ l 3 ) for l < l 0 (large angular scales), and then drops off again as a steep power law (C l ∝ l -2.8 ) for l > l 0 (small angular scales). The smallscale behavior of our models is somewhat different from that found by Polin et al. (2022) for their 2D axisymmetric models (C l ∝ l -3.5 ), or from that obtained by Warren & Blondin (2013) for 3D models that have ejecta with an exponential density profile (C l ∝ l -3.9 ). Nevertheless, our power-law scaling is too steep to be consistent with a turbulent cascade (Kolmogorov 1941), in which case the power law should vary as C l ∝ l -5/3 . We note that a power spectrum with scaling of l -3 would be consistent with a model where all scales shared the same eddy turnover time. We therefore confirm that RTI structures formed in young SNRs are not likely to exhibit an energy cascade.

Time Evolution

Polin et al. (2022) put forward the hypothesis that the power spectrum may instead be due to independent growth and saturation of RTI at all scales. We test this hypothesis by looking at the time evolution of the power spectra for a typical SNR model, shown in Figure 5. The spectrum at early times is relatively flat (except for multipole moments described earlier in Section 4), but gradually peaks at l 0 with time, with a steep power law (C l ∝ l -2.8 ) for large l or small scales. Since the power spectrum should scale in general as

where δv is the typical speed of eddies at wavenumber l, it follows that for our small-scale power law

0.9 0.1

where we have used the fact that the eddy size λ is inversely proportional to l and τ is the eddy turnover time. Hence, the eddy turnover time is almost the same for eddies of all sizes with l > l 0 . This is in contrast with Kolmogorov turbulence, where the eddy turnover time scales as τ ∼ l -2/3 . In the latter case, small-scale eddies live longer, allowing energy to be transferred to them from larger scales, where most of the energy initially resides. This is not possible for our results, since all eddies share roughly the same lifetime. Thus, we may not be seeing a turbulent cascade. Instead, each angular mode in our model could be experiencing its own independent net All curves peak at l/n ∼ 4, demonstrating the peak wavenumber l 0 is proportional to n. As these represent idealized self-similar solutions, these results apply to any time after the solution has converged to a statistically self-similar state. For a typical Type Ia SN, they apply to any time between 30 and 10 4 yr. Figure 5. Evolution of the velocity power spectrum for the n = 9, s = 0 solution. The power spectrum is evaluated for 12 logarithmically spaced intervals between t = 40 days and 30 yr (or t/t i = 40-10,000 in code units). The red and blue curves denote the earliest and the latest instances, respectively. An initially flat spectrum evolves into the peaked structure, indicating maximum net growth at the break frequency.

growth rate, dictated by the growth and saturation of RTI at that mode. The dominant mode l 0 therefore has the largest net growth rate. We also note that the net growth rate for modes l < l 0 is smaller. Fluctuations at these modes persist for much longer than at modes l > l 0 , where initial perturbations are drowned much faster. We find that the SNR has to expand by ∼40 times before developing a steady-state velocity power spectrum. The density spectrum takes much longer to reach a steady state, requiring expansions by ∼100 and ∼200 times for the l > l 0 and l < l 0 modes, respectively. These expansion factors can be compared to a typical young SNR, say Tycho's SNR, which has an age of ∼450 yr. It is estimated to have a radius of the order of a few parsecs. Assuming a white dwarf progenitor with radius 10 4 km for the Type Ia SN responsible for Tycho's SNR, we find an expansion factor of >10 10 for the radius, which should be sufficient for developing steady-state angular power spectra.

Clumpiness in the Ejecta

Numerical studies (Wang & Chevalier 2001;Orlando et al. 2012) have shown that anisotropies in supernova ejecta can significantly affect the behavior and size of Rayleigh-Taylor fingers and thus leave their imprint on the SNR structure. In fact, observations of clumpy metal distributions in some remnants suggest this is a likely case, e.g., SN 1006 (Winkler et al. 2014). In this section, we discuss the results of our investigation on the survival of density perturbations in the supernova ejecta of our models. We add the following density perturbations to the ejecta in our n = 9, s = 0 solution:

shows the effect of varying amplitudes (A) for perturbations at harmonics of low spatial frequency (k = 6) and high spatial frequency (k = 100). The behaviors of both the velocity and the density spectra at large wavenumber (l > l 0 ) are essentially independent of the nature of fluctuations in the ejecta density profile. The density spectra show a deviation from the fiducial (A = 0) spectrum at small wavenumbers (l < l 0 ), especially for A = 1.0. This is consistent with our observation of the time evolution of the power spectra, where we find that the small-wavenumber modes persist for longer durations. The velocity spectrum on the other hand, is only mildly affected, even for small l. It is interesting to note that even for small-scale anisotropies (k = 100), the fluctuations only show up at small wavenumbers. We do not yet have an explanation for this behavior. Overall, our findings would imply that the small-l component of the density spectrum is acutely sensitive to anisotropies in the ejecta.

Clumpiness in the CSM

The departure from spherical symmetry for many SNRs has also been attributed to a clumpy surrounding medium, as opposed to anisotropies in the explosion itself. In the case of Cas A, Sato et al. (2018) and Orlando et al. (2022) argue that the radially inward velocity in parts of the reverse shock and the offset between forward and reverse shock centers require interaction of the SNR with an asymmetric environment. Fang et al. (2018) propose that anisotropies in the surrounding stellar wind may be able to account for anisotropies seen in Tycho's SNR. Here we systematically investigate the effects of such clumpy CSM on our results. We again introduce density perturbations of the form ( 14), but in the CSM instead of the ejecta, and calculate solutions for n = 9, s = 0 as in Section 5. The values of δ and k remain the same as before. The resulting spectra are shown in Figure 7. We find the nature of the anisotropies in the CSM does not affect the velocity power spectra irrespective of the magnitude of density fluctuations. The density spectra are also affected very mildly, unlike the case of clumpy ejecta. The only noticeable departure is seen when the density perturbations are of the order of unity, as found previously by Polin et al. (2022).

Figure 8 shows the square of the density of the shocked plasma from a few different models, integrated along the line of sight. The images are examples of a first approximation of the thermal emissivity map of SNRs, assuming that emission is produced only by the shocked plasma and that the plasma is optically thin. The left half shows images for models with large-scale initial perturbations (with l = k = 6) and illustrates the difference between structures caused by large external clumps and those caused by large internal clumps. It can be seen the large-scale structures are only formed when the clumps are massive (bottom panels, left half). Also, the forward shock can be perturbed when there are large-scale massive clumps in the CSM. In contrast, no such large-scale structures or shock perturbations are seen in the right half of Figure 8, which shows thermal emission from models perturbed with small-scale clumps (l = 100). The dominant angular scale is the same for those models, and is solely determined by the density scale height.

Discussion

We apply the technique of power spectral analysis to 3D hydrodynamic numerical models of young supernova remnants to probe the nature of anisotropies they exhibit. We study the impact of seed anisotropies or clumps in the CSM for our models to verify the well-studied role of these CSM clumps in generating the turbulent structures in SNRs. In addition, we perform a systematic parameter study of anisotropies endemic to the SN explosion to determine the conditions under which they can impact the SNR structure. We confirm that the power spectrum of the spatial distribution of either density or radial velocity of young SNRs has the form of a broken power law, as found previously (Warren & Blondin 2013;Ferrand et al. 2019;Polin et al. 2022). The break (or peak) harmonic l 0 represents the angular scale where most of the power resides. In summary, we find the following characteristics for our 3D calculations:

1. The dominant angular mode (l 0 ) of the power spectrum is related to the power-law index (n) of the ejecta density profile as l 0 ∼ 4n. Thus, this is a direct diagnostic of the density profile of the outer layers of the SN ejecta, as pointed out by Polin et al. (2022). However, the dominant angular mode we obtain is different from that obtained by Polin et al. (2022) (l 0 ∼ 3n) or Warren & Blondin (2013) (l 0 ∼ 10n). More studies are required to find an agreement on what dominant angular mode is to be expected for observations of SNRs. 2. The power spectrum falls off as C l ∝ l -2.8 at large l (or small angular scales). This is somewhat shallower than the scalings found by Warren & Blondin (2013) and Polin et al. (2022). The small-scale behavior of the power spectrum is found to be very consistent across all of our models, irrespective of the size, strength, or location of seed perturbations introduced. 3. Consistent with previous findings, the power spectra of our 3D SNR models are minimally affected by perturbations or clumpiness in the surrounding medium. The only impact is seen at small wavenumbers (l < l 0 ) in the density spectrum, as in the case of clumpy ejecta, but for very strong perturbations. 4. In contrast with the previous point, we find that density power spectra of SNRs are quite sensitive to large-scale (l < l 0 ) asymmetries in the ejecta. They are found to impart a long-lasting imprint upon the density power spectrum. Fluctuations at small scales subside quickly. The velocity spectra are even less affected irrespective of the angular mode of fluctuations.

In summary, our study confirms known or suggested properties of the angular power spectrum of SNRs, and establishes their connections to intrinsic properties of the ejecta. Although most of these connections were pointed out by Polin et al. (2022) in the same form, 2D numerical studies of turbulence are not necessarily accurate descriptions of turbulence seen in nature, as mentioned earlier. Warren & Blondin (2013) identify the dominant angular scale in their 3D models as solely a function of a scaled age for a few different times, but without an overall prediction for how this power spectrum should evolve with time generally. Hence, it is difficult to establish a direct connection of the observed dominant angular scale for a given SNR to the ejecta properties based on this model. Our discovery thus provides a path to infer density profiles of ejecta interacting with the reverse shock for observations of young SNRs by constructing their angular power spectra, in addition to the prediction that the presence of low harmonic modes in the power spectra is a likely indicator of anisotropies endemic to the SN itself, rather than the CSM.

There are other, arguably more realistic choices for SN ejecta models, ranging from spherically symmetric analytical profiles to 3D profiles obtained from numerical solutions of SN explosions. A preferred choice for modeling Type Ia SN ejecta (Dwarkadas 2000;Warren & Blondin 2013) is the exponential density profile (Dwarkadas & Chevalier 1998), designed to approximate the W7 carbon deflagration model for thermonuclear SN explosions (Nomoto et al. 1984). 3D density profiles, obtained from numerical calculations of thermonuclear SN explosions (Seitenzahl et al. 2013;Fink et al. 2014;Tanikawa et al. 2018), are also used to model SNRs. One of the advantages in this case is that one can incorporate realistic large-scale anisotropies in the ejecta as a part of the initial conditions (Ferrand et al. 2019(Ferrand et al. , 2021)). Another point worth mentioning is that our models do not take into account the cooling down of the shocked region via cosmic-ray acceleration. A common way to mimic shock cooling in hydrodynamic calculations is to increase fluid compressibility by reducing the adiabatic index (Blondin & Ellison 2001;Warren & Blondin 2013). Ferrand et al. (2010) explicitly study the impact of cosmic-ray cooling on instabilities in SNRs using hydrodynamic models that include a kinetic prescription for particle acceleration at the shock front. The presence of magnetic fields in the SNR is another important factor that could potentially impact our results. For instance, it has been shown that magnetic tension can suppress Rayleigh-Taylor instability, depending on the field strength and orientation (Mac Low et al. 1994;Jones et al. 1996;Bucciantini et al. 2004;Fragile et al. 2005). On the other hand, Porth et al. (2014aPorth et al. ( , 2014b) ) find that magnetic fields have minimal effect on the turbulent structure in their numerical models for pulsar wind nebulae, due to strong magnetic dissipation and the inability of magnetic fields to suppress unstable modes that exchange fluid between the field lines (Stone & Gardiner 2007). A systematic examination of models that incorporate these features is thus required to understand their impact on our results as well as observations and will be pursued in a future study.

Our findings indicate that large-scale features in the density profile of young SNRs, being acutely sensitive to perturbations Figure 6. Velocity (left) and density (right) power spectra for models seeded with perturbations in the ejecta according to Equation (14). The velocity spectra are not heavily impacted by even the strongest of perturbations. The density spectra do show a deviation from the fiducial spectrum (for the n = 9, s = 0 model), but only at large angular scales (or small l) for perturbations of large amplitude. As these represent idealized self-similar solutions, these results apply to any time after the solution has converged to a statistically self-similar state. For a typical Type Ia SN, they apply to any time between 30 and 10 4 yr. in the ejecta, are better tracers of the inherent anisotropy of the explosion, and therefore are not predictable by Rayleigh-Taylor instability alone. The large-scale structures in observed remnants may be indicative of the angular scale of the traces of perturbative activity (such as formation of 56 Ni clumps) in the SN ejecta under consideration. On the other hand, small-scale features in the density and velocity field are robust to fluctuations, both in the SN ejecta and in the surrounding medium. Ferrand et al. (2021) come to a very similar conclusion. They find that the low l modes in the SNR power 14). The velocity spectra are not influenced by perturbations, as with anisotropic ejecta. The density spectra are only affected by the strongest perturbations. As these represent idealized selfsimilar solutions, these results apply to any time after the solution has converged to a statistically self-similar state. For a typical Type Ia SN, they apply to any time between 30 and 10 4 yr.

Figure 8. The square of the density of the shocked plasma integrated along the line of sight. These images were developed using our models with n = 9, s = 0, and different initial perturbations to the density field, assuming the plasma is optically thin and is heated to over 10 6 K in the shocked region. Left: perturbations for these models are large-scale (corresponding to l = 6) and are either internal or external to the SN explosion. The quantity A is as defined in Equation (14). It can be seen large structures persist in the SNR only if the clumps are massive (bottom panels), in addition to being large. Right: the initial density perturbations for these models are on much smaller length scales (l = 100). No clumps of the corresponding size are seen in the SNR morphology, irrespective of the mass. Therefore in each case, the dominant angular scale seen is purely a function of the density scale height. As these represent idealized self-similar solutions, these results apply to any time after the solution has converged to a statistically self-similar state. For a typical Type Ia SN, they apply to any time between 30 and 10 4 yr. spectra are dominated by asymmetries endemic to the SN explosion (used as 3D initial conditions for their SNR models), whereas large l (small angular scale) modes form due to RTI alone. We therefore expect the break wavenumber l 0 , and the power-law index for large values of l, to be largely consistent between SNRs, provided that extreme events such as jets, highly asymmetric explosions, or collisions with a companion do not overwhelm the small-scale structure. We note that even though l 0 may not be the peak wavenumber due to the presence of large-scale features, the power spectrum is nevertheless expected to have a break around l 0 . This motivates a search for dominant angular scales in the morphologies of observed SNRs. Some work has already been done in this direction. For example, power spectra have been developed from resolved X-ray images of Tycho's SNR using the azimuthal distribution of the rotation measure (Lerche & Caswell 1979) and measures of remnant radius (Warren et al. 2005). A different direction was taken by Shimoda et al. (2018), who construct correlation functions of the synchrotron emissivity maps of Tycho's SNR. Lopez et al. (2009aLopez et al. ( , 2009bLopez et al. ( , 2011) ) analyze resolved Chandra images of a sample of nearby SNRs using multiple techniques, including a 2D multipole expansion, two-point correlation analysis, and wavelet transforms. They demonstrate that corecollapse SNe exhibit more large-scale asymmetry (stronger low l modes) than Type Ia SNe in general. It is perhaps worthwhile to mention here that an approach similar to the wavelet transform, called the Δ-variance method, was developed by Arévalo et al. (2012) to extract power spectra for 2D images as well as 3D data sets. However, most of these works start with 2D images, which inherently suffer from projection effects and therefore may not be able to provide accurate three-dimensional information. This problem may be circumvented by deprojecting observations to obtain estimates of the 3D density and velocity distribution of the SN ejecta, or by using highresolution imaging and spectroscopy to directly estimate velocities (DeLaney et al. 2010;Williams et al. 2017;Seitenzahl et al. 2019). Therefore, more work on observations and analysis techniques is needed to learn about the underlying physical properties of SN explosions and progenitors by comparing power spectra of observations to those of theoretical models.

Appendix Numerical Convergence Studies

Here we describe the impact of spatial resolution and the duration of numerical calculations on our results. Figure 9 shows the velocity power spectrum (without additional smoothening) obtained from one of our fiducial models (n = 9, s = 0) at four different spatial resolutions. It can be seen that spectra for all resolutions agree with each other in their general shape and location of the break wavenumber. Moreover, they all have a distinct feature at l = 4. As previously mentioned, we identify this feature, also seen in Figure 2 as square-like patterns, with a numerical instability (called the carbuncle instability) that arises due to the shock (spherical here) not perfectly aligning everywhere with the grid geometry (Cartesian in this case). Nevertheless, the spurious mode at l = 4 introduced by this instability is seen to decrease with resolution. No other similar effect is seen on the power spectrum that decreases with resolution and that can therefore be identified with grid effects or the carbuncle instability. A similar phenomenon was observed by Fraschetti et al. (2010) for their 3D numerical models of SNRs using a Cartesian grid. They used adaptive mesh refinement to resolve the shock structure in the SNR and found a similar agreement for the shock structure between the lowest and highest levels of resolution.

The effect from the grid geometry on the power at different modes is the most prominent at early times, as exhibited by Figure 5. At t ∼ 40 days (or t/t i = 40), the spectrum is mostly flat except for features at low values of l (l  10), where it peaks. These features are indicative of grid effects. Even though initially these features have power roughly a couple orders of magnitude greater than power at any other scale, turbulent activity quickly introduces power at all modes (with the largest growth for l = l 0 ). The power spectrum is seen to reach its steady-state shape after nearly three decades of time have elapsed, and saturation of the power spectrum is attained after four decades. This is in agreement with Warren & Blondin (2013), who find that running calculations for at least four decades of time is required for the effects of the initial perturbations to get washed out by saturation of the instabilities.