© The Authors 2023

1 Introduction

Supernova remnants (SNRs) are considered a significant actor in the acceleration of galactic cosmic rays (CRs; Blandford & Eichler 1987). Observations of non-thermal X-ray synchrotron emission originating from the rims of SNRs (Koyama et al. 1995) led to the hypothesis of an acceleration of high-energy CR electrons at a shock front. As CR protons are hard to observe, measuring CR acceleration efficiency is one of the primary interests of this research field. The proximity between blast wave (BW) and contact discontinuity (CD) can be a probe of high-energy CR acceleration efficiency. Indeed, an efficient CR acceleration lowers the post-shock gas pressure, thus narrowing the distance between BW and CD. The first analysis of this distance established the BW as the main stage for CR acceleration (Warren et al. 2005). An average radius ratio of 1:0.93:0.70 for the respective BW, CD, and reverse shock (RS) was measured for the type Ia SNR, Tycho (SN 1572), which is nearly 450 yr old. The proximity between BW and CD (0.93) could not be explained by the one-dimensional (1D) hydrodynamic simulations performed at that time. This led to the current consensus that a significant fraction of the BW energy is lost in CR acceleration.

Since this first work, many studies pertaining the relative position of BW, CD, and RS in SNRs have been performed. Observations of Tycho (Warren et al. 2005; Cassam-Chenai et al. 2007), SN 1006 (Hamilton et al. 1986a; Cassam-Chenaï et al. 2008), Cassiopeia A (Gotthelf et al. 2001; DeLaney et al. 2004; Lee et al. 2014), and so on have led to similar conclusions. Similarly, more complex simulations have been performed to recreate the morphology of these SNRs. The first major process considered was the impact of multi-dimensional simulations on the problem. Ferrand et al. (2010), and Fraschetti et al. (2010) show that the development of interfacial instabilities at the CD, mainly Rayleigh-Taylor instability, leads to greater proximity between CD and BW than predicted with 1D simulations. However, this improvement was not sufficient to explain observations. Other parameters have also been studied. For instance, the nature of the supernova (SN) explosion model influences the initial structure and symmetry of the remnant, and these have an important impact on the simulated evolution of the young SNR. In particular, the asymmetry linked to the explosion remains in the remnant for its first few millennia of evolution. Other physical processes, such as the fission of the heavy elements or the presence of a magnetic field (in magnetohydrodynamic (MHD) simulations), also affect the evolution of the SNR but on a more localised scale (Orlando et al. 2021). For instance, the fission of the heavy elements leads to a local injection of energy, which results in a quicker expansion of the ejecta in local bubbles (Ferrand et al. 2019; Orlando et al. 2021). This logically depends on the distribution of these elements, which is asymmetric. The magnetic field tends to stabilise the CD, limiting the growth of instabilities, without significantly changing the relative position of the CD or BW.

To date, the best simulation results in reproducing the observed radius ratio were achieved by including the back reaction of CR acceleration on the dynamics of the fluids and shock. In particular, the energy loss at the shock front due to CR acceleration should decrease the shock compression ratio. This lower compression ratio results in closer proximity between CD and BW. In order to reproduce this phenomenon in hydrodynamic simulations, it became standard procedure to use an effective adiabatic index for the shocked interstellar medium (ISM; Ferrand et al. 2010; Fraschetti et al. 2010; Orlando et al. 2012). This effective adiabatic index depends on the particle-injection rate as well as on the history of the SNR, or more simply its age, as shown by Fig. 2 of Orlando et al. (2016). Therefore, the adiabatic index used in hydrodynamic simulations may evolve from 5/3 (typical unshocked ISM) to 4/3 or 3/2 at the BW front. This artificial correction leads to a higher shock compression and therefore to a thinner shocked ISM shell. While the results obtained through this method are comparatively closer to observations, they still do not reach the correct proportion (Ferrand et al. 2010).

Despite all these considerations of the initial state of the SNR and the physical mechanisms that may apply, few studies focus on the effect of the ISM. Usually, spherical symmetry is supposed for the density of the ISM, using either a constant or a power law, ρ ∝ r⁻², in the case of stellar wind (Ostriker & McKee 1988; Truelove & McKee 1999; Tang & Chevalier 2017). While such profiles are realistic from an average point of view, as shown in the study of Cassiopeia A by Lee et al. (2014), they do not reflect the reality of the local scales. As the ISM and circumstellar medium are turbulent (Ferrière 2020), local density variations are to be expected. This should affect the dynamics of the SNR both in terms of symmetry and BW-to-CD ratio. Indeed, the propagation of a shock wave inside turbulent medium results in amplification of the turbulence (Inoue et al. 2012, 2013) and energy loss (Chashei & Shishov 1997; Fang et al. 2020). Therefore, we can expect the BW to be affected by the turbulence of the ISM in a fashion similar to the feedback of CR acceleration. In addition to these hydrodynamic effects, some MHD effects should also exist, as the compression of a turbulent magnetic field by a shock wave tends to amplify it. We will study this effect in future works.

In this article, we present a study of the effect of a turbulent ISM on the evolution of young SNRs from a purely hydrodynamic point of view. To limit the effect of external factors on our study, we neglect most of the physics usually associated with SNRs (CR acceleration, asymmetry in the explosion, etc.) and use a simplified profile to initialise our simulations.

This article is organised as follows. In Sect. 2, we describe our simulation setup. This includes the simulation code we used, the equations that are solved, the initialisation of the turbulence, and the initial profiles of the SNR. In Sect. 3, we present our analysis method and the results of our simulations, which we compare to SNR observations. Finally, we draw our conclusions in Sect. 4.

2 Methods

2.1 Simulation code

The starting point of our study is Warren et al. (2005), which focuses on the nearly 450-yr-old type Ia SNR Tycho. These latter authors used a 1D Lagrangian code. All materials were described with a mono-atomic adiabatic equation of state with a constant adiabatic index (γ = 5/3). The simulations were initialised sometime after the SN explosion (~2 weeks) using profiles either from simulations of SN explosions or self-similar SNR expansion models. We should note that each SN explosion model used in Warren et al. (2005) results in a slightly different BW:CD:RS radius ratio, none of which reproduced the observations.

We based our study on the same kind of consideration. We employed the FLASH4 code (version 4.6.2) in a purely hydrodynamic configuration. We used ideal gas equations of state with a single constant adiabatic index equal to 5/3. Contrary to the above-cited study, we focus on the effect of a turbulent ISM on SNR evolution, and therefore use three-dimensional (3D) simulations. As a result, hydrodynamic instabilities develop at the CD and the BW may be deformed. The initial state of our simulations is also different from Warren’s, as our ISM is non-uniform.

Contrary to Warren et al. (2005) Lagrangian simulations, we employed a Eulerian code. This comes with its own set of advantages and drawbacks. Eulerian codes are globally simpler in 3D compared to their Lagrangian counterpart. For instance, there is no mesh tangling problem. However, the interfaces between fluids such as the CD are harder to define. In our simulation, we use a passive scalar initialised at the ejecta position to follow this interface. However, this passive scalar tends to diffuse numerically due to the Eulerian nature of the code. This leads to greater uncertainty over the position of the CD.

To highlight the effect of the turbulence on the evolution of the SNR, we performed two types of simulations: with a turbulent ISM, or with a uniform ISM, all other parameters being the same. The non-turbulent ISM is used as a reference.

2.2 Hydrodynamic equation

One of the characteristics of the evolution of an SNR is its rapid expansion. In 500 yr, a young SNR can see its diameter multiplied by a factor of 5000. This brings an important difficulty for a Eulerian simulation, where one should balance the domain size of the simulation with its actual resolution. As our study deals with the relative position of the RS and the CD to the BW with an expected variation of the order of a few percent (0.86 simulated in Hamilton et al. 1986b against 0.93±2 measured by Warren et al. 2005), we should at least ensure a spatial resolution of higher than 1% of the radius of the SNR at any time. This can be achieved by a direct approach in 1D, but the computational cost in 3D is prohibitively expensive.

The usual approach to keeping a small total number of cells would be to employ an adaptive mesh refinement (AMR) scheme (Berger & Colella 1989; Fryxell et al. 2000; Orlando et al. 2016). However, AMR is not well adapted to simulations of a turbulent domain, as derefinements smooth the turbulence. Furthermore, the numerical cost of an AMR scheme remains important for the late-time evolution, as the SNR spans a large part of the domain. For this reason, we employed the expanding frame approach theoretically developed by Poludnenko & Khokhlov (2007). This approach consists of a change in the referential frame to a co-expanding frame, which is non-inertial. This approach has already been employed to study SNRs in a uniform ISM with the RAMSES code (Fraschetti et al. 2010).

As the expanding frame is non-inertial, the Euler equation should be modified. In our work, we used a non-rotating non-accelerating expanding frame, with ideal gas equations of state. We set the scaling parameters of Poludnenko & Khokhlov (2007) to simplify the equation: α = ν = 3 and β = 1 in their notation. As a result, the following equations were numerically solved: $${\partial }_{\tau }\tilde{\rho }+\tilde{\nabla }\cdot \left(\tilde{\rho }\tilde{u}\right)=0$$(1) $${\partial }_{\tau }\left(\tilde{\rho }\tilde{u}\right)+\tilde{\nabla }\cdot \left(\tilde{\rho }\tilde{u}\otimes \tilde{u}\right)+\tilde{\nabla }\tilde{P}=0$$(2) $${\partial }_{\tau }\tilde{E}+\tilde{\nabla }\cdot \left(\tilde{u}\left(\tilde{E}+\tilde{P}\right)\right)=\frac{{\upsilon }_0}{1-{\upsilon }_0\tau }\left(2-3\left(\gamma -1\right)\right)\tilde{\rho }\tilde{e}$$(3) $$\tilde{P}=\left(\gamma -1\right)\tilde{\rho }\tilde{e}.$$(4)

This set of equations corresponds to the Euler equations after applying the following transformation: $$\tilde{r}=a^{-1}r$$(5) $$\tau ={\displaystyle {\int }_0^t\frac{\text{d}t\prime }{a^2\left(t\prime \right)}}$$(6) $$\tilde{\rho }=a^3\rho$$(7) $$a\left(t\right)=1+{\upsilon }_{0}t,a\left(\tau \right)=\frac{1}{1-{\upsilon }_0\tau },$$(8)

where υ₀ is a constant corresponding to the non-inertial frame expansion rate, r and t are the space-time coordinate in the inertial frame, and $\tilde{r}$ and τ are their equivalent in the expanding frame. From this transformation, we also have: $$\tilde{P}=a^{5}P$$(9) $$\tilde{e}=a^{2}e$$(10) $$\tilde{u}=au-{\upsilon }_{0}r$$(11) $$\tilde{E}=\tilde{\rho }\tilde{e}+\tilde{\rho }{\tilde{u}}^2/2,$$(12)

where ρ is the density, P the pressure, u the fluid velocity, and e the internal energy, all expressed in the inertial frame. The tilde values are their equivalent in the expanding frame.

Except for the velocity, this transformation does not change the initial profiles of the simulations, as a(t = τ = 0) = 1. Furthermore, only the energy-conservation equation is changed to include a source term. However, care must be taken as this transformation leads to some other modifications that are not conspicuous in the above equations. As there is a time-dependent renormalisation of mass, space, and time, all physical constants depending on these physical quantities are also time dependent. For instance, in FLASH4, the temperature of the fluids is calculated at each time step using the internal energy: $\tilde{e}=N_{\text{a}}\tilde{k}T/\left(\left(\gamma -1\right)A\right)$, with N_a being the Avogadro number, A the average atomic mass, and k the renormalised Boltzmann constant. As $\tilde{k}=a^{2}k$, the Boltzmann constant is time dependent.

The expanding frame also changes the domain boundaries. As the simulation domain expands from the inertial frame point of view, there is an influx of matter from the expanding frame point of view. This influx shapes the ISM in our simulations. Therefore, all the information related to the ISM and its turbulence is input into the domain ghost cells at each time step. Poludnenko & Khokhlov (2007) highlights that this inflow can result in noise-like features on the velocity field when there is a mismatch between the velocity inputs into the ghost cells and the velocity as extrapolated from the simulation domain. Such features are negligible in our simulations as the expansion rate is constant (no noise was observed in our simulations).

2.3 Turbulence

As we employ an expanding frame to maximise the size of the SNR compared to the simulation domain, the ISM flowing in at the boundary interacts quickly with the BW. The use of a stirring force to generate turbulence is therefore insufficient, as time is required for the turbulence to appear. For this reason, we initialised both our ISM and its inflow with turbulent fluctuations.

In our simulations, only the density field of the ISM is initialised to follow a turbulence spectrum as observed in Armstrong et al. (1995). The velocity field, which should also present local fluctuations, is taken to be null in the inertial frame. This choice is not realistic as there is a correlation between density and fluid velocity. However, the exact details of this correlation are still unknown, and so the initialisation of both fields by hand would also bring unphysical conditions. We considered a constant temperature for the ISM, and we calculate its pressure using the perfect-gas equation of state.

We initialised the density in the inertial frame as: $$\rho \left(r\right)={\rho }_0\left(r\right)\left(1+\delta \overline{\rho }\left(r\right)\right)$$(13) $$\text{with\hspace{1em}}\delta \overline{\rho }\left(r\right)={\displaystyle \sum _{k\ne 0}\sqrt{C{k}^{-\alpha }}\sin \left[2\pi \left(k\cdot r_L+\text{Φ}\left(k\right)\right)\right],}$$(14)

where the density, ρ, is composed of its average, ρ₀, and dimensionless fluctuating part, $\delta \overline{\rho }$. The coefficient of the turbulent power law, α, is taken to be equal to 11/3 in order to correspond to the 3D Kolmogorov spectrum (Kolmogorov 1991). Here, r_L is the coordinates renormalised to the field periodicity: r_L = {x/L_x; y/L_y; z/L_z}. k ∈ ℤ³ is a wavevector with a norm k. In our simulations, we limit k to the subset 〚−K; K〛 for numerical reasons (here K = 128). Φ is a phase, which we randomly initialised.

The implementation of the above formula by a direct calculation is numerically expensive for a 3D domain as the number of operations scale as 𝒪(K³N³), with N being the number of cells per direction per simulation block. Even if the initialisation hurdle is overcome, similar calculations are needed at each time step to fill the ghost cells. To limit the numerical cost of turbulence initialisation, we based our rationale on the signal-processing principle, according to which the quantity of significant data output is equal to the quantity of data input. Accordingly, only a number of points equal to the number of input phases is significant in defining the density variation. Thus, we solve Eq. (14) using a fast Fourier transform (FFT) algorithm for a 𝒪(K³ log K) complexity. As a result of the FFT calculation, we obtain the density variations on a grid of K³ cells, which we hereafter refer to as FFT grid and FFT cells, respectively. To obtain the density value on the simulation grid, we linearly interpolate each simulation cell with its neighbours on the FFT grid based on their respective spatial coordinates.

When using an expanding frame, further considerations are required to generate a turbulent inflow of matter. Indeed, the size of the cells varies during the simulations, including the ghost cells. As a result, a simple linear interpolation using the eight neighbouring FFT cells cannot be used. Indeed, as the simulation cells become larger due to the expanding frame, they can contain more than eight FFT cells at late times. This is also due to the fact that the size of the FFT grid cells is constant. The density of a specific simulation cell should therefore be the average of the density of each FFT cell it contains, or more precisely a volume-weighted average of the intersection of all FFT cells with this specific simulation cell. We achieved this average through the integration of the FFT cell values, which we performed before the FFT calculation. Thus, Eq. (14) is reformulated as: $$\delta \overline{\rho }\left(r\right)={\displaystyle \sum H\left(k\right){\text{e}}^{2\pi ik\cdot r_L}}$$(15) $$H\left(k\right)=H_0\left(k\right)H_1\left(k\right)$$(16) $$H_0\left(k\right)=-\frac{i}{2}\sqrt{C{k}^{-\alpha }}\left({\text{e}}^{2\pi i\text{Φ}\left(k\right)}-{\text{e}}^{-2\pi i\text{Φ}\left(-k\right)}\right)$$(17) $$H_1\left(k\right)=\mathrm{sinc}\left(\pi {\text{Δ}}_x{k}_x/L_x\right)\mathrm{sinc}\left(\pi {\text{Δ}}_y{k}_y/L_y\right)\mathrm{sinc}\left(\pi {\text{Δ}}_z{k}_z/L_z\right),$$(18)

where H(k) corresponds to the FFT coefficient with H₀ directly related to Eq. (14) and H₁ the consequence of the integration. Here, sinc is the cardinal sine defined as sinc(x) = sin(x)/x, and Δ_x (respectively Δ_y, Δ_z) correspond to the width of the cells in the x direction (respectively y, z) in the inertial frame.

From an algorithmic standpoint, the FFT coefficient and FFT should be recalculated at each time step. The use of an FFT-type algorithm limits the cost of this coherent implementation of turbulence in expanding-frame simulations.

The above implementation results in a correct spectrum from the standpoint of the Kolmogorov theory, but also in some unphysical values. Indeed, $\delta \overline{\rho }$ can reach values of lower than −1. For this reason, we defined a lower cutoff value, 1 × 10⁻³⁰ g cm⁻³, to which we artificially set all lower density values.

The strength of the turbulence is set with C; however, a singular value of C does not hold a significant meaning without consideration of the distribution of $\delta \overline{\rho }$. This distribution is nearly Gaussian, has an average value of 0 and its standard deviation, σ, depends on C and the 256³ random phases. In this article, we mainly set C such that σ = 0.5, as this value is realistic in this situation. Indeed, given a density perturbation of the order of the mean density at the driving scale (L_drive ~100 pc), and the turbulence scaling as ℓ^1/3, with ℓ being the spatial scale, we should have density variations of the order of half the average density at around 10 pc (nearly the size of our FFT grid). We also studied some smaller values of σ to see the effect of the strength of the turbulence.

2.4 Simulation parameters

In our simulations, we work with a Cartesian domain and use the unsplit hydrodynamic solver of FLASH4. We employ an HLLC-type Riemann solver with a second-order spatial reconstruction and a minmod-type slope limiter. The simulations are performed on a cubic domain, with an initial length of 2 × 10⁻³ pc. The majority of the simulations were performed with 300³ cells, but some were redone with higher resolution, namely 520³ cells. The simulation domain contains one-eighth of the SNR (an octant). As a result, half of our boundary conditions are reflective (SNR side), while the other half are inflow for the expanding frame (ISM side). The use of reflective boundary conditions allows higher resolution with no additional computational cost. However, artefacts may be produced near these reflective boundaries.

To initialise the SNR, we use a 1D profile supposing a spherical symmetry. The profile is piston-like, consisting of radially linear velocity, going from 0 cm s⁻¹ at the centre to 1.2 × 10⁹ cm s⁻¹ at the CD, and a constant density of ρ ~9.0 × 10⁻¹⁵g cm⁻³. The SNR, thus initialised, has a radius of 4.16 × 10¹⁵ cm, a mass of 1.37 M_⊙, and kinetic energy of 1.16 × 10⁵¹ erg. This corresponds to Tycho 10⁶ s after its explosion according to the model DDTc of Badenes et al. (2003, 2005), which was used as the initial state in the simulations of Warren et al. (2005). The profile we used is simpler than the one resulting from SN explosion simulations, which typically have anon-constant density. As such, our simulations are not meant to recreate the exact state of Tycho, but to help us to understand the effect of the ISM turbulence on the evolution of such an SNR. The ISM is initialised with no velocity and an average density of 1 × 10⁻²⁴ g cm⁻³. The temperature of the ejecta and the ISM is initially set to 10 000 K, which is negligible in comparison to the SN kinetic energy. These values for the density and temperature of the ISM correspond to a warm neutral medium.

The expanding rate, υ₀, was set by taking into account two constraints. First, the whole SNR, with the BW included, should remain in the simulation domain during the whole simulation. This adds a lower limit to the expanding frame velocity. Second, υ₀ should not be too large, as the SNR would then contract quickly. This contraction would result in a poorly resolved SNR for the late times. In our simulations, we set υ₀ to 2.7 × 10⁻⁷ s⁻¹.

We initialise the turbulence with sets of 256³ random phases. We use multiple sets to limit the impact a single one would have on our conclusions. We set the periodicity of the turbulence pattern to a single value, L_x = L_y = L_z, of the order of the size of the simulation box at the end of the simulation. This limits the smoothing and disappearance of the density fluctuations, which happen when the cell size becomes of the order of L_x. We set the periodic length to 8.1 pc for a 8.3 pc simulation domain after 500 yr.

3 Results and analysis

3.1 Simulations results

Figure 1 shows density slices in the YZ plane of the simulations corresponding to an SNR of 450 yr old. In this figure, the displayed values are taken in the inertial frame. As can be seen here, the presence or absence of turbulence has a low impact on the actual density and apparent radius of the remnant. This is not surprising as, on average, with or without turbulence, the situation is the same. However, a major difference exists between the two cases. Without turbulence, both the BW and RS appear almost spherical. This is usual in astrophysics as the theoretical answer to a punctual energy perturbation is a spherical shock wave. With turbulence, this spherical symmetry is broken for both BW and RS. This is a direct consequence of the variation of the shock velocity with the density of the medium where it propagates. When the BW encounters an over-density (respectively under-density) of the ISM, its velocity becomes lower (respectively higher). Thus, the shock front gets deformed, keeping an imprint of the encountered variations of density. As the BW interacts directly with the density fluctuations of the ISM, its deformations are the most prominent; whereas the density fluctuations seen by the RS are dampened, as they result from the interaction between ejecta and the ISM. The RS deformations are therefore more subdued.

Special attention must be paid to the CD as it is not spherical even in absence of turbulence. This deformation of the CD is not surprising; it is the result of the growth of the Rayleigh-Taylor instability (RTI). This hydrodynamic instability develops at the interface between fluids when the gradient of density opposes the inertial pseudo-force in the reference frame of the interface, and appears here due to the deceleration of the CD. This instability cannot be produced in 1D simulations and is therefore a source of divergence with Warren et al. (2005). In the case without turbulence, the RTI growing at the CD is relatively small and regular. In the turbulent case, the structures created are globally larger and more irregular. A visualisation of the whole simulation shows that deformation of the average position of the interface due to the turbulent variation of density leads to merging between neighbouring spikes or bubbles; thus, the turbulence affects more than simply the average position of the CD. We should note that the Richtmyer–Meshkov instability, which should also grow at the CDs of young SNRs, does not appear in our simulations. This is due to the initialisation of our simulation with the BW outside of the ejecta.

Fig. 1 Density slice taken 450 yr after SN explosion. The slices are performed on the YZ plane (x = 0) for both types of simulation: without turbulence (a), and with turbulence (b). The spatial coordinates and density values displayed here are taken in the inertial reference frame. Here, the main effect of the turbulence is the loss of the spherical symmetry of the SNR and a change in the instability wavelength and amplitude at CD.

Fig. 2 Radial average profiles 450 yr after the SNR explosion. The radial profiles of the density, ejecta, and norm of the velocity are displayed in the inertial coordinates both with (plain curve) and without (dashed curve) turbulent ISM. These profiles are averaged over all angles of the octant.

3.2 Analysis method and interpretation

Figure 2 shows the density, velocity, and ‘ejecta’ radial profiles associated with Fig. 1; they are obtained by performing an average over all angles. Here, the ejecta field corresponds to a passive scalar initialised at the position of the ejecta material. This passive scalar is transported by the fluid and is thus subject to both convection and numerical diffusion, but not to physical diffusion. A value of 1 corresponds to 100% ejecta and 0 to 100% ISM, the value in between showing a mix between the two materials. In the absence of turbulence, the 1D profiles in Fig. 2 give us a good idea of the positions of the three interfaces we are interested in. The discontinuities in velocity and density correspond to the BW and RS. The CD position is usually taken at the local density minimum between both shocks (Fraschetti et al. 2010); in our case, this corresponds to an ejecta’s value of 0.14. Hereafter, we use 0.1 as the value for the CD position. In the presence of turbulence, the loss of spherical symmetry leads to a broadening of all structures and a smooth transition in density at the CD. The shocks might appear broader in this figure, but this is due to the spatial average and not to a physical broadening of these shocks.

To perform a quantitative study of our simulations and compare our results to Warren’s observations, we have to consider the relative position of BW:CD:RS in each spatial direction. For this, we project our Cartesian simulation on a Spherical covering grid and then find the radius of the BW, CD, and RS for each solid angle. Figure 3 shows a line-out of the pressure and the ejecta for an arbitrarily chosen solid angle. The two jumps in pressure correspond to the BW and RS, and a value of 0.1 in the ejecta field corresponds to the CD. In our analysis, we used the pressure field to find both BW and RS, as the simplicity of this profile makes it easier to automate their measurement. We should add that in some situations, the position of each interface was not unique for some solid angles. This happens when the tip of the RTI spikes broadens in the non-linear phase of the instability. We removed these points from the following analysis as their interpretation is ambiguous. However, these points only have a small effect on the following analysis. For instance, for the 450-yr SNR with a turbulent ISM, the ratio of radius between CD and BW varies between 0.874(32) when keeping these points and 0.872(30) when removing them.

From the analysis method described above, we can associate each time step of the simulations with the distributions of BW, CD, and RS radii. Figure 4 shows the distributions associated with the Fig. 1 simulations. As can be seen, without turbulence, the BW and RS both have a sharp distribution (small standard deviation ~3 × 10¹⁶ cm for the BW) and are therefore spherical. The CD shows a broad distribution, with a standard deviation of 3.6 × 10¹⁷ cm for a mean value of 1.021 × 10¹⁹ cm. This is due to the RTI growth. Here the standard deviation of the CD distribution is a good indication of the scale of the width of the mixing zone, which is a source of uncertainty in astronomical observations. In the turbulent case, all distributions are broader because of the loss of spherical symmetry.

By normalising the interface radius by either the BW or the RS radius, we can partially compensate for the loss of spherical symmetry. The resulting distributions shown in the right panels of Fig. 4 are still broader in the turbulent case as the BW is more sensitive to ISM turbulence than the RS or CD, as explained above.

Fig. 3 Radial profiles for an arbitrary solid angle 450 yr after SNR explosion. The pressure and density profiles of the simulation with turbulence are displayed. The discontinuities in the pressure profile correspond to the BW and RS. The CD is defined as the point where the ejecta crosses 0.1.

3.3 Aspect ratio and RTI

Considering the normalised distribution obtained at each time step, we can calculate the temporal evolution of the ratios of BW:CD:RS and compare them to Warren’s observation. Figure 5 shows the temporal evolution of the average of the distribution of each ratio. The coloured zones correspond to the standard deviation of each distribution, which is due to RTI growth and/or partially compensated loss in spherical symmetry. The hatched zones correspond to the ratios observed by Warren. In a similar fashion to Warren et al. (2005), the observed BW:CD and BW:RS ratios are not reproduced simultaneously by our non-turbulent simulation. This is also the case for simulations with turbulent ISM when working only with the averages; however, the conclusion is reversed when taking into account the standard deviation. From these results, we can draw two conclusions. First, the ISM turbulence does not significantly change the radius ratios on average, and so in purely hydrodynamic simulations, turbulent density variations of the ISM and the resulting post-shock turbulence do not significantly change the effective adiabatic index. However, with turbulence, the RTI growing at the CD becomes large enough to reproduce the observed ratios. This brings into question the true nature of Warren’s observations. If the observations made by these authors are of sufficiently high resolution to see the deformation of the interface due to RTI, and if their analysis method correctly compensates for the 3D nature of Tycho, then the effect of the ISM turbulence alone is insufficient. However, if the resolution of the observation is too low, a bias in the measurement of the CD radius might appear, shifting the measured position of the CD in the direction of the RTI spikes. In this case, the ISM turbulence can explain the observation.

Figure 6 shows the temporal evolution of the CD standard deviation after spherical compensation by the RS radius (the CD radius obtained for each solid angle is divided by the equivalent RS radius and multiplied by the average RS radius $\langle \frac{r_{\text{CD}}}{r_{\text{RS}}}\rangle \langle r_{\text{RS}}\rangle$). As can be seen in this figure, the RTI appears to grow more quickly in the presence of turbulence. The temporal evolution of the CD average position is similar in the two cases, with a slightly higher deceleration in the absence of turbulence. Thus, the difference in RTI growth is not due to differences in the acceleration profiles of the interface, but mainly to the deformation of RTI (change in wavelength, merging of structures, etc.) and to the change in ambient medium (variation of Atwood number due to the varying ISM density and zone of higher and lower pressure due to the deformation of the shock).

This enhancement in RTI growth is reminiscent of previously performed 2D studies of the interaction of young SNRs with clumpy circumstellar medium (Jun et al. 1996). In these studies, the interaction of the BW with spherical clumps of over-densities present in the ISM results in the creation of eddies. Counter-rotating eddies form flow channels that boost the kinetic energy of well-placed RTI fingers, and thus their growth. Our simulations show that the average vorticity of the shocked ISM is four times higher in the turbulent case than in the non-turbulent one. However, rather than being well distributed in the shocked shell, the high-vorticity regions are concentrated along the RTI spikes at CD. This is particularly clear for the non-turbulent case where the gap between BW and CD is devoid of vorticity, while high-vorticity zones clutter the RTI. This suggests that rather than being created by the interaction between the BW and the density perturbations, the eddies are due to the growth of the RTI (Bian et al. 2020; Rigon et al. 2021). As mentioned above, in our case, the perturbations of density influence the RTI mainly through the changes in Atwood number and through the deformation of the BW, which results in high- and low-pressure regions (cf. mechanism of the Richtmeyer–eshkov instability). Contrary to the findings from studies of clumps, in the presence of turbulence, the changes in density are progressive (clumps have an interface with the ISM) and the post-shock eddies might not align in 3D as well as in 2D.

3.4 Discussion and limits

The turbulence present in our simulations enhances the growth of RTI and leads to BW:CD:RS ratios comparable to those seen in Warren’s observations (Warren et al. 2005). This conclusion does not depend on the specific simulation shown in this article. We performed additional simulations with different random sets of phases for the turbulence (see Appendix for additional curves). In all of them, the RTI growth is enhanced by the turbulence and there is a temporal range where both CD:BW and RS:BW ratios are reproduced simultaneously. However, the actual temporal range may be shifted (with some as early as 370 yr), and the RTI growth, while keeping the same concavity, slightly differs from case to case.

However, we should note that the strength of the turbulence is an important parameter. For instance, for σ = 0.05, and therefore a ten-times-lower density perturbation, the simulation results are close to those obtained in the simulation without turbulence. In particular, both the BW and RS are almost spherical. The standard deviation of the BW radius corresponds in this case to 0.9% of the average BW radius, compared to 0.2% without turbulence. Furthermore, both BW:CD and BW:RS ratios are not reproduced simultaneously, as in the non-turbulence case.

As in any other numerical study of the young SNR, the initial profile of the remnant plays an important role in its dynamics, even when neglecting all initial asymmetry. Therefore, one might question whether or not our conclusion still holds true with more realistic initial profiles. To attempt to answer this question, we performed a pair of simulations, with and without turbulent ISM, starting from the DDTc radial profile used in Warren et al. (2005). In particular, this profile presents a drop in density near the edge of the remnant (CD) and two pressure spikes near its centre. The dynamic of the resulting remnants is different from our piston case, and requires a higher expansion rate (9.2 × 10⁻⁷ s⁻¹). However, similar conclusions on the RTI growths and radii ratios can be drawn. More precisely, in the simulation that we performed based on a DDTc profile, the BW:CD:RS ratios could not be recreated on average (without taking the RTI into account) but are obtained only in the turbulent case when taking into account the standard deviation. In the specific simulation performed, both ratios were obtained simultaneously between 332 and 423 yr.

Fig. 4 Spatial distribution of the SNR radius. Histogram of the RS (blue), CD (orange), BW (green) radius both without (a) and with (c) turbulence obtained 450 yr after SN explosion (simulations displayed in Fig. 1). To compensate for the spatial variation of the BW radius appearing in the turbulent case, the RS and CD radii are normalised using their respective BW radius for each solid angle both without (b) and with (d) turbulence.

Fig. 5 Time evolution of the various ratios of BW, CD, and RS radii. The temporal evolution of the CD (blue) and RS (orange) radii normalised by the BW radius is displayed for 500 yr both for the simulation without (left) and with (right) turbulence. The plain curve corresponds to the average over all solid angles of the distribution of the ratio of the radii, and the coloured zones to their standard deviation. The hatched areas correspond to the radius ratios as observed in Warren et al. (2005).

Fig. 6 Time evolution of the standard deviation of the CD after spherical compensation. The values displayed here correspond to the standard deviation of the distribution of CD/RS radii multiplied by the average RS radius. The RTI at the CD, and therefore the mixing zone, shows faster growth in a turbulent ISM (orange) than in a non-turbulent one (blue).

3.5 Qualitative comparison to observations

Our simulations can also be qualitatively compared to SNR observations. While such comparisons are more limited as they are based on subjective factors, due to the lack of metrics, they can still provide information as to the global trends.

A direct observation of the results of our simulations, such as the one displayed in Fig. 1, already contains a multitude of information on the remnants and more precisely on their morphology. We mention above that the remnant stays almost spherical with the exception of the CD in the non-turbulent case. To this, we should add that, while RTI grows at the CD, most of it remains at a limited size. As a result, there is a gap between RTI spikes and the BW. On the other hand, in the presence of a turbulent ISM, the remnant is deformed as a whole; the RTI is less regular, and we can observe long deformed spikes that reach the BW and even push it. This results in protrusions of the BW, which do not exist in the non-turbulent case (see Appendix). These deformations (loss of spherical symmetry and protrusions) can be found in real systems such as Tycho. This could be seen as further proof of the importance of ISM perturbation on SNR evolution.

We should also mention the limits of our simulations. It is well known that the BW of an SNR is non-collisional, because the conditions of the plasma estimated at this interface suggest a mean-free path of the order of the radius of the SNR or even more. Our simulations are based on a hydrodynamic code FLASH4. As such, all shocks present here are collisional. This qualitative difference between code and observation might not be obvious here, but it constitutes a major limitation of our work even if the trends are the same (deformation of the BW, enhancement of RTI growth).

3.6 Comparison to observed velocities

In a recent analysis of Tycho’s emission, the radial profiles of the line-of-sight (LoS) velocity were measured based on the Doppler shift of the emission lines of the heavy elements of the SNR. From the centre of the remnant to its edge, the profile of the LoS velocity drops, approximately following a $\sqrt{\left(R-r\right)/\left(R+r\right)}$ profile (with R the radius of the external ejecta and r the position from the centre), and then plateaus near the edge of the SNR. While the first part can be understood and modelled by simple spherical shell-like emission, the plateau is controversial as it cannot be recreated by this simple model. It has been suggested that this plateau feature may be due to the interaction of the SNR with a dense medium.

By post-processing our simulations, we can study the effect of turbulence on the radial profile of the LoS velocity. Figure 7 shows the radial profiles of the different simulations (with or without turbulence) based on two different methods to calculate the density-averaged LoS velocity, that is, considering either the shocked shell or the shocked ejecta as the emitting medium. The shocked shell is almost equivalent to the results of a simplified model, where a uniform spherical shell emits light. However, we believe the shocked ejecta to be more relevant. As opposed to the shocked shell, this latter contains the heavy elements producing the line emission, whose Doppler shifts can be observed. The results shown in Fig. 7 use both idealised values and real values for the velocity and density. The idealised values suppose a constant radial velocity and a constant density. The real values use the results of the simulation.

As can be seen in the figure, idealised and real values bear similar results, with noisier profiles for the real values. This means the exact velocity and density distributions are close to their idealised value and that they have a low impact on the observed LoS velocity. On the other hand, the nature of the emitters has an important impact on the profile. For instance, in the case of idealised values with the shocked shell and without turbulence, the LoS velocity profile obtained corresponds to the results of a uniform spherical shell model, which shows that the morphology of our simulated shocked shell does not depart much from a spherical shell. However, this is inconsistent with the observations, as there is no plateau at the edge of the profile. On the other hand, the radial profile is completely changed when considering only the shocked ejecta. In particular, a plateau appears at the SNR edge in the turbulence case. This plateau is an indirect result of the RTI growth, which shapes the ejecta boundary. The regularity of this instability in the non-turbulent case leads to the quasi-linear drop in LoS velocity. In the case of turbulence, the loss of regularity leads to the superimposition of large and small spikes, and as a result, a plateau may be created.

From the profiles we obtained, we can conclude that ISM turbulence can lead to the creation of a plateau and that the LoS velocity observed corresponds to the ejecta. In this case, the beginning of the plateau almost coincides with the base of the bubble and its end with the tip of the RTI spikes. As such, the observation of this plateau might be a good measure of the RTI growth and resulting mixing zone at the CD.

Fig. 7 Line-of-sight (LoS) velocity 450 yr after SN explosion. The average density-weighted LoS velocity is displayed both for the nonturbulent (a, b) and turbulent (c, d) cases. The calculation of the LoS velocity uses either the shocked shell contained between BW and RS (a, c) or the shocked ejecta contained between CD and RS (b, d). Both idealised values (blue) and simulated values (orange) of the velocity and density are considered for the calculation.

4 Conclusions

In conclusion, we performed 3D hydrodynamic simulations of the evolution of a young SNR, initialised on conditions designed to correspond to Tycho. These simulations employ only a constant adiabatic index, contrary to the usual SNR simulation, which uses an effective adiabatic index to take into account the effect of CR acceleration. The main aim of these simulations is to explore the effect that turbulent-like density variations of the ISM would have on SNR dynamics.

We show that a turbulent ISM only has a small effect on the average position of BW, RS, and CD. However, this turbulence leads to a loss of spherical symmetry in the remnant and a higher RTI growth at the CD. This faster-growing RTI might explain the aspect ratio of BW:CD:RS measured in Warren et al. (2005), and also appears to be responsible for the plateau observed on the radial profiles of the LoS velocity. Therefore, measurement of this plateau is a good indication of the width of the mixing zone, which is linked to RTI growth.

Acknowledgements

G. Rigon’s work was performed as an international research fellow of the Japan Society for the Promotion of Science (postdoctoral fellowship for research in Japan (Standard)). This work is supported by Grants-in-aid from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan (20H01944). Simulations were performed on the XC50 supercomputer of the Center for Computational Astrophysics of the National Astronomical Observatory of Japan, under XC-trial and XC-B type proposals. Furthermore, the authors acknowledge Tomoaki Kasuga for his work on the LoS velocity of Tycho.

Appendix A Deriving the equations of turbulence

As the relation between equation (14) and the set of equations from (15) to (18) might not be obvious, we describe the details of this transition below. $$\delta \overline{\rho }\left(r\right)={\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}\sqrt{C{k}^{-\alpha }}\sin \left[2\pi \left(k\cdot r_L+\text{Φ}\left(k\right)\right)\right]}$$(A.1) $$={\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}\sqrt{C{k}^{-\alpha }}\left(e^{2i\pi \left(k\cdot r_L+\text{Φ}\left(k\right)\right)}-e^{-2i\pi \left(k\cdot r_L+\text{Φ}\left(k\right)\right)}\right)/2i}$$(A.2) $$={\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}\sqrt{C{k}^{-\alpha }}{e}^{2i\pi k\cdot r_L}\left(e^{2i\pi \text{Φ}\left(k\right)}-e^{-2i\pi \text{Φ}\left(-k\right)}\right)/2i}$$(A.3) $$={\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}{H}_0\left(k\right)e^{2i\pi k\cdot r_L}}.$$(A.4)

In this form, we can recognise a discreet Fourier transform with Fourier coefficients equal to $\sqrt{C{k}^{-\alpha }}\left(e^{2i\pi \text{Φ}\left(k\right)}-e^{-2i\pi \text{Φ}\left(-k\right)}\right)/2i$. An FFT-type algorithm can be used to solve it.

Now, let us consider the average density, 〈ρ〉, of a simulation cell with centre (x₀, y₀, z₀) and size (Δ_x, Δ_y, Δ_z). $$\langle \rho \rangle \left(x_0,y_0,z_0\right)={\displaystyle {\int }_{x_0-{\text{Δ}}_x/2}^{x_0+{\text{Δ}}_x/2}{\displaystyle {\int }_{y_0-{\text{Δ}}_y/2}^{y_0+{\text{Δ}}_y/2}{\displaystyle {\int }_{z_0-{\text{Δ}}_z/2}^{z_0+{\text{Δ}}_z/2}\frac{\rho \left(x,y,z\right)}{{\text{Δ}}_x{\text{Δ}}_y{\text{Δ}}_z}\text{d}x\text{\hspace{0.17em}d}y\text{\hspace{0.17em}d}z}}}$$(A.5) $$={\rho }_0\left(1+{\displaystyle {\int }_{x_0-{\text{Δ}}_x/2}^{x_0+{\text{Δ}}_x/2}{\displaystyle {\int }_{y_0-{\text{Δ}}_y/2}^{y_0+{\text{Δ}}_y/2}{\displaystyle {\int }_{z_0-{\text{Δ}}_z/2}^{z_0+{\text{Δ}}_z/2}\frac{\delta \overline{\rho }\left(x,y,z\right)}{{\text{Δ}}_x{\text{Δ}}_y{\text{Δ}}_z}}\text{d}x\text{\hspace{0.17em}d}y\text{\hspace{0.17em}d}z}}\right)$$(A.6) $$={\rho }_0\left(1+{\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}\frac{H_0\left(k\right)}{{\text{Δ}}_x{\text{Δ}}_y{\text{Δ}}_z}}{\displaystyle {\int }_{x_0-{\text{Δ}}_x/2}^{x_0+{\text{Δ}}_x/2}{\displaystyle {\int }_{y_0-{\text{Δ}}_y/2}^{y_0+{\text{Δ}}_y/2}{\displaystyle {\int }_{z_0-{\text{Δ}}_z/2}^{z_0+{\text{Δ}}_z/2}{e}^{2i\pi \left(k_{x}x/L_x+k_{y}y/L_y+k_{z}z/L_z\right)}\text{d}x\text{\hspace{0.17em}d}y\text{\hspace{0.17em}d}z}}}\right)$$(A.7) $$={\rho }_0(1+{\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}\frac{H_0\left(k\right)L_x{L}_y{L}_z}{{\left(2i\pi \right)}^3{\text{Δ}}_x{\text{Δ}}_y{\text{Δ}}_z{k}_x{k}_y{k}_z}}\left(e^{2i\pi k_x\left(x_0+{\text{Δ}}_x/2\right)/L_x}-e^{{}^{2i\pi k_x\left(x_0-{\text{Δ}}_x/2\right)/L_x}}\right)$$(A.8) $$\left(e^{2i\pi k_y\left(y_0+{\text{Δ}}_y/2\right)/L_y}-e^{2i\pi k_y\left(y_0-{\text{Δ}}_y/2\right)/L_y}\right)\left(e^{2i\pi k_z\left(z_0+{\text{Δ}}_z/2\right)/L_z}-e^{2i\pi k_z\left(z_0-{\text{Δ}}_z/2\right)/L_z}\right))$$(A.9) $$={\rho }_0(1+{\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}\frac{H_0\left(k\right)L_x{L}_y{L}_z}{{\left(2i\pi \right)}^3{\text{Δ}}_x{\text{Δ}}_y{\text{Δ}}_z{k}_x{k}_y{k}_z}{e}^{2i\pi k_x{x}_0/L_x}\left(2i\text{\hspace{0.17em}sin\hspace{0.17em}}2\pi k_x{\text{Δ}}_x/2L_x\right)}$$(A.10) $$e^{2i\pi k_y{y}_0/L_y}\left(2i\sin 2\pi k_y{\text{Δ}}_y/2L_y\right)e^{2i\pi k_z{z}_0/L_z}\left(2i\sin 2\pi k_z{\text{Δ}}_z/2L_z\right))$$(A.11) $$={\rho }_0\left(1+{\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}{H}_0\left(k\right)e^{2i\pi k\cdot r_L}\text{sinc}\left(\pi k_x{\text{Δ}}_x/L_x\right)\text{sinc}\left(\pi k_y{\text{Δ}}_y/L_y\right)\text{sinc}\left(\pi k_z{\text{Δ}}_z/L_z\right)}\right)$$(A.12) $$={\rho }_0\left(1+{\displaystyle \sum _{k\in {〚-K;K〛}^3;k\ne 0}{H}_0\left(k\right)H_1\left(k\right)e^{2i\pi k\cdot r_L}}\right).$$(A.13)

We retrieve here the formula used in this article to generate the turbulence. The H₁ part depends on the size of the simulation cell and therefore evolves with time.

Appendix B Additional plots

Here are some additional plots showing the effect of the turbulence on the BW:CD:RS radii ratios (figures B.1a, B.4a) as well as the RTI growth (figures B.1b, B.4b). The number that indexes the turbulence plot corresponds to a different set of random phases used to initialise the turbulence in the simulation.

In addition, the projection of the density for different simulations is displayed in figure B.5.

Fig. B.1 Effect of the turbulence random phase set on the results of the simulations. The simulations are based on the piston-like profile (the same initial profile as the main article). As can be seen, the specific set of random phases influences the global evolution of the SNR without changing the global trend (RTI growth and radius ratios).

Fig. B.2 Effect of the turbulence phase of a low-temperature ISM on the evolution of an SNR. Here the simulation uses the same parameters as the article: a piston-like profile for the SNR and the same sets of phases for the turbulence. Only the temperature of the ISM differs, here 10 K. The ISM temperature has a limited effect on the system.

Fig. B.3 Effect of the strength of the turbulence on the evolution of an SNR. Here the simulation uses the same parameters as in the article, except for the strength of the turbulence, which includes a lower-σ value (0.05) than the one used in the article (0.5). The strength of the turbulence is important for our conclusions.

Fig. B.4 Effect of the initial SNR profile on the simulation results. Contrary to the rest of the article, the simulations here are based on the DDTc profile taken from Warren et al. (2005). While the specific dynamic differs from our piston-like profile, the effect of the ISM turbulence on the RTI growth and on the radius is the same as in our piston model.

Fig. B.5 Density projection. Projection of the density of the SNR (ejecta and shocked ISM) according to the three spatial axes (columns) taken 450 yr after SNR explosion. Here, four different situations are represented (one per line). The ‘no turbulence’ case corresponds to the simulation with a uniform ISM. The ‘low turbulence’ case includes turbulence and uses the lower-σ value (0.05). This case does not differ significantly from the ‘no turbulence’ one. Both ‘turbulence 1’ and ‘turbulence 2’ employ a turbulent ISM with σ = 0.5, but each uses a different set of random phases. The BW appears wavy in these cases, which is closer to the morphology of the BW observed in Tycho.

References

All Figures

Fig. 1 Density slice taken 450 yr after SN explosion. The slices are performed on the YZ plane (x = 0) for both types of simulation: without turbulence (a), and with turbulence (b). The spatial coordinates and density values displayed here are taken in the inertial reference frame. Here, the main effect of the turbulence is the loss of the spherical symmetry of the SNR and a change in the instability wavelength and amplitude at CD.

In the text

	Fig. 2 Radial average profiles 450 yr after the SNR explosion. The radial profiles of the density, ejecta, and norm of the velocity are displayed in the inertial coordinates both with (plain curve) and without (dashed curve) turbulent ISM. These profiles are averaged over all angles of the octant.
In the text

	Fig. 3 Radial profiles for an arbitrary solid angle 450 yr after SNR explosion. The pressure and density profiles of the simulation with turbulence are displayed. The discontinuities in the pressure profile correspond to the BW and RS. The CD is defined as the point where the ejecta crosses 0.1.
In the text

Fig. 4 Spatial distribution of the SNR radius. Histogram of the RS (blue), CD (orange), BW (green) radius both without (a) and with (c) turbulence obtained 450 yr after SN explosion (simulations displayed in Fig. 1). To compensate for the spatial variation of the BW radius appearing in the turbulent case, the RS and CD radii are normalised using their respective BW radius for each solid angle both without (b) and with (d) turbulence.

In the text

Fig. 5 Time evolution of the various ratios of BW, CD, and RS radii. The temporal evolution of the CD (blue) and RS (orange) radii normalised by the BW radius is displayed for 500 yr both for the simulation without (left) and with (right) turbulence. The plain curve corresponds to the average over all solid angles of the distribution of the ratio of the radii, and the coloured zones to their standard deviation. The hatched areas correspond to the radius ratios as observed in Warren et al. (2005).

In the text

	Fig. 6 Time evolution of the standard deviation of the CD after spherical compensation. The values displayed here correspond to the standard deviation of the distribution of CD/RS radii multiplied by the average RS radius. The RTI at the CD, and therefore the mixing zone, shows faster growth in a turbulent ISM (orange) than in a non-turbulent one (blue).
In the text

Fig. 7 Line-of-sight (LoS) velocity 450 yr after SN explosion. The average density-weighted LoS velocity is displayed both for the nonturbulent (a, b) and turbulent (c, d) cases. The calculation of the LoS velocity uses either the shocked shell contained between BW and RS (a, c) or the shocked ejecta contained between CD and RS (b, d). Both idealised values (blue) and simulated values (orange) of the velocity and density are considered for the calculation.

In the text

	Fig. B.1 Effect of the turbulence random phase set on the results of the simulations. The simulations are based on the piston-like profile (the same initial profile as the main article). As can be seen, the specific set of random phases influences the global evolution of the SNR without changing the global trend (RTI growth and radius ratios).
In the text

	Fig. B.2 Effect of the turbulence phase of a low-temperature ISM on the evolution of an SNR. Here the simulation uses the same parameters as the article: a piston-like profile for the SNR and the same sets of phases for the turbulence. Only the temperature of the ISM differs, here 10 K. The ISM temperature has a limited effect on the system.
In the text

	Fig. B.3 Effect of the strength of the turbulence on the evolution of an SNR. Here the simulation uses the same parameters as in the article, except for the strength of the turbulence, which includes a lower-σ value (0.05) than the one used in the article (0.5). The strength of the turbulence is important for our conclusions.
In the text

	Fig. B.4 Effect of the initial SNR profile on the simulation results. Contrary to the rest of the article, the simulations here are based on the DDTc profile taken from Warren et al. (2005). While the specific dynamic differs from our piston-like profile, the effect of the ISM turbulence on the RTI growth and on the radius is the same as in our piston model.
In the text

Fig. B.5 Density projection. Projection of the density of the SNR (ejecta and shocked ISM) according to the three spatial axes (columns) taken 450 yr after SNR explosion. Here, four different situations are represented (one per line). The ‘no turbulence’ case corresponds to the simulation with a uniform ISM. The ‘low turbulence’ case includes turbulence and uses the lower-σ value (0.05). This case does not differ significantly from the ‘no turbulence’ one. Both ‘turbulence 1’ and ‘turbulence 2’ employ a turbulent ISM with σ = 0.5, but each uses a different set of random phases. The BW appears wavy in these cases, which is closer to the morphology of the BW observed in Tycho.

In the text