© ESO 2018

1 Introduction

The study of the structure of supernova remnants (SNRs) is a difficult task because SNRs experience various physical processes during their evolution. From the beginning of its spherical propagation, the blast wave (BW) created by the supernova (SN) pushes away the surrounding material, that is, the circumstellar medium and/or interstellar medium (ISM), and a dense shell forms just behind the BW front while a nearly empty but very hot bubble is created in the inner, central region, of the SNR. As the shell density is very large compared to the density of the surrounding medium and as the shell is strongly decelerated, such a configuration is unstable against hydrodynamic instabilities. For example, in the 1970s, Chevalier (1974) showed that Rayleigh–Taylor instabilities (RTI) might have played a role in the apparition of filaments in the outer region of the Crab Nebula. The existence of RTI is supported from observations of SNRs where complex and messy structures are in evidence.

In the 1980s another instability was proposed analytically by Vishniac (1983) and Ryu & Vishniac (1987) in the context of expanding spherical BW in SNRs. The former study was carried out for an isothermal SNR, that is, one in which the adiabatic index γ is γ = 1 (this case is called “radiative” by Vishniac – see below) while in the latter an adiabatic behavior, one in which γ satisfies γ > 1, was assumed by the authors. Nevertheless, the instability was shown to occur in both situations and it was suspected to explain the fragmentation of the front of the propagating BW. In the following, we will refer to this instability as the Vishniac instability (VI), although this mechanism is called “overstability” by Vishniac (1983).

The stability of BWs is therefore a topic that has received attention for several decade. In spite of the numerous efforts that have been already made, a topical issue remains because some areas need still to be clarified. This paper is the continuation of several steps achieved in our hierarchical study of the development of the VI and presents new numerical investigations. They take into account a different approach concerning the BW expansion in the radiative regime; radiative cooling is included in our model by energy losses, therefore the adiabatic index γ does not need to approach one in order to mimic radiative losses, contrary to the Vishniac (1983) theory.

Indeed, in his paper, Vishniac considers that the VI arises in radiative SNRs where both their structure (high density regions) and dynamics (strong deceleration) are altered by cooling effects. Actually, the authors do not include any cooling in their studies (Vishniac 1983; Ryu & Vishniac 1987). They used an adiabatic index γ going to one in order to get a strong compression in the shocked material. Such behavior arises when radiation escapes from the downstream flow (Draine & McKee 1993). Indeed, the compression given by the ratio of “density in the downstream flow” to “density in the upstream flow” in the vicinity of a high Mach number BW front is . It is obvious that by decreasing γ, the compression C increases. Nevertheless, we will demonstrate that the evolution of all quantities in the BW is not equivalent in both situations, either γ → 1 or really modeling radiative losses by a cooling function. Therefore, in the present work we are going to numerically revisit the VI in a more realistic context.

This paper is organized as follows. In Sect. 2, we give a general overview of the VI to remember the general mechanism and most important results. In Sect. 3, the HADES (Hydrodynamique Adaptée à la Description d’Ecoulements Supersoniques) two-dimensional (2D)-numerical code used for this work is presented together with the initial and boundary conditions implemented for the numerical simulations. Simulations of the VI are performed with no cooling in Sect. 4. We explain that even if they can look similar to previous work (Michaut et al. 2012), they are different and closer to the theoretical study (Ryu & Vishniac 1987). Section 5 deals with the evolution of an unperturbed SNR (no VI instability is triggered), which experiences cooling, and the effects of radiative losses on the SNR dynamics are discussed. Then, the coupling between the VI and the cooled SNR is studied in Sect. 6. The highly nonlinear behavior of the perturbation is described and analyzed. These results are completely new and reveal an unknown behavior of the VI Finally, the various situations (VI only, cooling only, and VI with cooling) are compared in Sect. 7 and our conclusion is also given.

2 Brief description of the Vishniac Instability

Let us examine the way the VI happens for a very thin shell. This corresponds to the initial work by Vishniac (1983). Due to the high thermal pressure in the hot central bubble of the SNR, the inner surface of the outer shell is pushed by a thermal force, F _th, which is always perpendicular to this surface (see Fig. 1). At the same time, while expanding, the shell accelerates the ISM that is initially at rest. This linear momentum transfer from the SNR to the ISM can be seen as a pressure, called “ram pressure”, exerting a deceleration force (“ram” force, F _ram) on the outer surface of the shell. The ram force is always colinear to the direction of expansion of the SNR, that is, it is always radial. If the SNR is perfectly spherical, these two forces (thermal and ram) are strictly colinear, and of opposite directions (their sum is of course not zero). However, if for any reason the SNR is no longer spherical due to any local perturbances (heterogeneities in the ISM, for example), these two forces are no longer colinear and their vector sum has a transverse component (see Fig. 1) that generates fluid motion (white arrows) along the shell from the extended front zones (A), called “hills” hereafter, to the receded shell zones (B), called “valleys” hereafter.

As the material flows from the hills to the valleys (we consider a perturbation with several modulations), the latter becomes more massive while the former loses some mass, so that the hills are more decelerated (due to their decreased inertia) than the valleys, the inertia of which have increased. This displacement of matter implies that the purely radial component of the velocity of a valley will be larger than the velocity (purely radial too) of a hill and, soon after a valley will overtake the average (over the angle) spherical position of the outer border of the shell, while the hill recedes with respect to this border. As a consequence, a valley (resp. a hill) will become a hill (resp. a valley) and oscillations take place in the shell. As mentioned in Sect. 1, this process was called “overstability” by Vishniac (1983) and according to him it might lead to the fragmentation of the shell provided that γ obeys γ < 1.2 (Vishniac 1983; Ryu & Vishniac 1987). Actually, other analytical approaches (Vishniac 1983; Ryu & Vishniac 1987; Kushnir et al. 2005; Sanz et al. 2016) show that for such low γ a domain of eigen modes l exists (see the definition of l in Sect. 3) for which a shell perturbation at the BW front should grow because the real part, s_r (γ, l), of the complex instability growth rate, s(γ, l), is positive, Re(s) > 0 (see further in the next sections and in Sanz et al. 2016).

The current state of our previous numerical investigations can be summarized as follows. First numerical works by Cavet et al. (2009, 2011), in both planar and cylindrical geometries, have fulfilled the conditions for the apparition of the VI depending on several initialization parameters. Nevertheless, the main conclusion is that in spite of an oscillating behavior, after one or two oscillations, that is, 10–30 kyr, the perturbation decays and no amplification is observed even for a low adiabatic index. Afterwards, Michaut et al. (2012) demonstrated the same behavior: triggering the instability, involving the overstability process and then vanishing much later, even in spherical geometry and using a low adiabatic index. We remind readers that all these calculations were performed with an adiabatic index γ = 5∕3 in the ISM and the SNR’s interior. However, in the thin shell, γ has been taken to vary in the range , in order to mimic the effect of radiative losses. This choice was made in order to approach the SNR expansion physics well enough. Using this strategy, all our previous numerical simulations have demonstrated that for low adiabatic indices (γ → 1) in the thin shell, the perturbation grows faster than for γ = 5∕3, at early times, as expected by theory.

Fig. 1 Sketch of the mechanism of the VI done in parallel-plane geometry illustrating a small part of the shell. The structure of the SNR is modeled by a hot bubble (left part) and a distorted thin dense shell (middle) that expands in the ISM (right) at the horizontal velocity V s. In this perturbed configuration, the thermal force, Fth, and the ram force, Fram, are not colinear, so that matter flows (white arrows) from the extended front zone (A), called “hill”, to the receded shell zone (B), called “valley”.

3 Methods and physical problem setup

In this paper, we perform axisymmetric numerical simulations of SNRs with the 2D radiation hydrodynamic eulerian code HADES (Nguyen 2011; Michaut et al. 2017). The equations of the model (Michaut et al. 2011; Busschaert et al. 2013) describing the evolution of SNRs are as follows: (1)

where ρ(r, t), v(r, t), p(r, t), I, and E(r, t) are respectively the fluid mass density, velocity, pressure, identity tensor, and total energy (thermal plus kinetic) density, and where Λ (r, t) is the cooling function detailed further. Using the closure relation (2)

where γ is the adiabatic coefficient, the system has only three independent variables ρ, v, and E.

The simulations are performed in a cartesian 2D mesh with the multi-processor version of HADES on 144 computer cores and the geometry used in the simulations is explained in Fig. 2. Any point P can be represented by its cartesian coordinates, x and y, or by its polar coordinates, r and θ, in the plane (x, y). We calculate the equation system (1) in a quarter of the disk (domain 1) taking into account directly both symmetries: an equatorial symmetry with respect to the y-axis to deal with half a disk (domains 1 + 2) and a 2π-rotation around the horizontal x-axis (domain 3) in order to reconstruct a full sphere. These symmetries are simulated by adequate boundary conditions and source terms (Michaut et al. 2017). The code simulates therefore a three-dimensional (3D) object with axisymmetry about the x-axis.

All our simulations are made in two steps: the first creates the Sedov–Taylor BW, the second triggers the VI and/or cools the BW. The gas composing the ISM is ionized hydrogen of mass density ρ_ISM ≈ 10⁻²⁰ kg m⁻³, which corresponds to about 10⁷ part m⁻³ ≈ 10 part cm⁻³ at a pressure of P_ISM ≈ 10⁻¹² Pa as T_ISM ≈ 10⁴ K.

Fig. 2 Geometry of the configuration: the numerical domain 1, a quarter of disk (), undergoes an equatorial symmetry with respect to y () defining domain 2, and a rotational symmetry around the x-axis () giving domain 3. It is representative of a 3D axi-symmetric object.

3.1 First step: blast wave evolving until 3 kyr

We initialize our simulations with an instantaneous thermal energy deposit, E_SN, in a small spherical central region in order to obtain a SNR evolving in the Sedov–Taylor regime. To preserve sphericity in the cartesian mesh, it is four times refined with d x =dy = 0.4 × 10¹⁴ m for an initial sphere of radius R_dep = 4 × 10¹⁵ m. This initial energy is E_SN = 10⁴⁴ J, which is a value we canobserve in SNRs (Utrobin & Chugai 2008, 2013; Fink et al. 2014). For these simulations, we solve the system (1) without cooling, that is, Λ = 0, because the SNR is supposed to be too optically thick to experience radiation losses (Woltjer 1972).

3.2 Second step: blast wave undergoing perturbations and/or cooling

By projecting the data resulting from the first step, at t₀ ≈ 3 kyr, into a coarse grid composed of 5400 × 5400 squared mesh cells, the simulations are initialized. The numerical domain is about 8.6 × 10¹⁷ m ≈ 28 pc in each direction x and y. The spatial step is therefore dx = dy = 1.6 × 10¹⁴ m and the resolution of these simulations is two times better than those from Cavet et al. (2009, 2011) and Michaut et al. (2012). An investigation of the mesh convergence is given in Appendix A.

To study the VI as presented in Sects. 4 and 6, a spatial perturbation is introduced at t₀, in order to trigger an instability. This time around 3 kyr is of the order of magnitude for which deviations from sphericity in the shape of the SNR might occur (for instance, due to density heterogeneities in the surrounding medium; Woltjer 1972). The perturbation is obtained by replacing the fluid quantities at a given coordinate (r,θ) with the corresponding quantities at the position (r + δr, θ), where δr is the spatial perturbation of r now depending on θ: (3)

with l the eigenmode and A the amplitude of the perturbation. In this way, any quantity in the flow is replaced with a new value that corresponds to the perturbed quantity given by (4)

where q stands for any quantity in the flow, density ρ, velocity v, pressure p, and energy E, displayed in system (1). As a consequence, the whole flow, not only the outer shell, experiences the perturbation like in the theoretical studies by Ryu & Vishniac (1987), Kushnir et al. (2005), and Sanz et al. (2016).

In this study, we have chosen l = 24 and A = 2.5 × 10¹⁵ m. These parameter values are taken in keeping with our previous numerical study presented in Michaut et al. (2012) and allowing a growth of the instability according to Ryu & Vishniac (1987). In Sect. 4, where γ = 1.1, A represents approximately 2% of the SNR radius at t₀ (R(t₀) ≈ 1.2 × 10¹⁷ m) and 8% of the initial wavelength (λ = 2πR(t₀)∕l = π × 10¹⁶ m). In Sect. 6, where γ = 5∕3, A represents approximately 1.4% of the SNR radius at t₀ (R(t₀) ≈ 1.8 × 10¹⁷ m) and 5% of the initial wavelength (λ = 2πR(t₀)∕l ≈ 4.7 × 10¹⁶ m).

In the further simulations, including the cooling that is linear in the density ρ and without temperature dependence (see details in Sects. 5 and 6), the cooling function Λ in the system (1) is simply turned on. Including radiative cooling during the evolution of SNRs is expected to model the radiative phase (McKee & Ostriker 1977; Cioffi et al. 1988). However, the analytical form of the cooling function is not obvious and several authors have proposed different expressions (Woltjer 1972; Raymond et al. 1976; Bertschinger 1986; Lequeux 2005). After several tests, we have chosen to use a cooling function, Λ, of the form , where Λ ₀ is a constant. Radiation losses of the form Λ ∝ ρ² would have been more relevant to account for line cooling but the simple form Λ ∝ ρ allows for faster numerical runs (as nonlinearities are less stiff than for Λ ∝ ρ²) although each run takes several days on a massively parallel device (see below). Actually, the ρ² case has been studied too but in less detail (see Sect. 5).

An optimal value of Λ₀ is therefore required. Fortunately, it is straightforward to fulfill this condition comparing the total energy, E_cool (t), lost by the SNR due to cooling until time t with the energy, E_SN, released by the SN. At a time t and in the total SNR volume, we have

because of the linear form of Λ = Λ₀ρ, with M_SNR the total mass of the SNR. When the SNR is tens of thousands of years, all energy is dissipated and Eq. (7) provides the order of magnitude of the cooling parameter, Λ₀ = 0.1 W kg⁻¹ (more details in Sect. 5).

Numerical simulations are run correctly if the local temperature does not drop below the ISM temperature (T ≥ T_ISM ≈ 10⁴ K), which we call the cut-off temperature. Actually, the drop in is due to both the cooling function Λ and the SNR expansion. Once T has decreased down to T_ISM, we keep T = T_ISM by neglecting the role of the adiabatic expansion. This approximation consists in a kind of heating but the involved amount of energy is small and we believe it does not change significantly the SNR evolution.

We have to bear in mind that Vishniac (1983) and Ryu & Vishniac (1987) did not include a cooling function in their theory, but they chose a low adiabatic index γ → 1 with the idea of mimicking the cooling. Indeed, the energy losses by radiation, that is, by cooling function, lead to the matter compression magnification that is almost the same effect as in decreasing γ, as explained above (Sect. 1).

4 Vishniac instability in SNR without cooling

In order to compare the evolution of a perturbation during the expansion of a SNR closer to the previous analytical studies by Vishniac (1983) and Ryu & Vishniac (1987), we have decided to set the adiabatic index to γ = 1.1 everywhere with no cooling function. Indeed, these analytical studies predict a positive value for the real part, Re (s) = s_r, of the perturbation growth rate,s(γ, l), for such a value of γ: that is, the perturbation is unstable. This approach is quite different from our previous simulations (Michaut et al. 2012) in which γ = 1.1 only in the thin shell and leads to an evanescent deformation after a long period (~ 30 kyr). In addition, these new simulations will be compared with the second part of this work including a cooling function.

We introduce the perturbation as explained above in Sect. 3, Eq. (4). We let it evolve during t ≈ 200 kyr in order to see the transition from the linear regime to the nonlinear one. As the amplitude A of the sinusoidal disturbance is very smallcompared to R(t₀), we have decided to show in Fig. 3, and for all further density maps, a zoom of the thin shell along the axis θ = π∕4. As a consequence, even if the mode number is l = 24, that is, six oscillations by a quarter of the disk, we can see only two of them.

Figure 3 shows the evolution of the density in the outer shell at four different moments and exhibits the overstability mechanism by the exchanges between hill and valley positions (the radius r increases in the four panels). Indeed, in Fig. 3a, the initial hill top is located at about (x, y) ≈ (2.9, 2.9), whereas the two valley bottoms are at (x, y) ≈ (2.45, 3.15) and mutually at ≈ (3.15, 2.45). Using a multi-layer colorbar expressed in log-scale, the stratification and the compression of the shell is highlighted. A little later at 4 kyr, Fig. 3b shows the tangential component (perpendicular to r, i.e. the shock motion direction) of the fluid velocity.

The complexity of the tangential-velocity structure benefits the VI by producing a net transport of mass along and at the edge of the outer shell as explained in Sect. 2. Consequently, the fluxes of matter create an overdensity in the valleys located around (2.65, 3.50) and around (3.50, 2.65) in Fig. 3b.

After 4 kyr the two initial valleys mentioned in Fig. 3a are accelerated and move forward while the initial hill recedes with respect to the unperturbed case as explained in Sect. 1. The first inversion of the hills and valleys (not shown in Fig. 3) happens before 10 kyr. Consequently, the positions of hills and valleys have been exchanged (see Fig. 3c) and the mechanism can restart. Indeed, a second inversion takes place at around 20 kyr.

Deeper in the SNR, the transverse movements have opposite directions (see Fig. 3b) due to the fact that the pressuregradient in this region is no longer radial as in the purely spherical case. We remind readers that ρ, v, and p have been perturbed in the whole configuration and the pressure gradient is no longer isotropic. At t ≈ 53 kyr (see Fig. 3d), the perturbation is still present, and the oscillation time has increased substantially. The SNR still has transverse velocities at the shock front and the perturbation is still well visible and clearly present. The behavior for a uniform γ = 1.1 is therefore new because at t ≥ 50 kyr, the shell perturbation as well as transverse velocities still occur contrary to our previous numerical results (Cavet et al. 2011; Michaut et al. 2012).

One way to measure the evolution of the instability is to calculate the mass variation, Δ M_i, for each type of region, initial valley or hill. A sketch of these regions is shown in Fig. 4 for l = 24, that is, six wavelengths per quarter. We calculate the mass M of each angular sector corresponding respectively to an initial hill or valley. Obviously, each sector covers half a wavelength (λ∕2) of the perturbation. So, we have 2l regions labeled by i = 1, …, 2l in the whole sphere, that is, l∕2 regions in the domain 1 (Fig. 2). With our conventions the initial hill regions are marked with an odd number whereas the initial valley regions have an even number. The mass expression, for each region i, is (8)

where R_SNR(θ, t) is the shock front location of the SNR.

Then, we compare it to the mass M_reg,i that these same regions should have without perturbation: (9)

We have to bear in mind that the masses M_reg,i are not the same for all regions due to the spherical geometry.

The mass variations ΔM_i are then defined by (10)

Figure 5 shows the quantity in percent as a function of time, where the index i varies from 5 to 9. We see that the amplitude of Δ M_i is all the greater given that i increases. The inversions are located where ΔM_i = 0, we cannot see the first inversion occurring before 10 kyr, and the second and the third inversions happen at respectively t ≈ 20 kyr and t ≈ 80 kyr, as shown in Fig. 5. We note that the ΔM_is change their sign between two inversions and they are positive in the range [10, 20] kyr and after 80 kyr for the valleys (i =6, 8) since they are negative in the range [20, 80] kyr. As expected, the opposite behavior takes place for the hills (i = 5, 7, 9). These oscillations exhibit the VI mechanism.

As a result, SNR’s of adiabatic index γ = 1.1 are much more unstable than for γ = 5∕3 in agreement with Vishniac (1983) and Ryu & Vishniac (1987). The oscillations of the hills and the valleys are well observed, even over long periods, which was not the case in the previous study (Michaut et al. 2012). Our work is closer to the analytical assumptions, in which the adiabatic index is uniform. However, even if we see the oscillations for a very long time as predicted, we do not observe a growth of the perturbation leading to a fragmentation of the SNR.

Fig. 3 Evolution in time and in space of mass density of shock front region along a radius; density is given by color table in 10−20 kg m−3; panel a: at t = t0 = 3 kyr, panel b: at t ≈ 4 kyr and highlighting transverse velocity with arrows, panel c: at t ≈ 18 kyr, panel d: at t ≈ 53 kyr.

5 Unperturbed SNR with cooling

In this section, we study the evolution of a SNR experiencing radiative losses (cooling) but not submitted to any disturbances at all. Cooling is expected to produce a compression of matter behind the BW front and in this section we are going to try to understand the effect of cooling on the structure of the SNR.

As explained in Sect. 3, we have studied the effect of different cooling functions, Λ =Λ₀, Λ =Λ₀ρ and Λ =Λ₀ρ². As may be guessed, the first cooling law does not lead to representative structures (absence of high compression), because the initial distribution of matter is not accounted for in this uniform, oversimplified, form of Λ. The simulations have shown that the second and the third cooling laws lead both to dense structures (see further), and in addition, they give comparable results. As a consequence, we have selected the second law (Λ ∝ ρ) in this publication because the simple linear form helps in the estimation of the value of Λ ₀ in an analytical way as mentioned in Sect. 3.

Indeed, for a radius R(t) given by the Sedov–Taylor law (Sedov 1946, 1959; Taylor 1950), , and using Eq. (7) we can rewrite the amount of radiated energy, E_cool(t), at time t as

This formula allows us to find the order of magnitude of the value of the cooling constant Λ ₀ once the relevant time for which the dynamics of the SNR is altered by the cooling. Following Cioffi et al. (1988), we consider that a substantial fraction of E_SN has been radiated away at t ≈ 1 00 000 yr, and setting E_cool ≈ E_SN, we obtain Λ₀ ≳ 0.01 W kg⁻¹. Actually, we have tested this value in numerical simulations and it came out that the compression of matter due to cooling was too slow and not sufficiently effective in contradiction with the observations of dense filaments in SNRs. For these reasons, the simulations presented in this work have been performed for several higher values of Λ ₀ = W kg⁻¹.

We show in Fig. 6 the density, temperature, pressure, and velocity profiles respectively at t = 33 kyr for these values of Λ₀, for γ = 5∕3. The larger the parameter Λ₀, the thinner and denser the shell (Fig. 6a). As expected, the highest compression (C ≈ 65) is obtained for Λ ₀ = 0.1 W⋅kg⁻¹ whereas the formula C = (γ + 1)∕(γ − 1), relevant for no cooling, gives C = 4 (see Λ₀ = 0 in Fig. 6a). We also notice that the higher the parameter Λ ₀, the more the SNR expansion is slowed (see Fig. 6). The hydrodynamic profiles are deeply changed in comparison with the Sedov–Taylor profiles (Λ₀ = 0).

Due to the cooling, three regions can be identified in the SNR, from the center to the outer border. This decomposition will be useful to understand the evolution of a SNR experiencing both cooling and a perturbation (see Sect. 6).

• 1

  The first region corresponds to a central inner bubble where the temperature remains very high (see Fig. 6b) due to the low mass density and, therefore, little cooling. The pressure is almost uniform as seen in Fig. 6c butat the outer edge of this zone the pressure (and also the temperature) drops drastically (the profiles become almost vertical) to a minimal value. This edge, which is not the shock front, is located at about r ≈ 12 pc for Λ₀ = 0.1 and 0.075 W kg⁻¹, r ≈ 12.5 pc for Λ₀ = 0.05 W kg⁻¹, and r ≈ 13.5 pc for Λ₀ = 0.03 W kg⁻¹ in Fig. 6c. In opposition to the Sedov–Taylor profile, the pressure gradient is always negative in this central region and the material is strongly accelerated outwardly. This process generates the velocity peak observed in Fig. 6d and the maximum velocity is no longer located at the shock front as in the Sedov–Taylor case.

• 2

  The second region corresponds to the shell located in between the outer edge of the first region (hot bubble) and the front of the BW as evidenced by the mass density peak. In this region, the SNR is cold because the temperature has fallen to the lowest value T_ISM due to the cooling, and for each value of Λ₀ the pressure gradient is positive as can be clearly seen in Fig. 6c.

• 3

  Finally, the region of the BW front corresponds to the third zone where a pressure peak arises together with the mass density peak. In this area, the temperature can be high (as observed in Fig. 6b) because the gas has been shocked recently and has not yet had enough time to cool significantly.

As a result, radiative losses strongly affect the structure of the SNR. For instance, in the Sedov–Taylor case, pressure increases monotonically from the center (r = 0) to the shock front R(t), whereas with cooling, it decreases from r = 0 to a minimum value located roughly at r = 90% R(t), and then the pressure grows to its maximum value at the shock front position R(t). The cooled material (second region) moves faster than the gas right behind the shock front and, therefore, this material catches up with the front, generates a huge compression, and forms a very thin and dense shell at the BW front, as expected. The radius of the dense shell is r ≈ 13.1 pc with C ≈ 65 for Λ ₀ = 0.1 W kg⁻¹ (Fig. 6a).

Moreover, Fig. 6d shows that the magnitude of the velocity grows with increasing values of Λ ₀. This property arises because in the second region the pressure gradient increases with Λ ₀. Such a structure is representative of the radiative stage of a SNR (compare, for instance, Fig. 6c to Fig. 3 in Cioffi et al. 1988) and equivalent to the description by Chevalier (1974) and Falle (1975, 1981). In these papers, the authors study the transition from the Sedov–Taylor stage to the radiative regime and the thin shell formation. They have a more realistic cooling function than us and let the SNR evolve from the adiabatic phase with cooling. They show that the radiative losses can generate additional shocks inside the SNR and drive a phenomenon called catastrophic cooling. In this paper, we have not examined the very detailed inner structure of the SNR, however, it is clear that our numerical simulations exhibit also inner structures that evolve with different pressures and velocities, and, as a consequence, collide with each other. These collisions lead to local growth of temperature, but since the cooling function we have chosen does not depend upon T, we do not really reproduce catastrophic cooling as originally stated by Falle (1975, 1981).

Energy losses produce an SNR that is much more decelerated than in the Sedov–Taylor case. Figure 7 shows the evolution of the radius, R(t), of the SNR using Λ₀ = 0.1 W kg⁻¹.

Definingthe exponent n such that R(t) ∝ tⁿ, the expansion initially follows the Sedov–Taylor law (n = 2∕5 = 0.4), very accurately from t ≈ 3 kyr (log t ≈ 0) to t ≈ 10 kyr (log t ≈ 0.5) as illustrated by the straight dashed line of slope 2/5 in the log–log plot. Then, a stronger deceleration of the shock front is visible from t ≈ 25 kyr (log t ≈ 0.9) and for t ≥ 35 kyr (log t ≥ 1.05) the deceleration rate n becomes constant and is n = 0.3. This is shown by the asymptotic dotted straight line with slope 0.3 in Fig. 7. The stabilization to the value n = 0.3 is in agreement with the numerical simulations by Cioffi et al. (1988), in which the cooling function comes from atomicphysics calculations of energy losses in the ISM, whereas we have found this law for both Λ ∝ρ and Λ ∝ρ².

The energy E_cool radiated by the SNR with Λ ₀ = 0.1 W⋅kg⁻¹ is exhibited in Fig. 8 for E_SN = 10⁴⁴ J. The figure shows that 90% of the initialenergy E_SN is radiated away over the whole evolution over 54 kyr. Moreover, about half of the total radiated energy has been lost over the first 15 kyr. At that moment, the coldest region (second region defined before with T =T_ISM) is already formed, and therefore this cold gas does not radiate anymore. The inner hot bubble (first region) still contributes to the radiative losses, and the matter right behind the shock front (third region) as well, so that the radiated energy continues to grow smoothly.

Fig. 4 Illustration of regions in which mass variation is calculated during the “overstability” mechanism; panel a: corresponding regions in the unperturbed case, panel b: regions in the perturbed case. Each region, numbered by , should be understood as a volume resulting from a 2π angular rotation around the x-axis.

Fig. 5 Mass variation ΔMi (%) by region depending on time. The evolution shows oscillations typical of the VI The regions are numbered by increasing index, i, with respect to the polar angle θ (see the main text). The amplitude of ΔMi is approximately constant with time for each region. We notice that the amplitude of Δ Mi grows with i.

Fig. 6 Comparison of profiles with γ = 5∕3 in the SNR at t ≈ 33 kyr for various values of Λ0 with the same color code, W kg−1; panel a: density, panel b: temperature, panel c: pressure, and panel d: velocity profiles. The case Λ 0 = 0 corresponds to the Sedov–Taylor case.

6 Vishniac instability in SNR with cooling

In this section, we present simulations in which we study both the cooling and the VI for γ = 5∕3 as in Sect. 5. We set the perturbation as defined in Sect. 3 and we use the previously studied cooling function, Λ =Λ₀ρ, with Λ ₀ = 0.1 W kg⁻¹ (see Sect. 5).

Simulations show that the perturbation undergoes attenuation from the beginning to t ≈ 38 kyr (see Fig. 9a). As γ = 5∕3, the amplitude of the perturbation decreases quickly at the beginning of the computation, when the cooling has not yet produced a strong effect on the SNR.Nevertheless the VI generates overdensities and structures behind the valleys, even if the oscillations have a decreasing amplitude. At t≈ 38 kyr, the thin dense shell is formed, and the modulation of the shock front has nearly totally vanished: the SNR is quasi-spherical.

According to the second zone as described in Sect. 5, the perturbed matter inside the SNR is cooled and accelerated and it reaches the shock front and reaccelerates the dense and thin shell. The traces of the early alterations in the SNR structure by the VI (Fig. 9a) are then transported from inside the SNR to the BW front and they produce a new perturbation of the shell, which can be seen in Figs. 9b–d.

In Fig. 9a, the initial wavelength λ = 2πR∕24 has been plotted and the valleys located in A and B (these two points are followed during the expansion in Fig. 9) begin to grow and become new bumps. In parallel, in Fig. 9b, the formation of an overdensity is observed at point C (this point is also followed during expansion). This overdensity is reminiscent of an older structure formed earlier (see Fig. 9a) and in Figs. 9c and d, C becomes also a bump. As a result, the inner material that reaches the BW front triggers a new perturbation with λ′ = λ∕2, i.e. l′ = 2l. The amplitude of the modulation then grows, as shown in Figs. 9c and d, and the SNR experiences the VI This is exhibited in Fig. 10 where points A, B, and C are also drawn and where the transverse mass flux ρv_T is mapped (v_T is the velocity component perpendicular to the radial direction). This flux of mass is evidence of matter transport from hills to valleys.

Figure 11 presents the mass variations for the angular sectors i = 5 to i = 9, defined in the same way as in Sect. 4. As explained previously, one can see the mass variation oscillations from t ≈ t₀ ≈ 3 kyr to t ≈ 50 kyr. Similar to the case with no cooling for γ = 5∕3 (see the end of Sect. 4), the VI is attenuated, the oscillations are not strong enough to change the excess of mass of the regions i = 6 and i = 8 (i.e., lack of mass in the regions i = 5, i = 7, and i = 9) into a lack of mass. After t ≈ 50 kyr, the ΔM_is decrease almost linearly with time and the curves intersect each other at t ≈ 80 kyr as they cross the horizontal straight line ΔM = 0. Later, the ΔM_is re-increase a little but this behavior is not representative of the actual configuration because a new perturbation of mode l′ =2 × l = 48 begins to grow at t ≈ 40−50 kyr, as explained before, and the angular sectors are too wide to account for the flow alterations due the perturbation with a two – times – shorter wavelength, λ′ = λ∕2.

This new growing perturbation is highlighted in Fig. 12, which shows the mass variation Δ M_i calculated for sectors with angular extension two times smaller than in Fig. 11. For instance, in Fig. 4b the initial sector i = 7 is split into two new subsectors, i = 7 and i′ = 7. It follows that instead of being divided into two angular sectors, each initial wavelength λ is now decomposed into four sectors of identical angular extensions. It can be seen that from t ≈ 3 kyr to t ≈ 30 kyr for l = 24, the mass variation in a sector depends on its location: an excess (resp. a lack) of mass is obtained for the sectors i′ =6 and i′ =8 (resp. i′ = 5 and i′ = 7), while for i = 5, 6, 7, and 8, we have ΔM ≈ 0 over this range of time. In contrast to Fig. 11 where the magnitude of the Δ M_is decreases after t ≈ 30 kyr, the magnitude of the ΔM_is grow until t ≈ 60 kyr in Fig. 12. Indeed, later the ΔM_is oscillate and these variations are evidence that the VI has taken place in the SNR modulated by a new disturbance of mode number l′ = 2 × l = 48. Moreover, it should be pointed out that since for i′ = 5 and 7 (resp. i′ = 6 and 8), we have Δ M < 0 (resp. Δ M > 0) for t < 30 kyr, the sectors i′ = 5 and 7 (resp. i′ = 6 and 8) correspond initially to valleys (resp. hills), but for t > 60 kyr the evolutions of these four structures become very similar (the four curves are rather close in Fig. 12) because they have given rise to the new overdensities in Figs. 9a–c. On the other hand, each of the sectors i = 5, 6, 7, and 14 is located between two overdensities and these sectors give birth to the new valleys. Finally, it is seen in Fig. 12 that from t ≈ 60 kyr to t ≈ 140 kyr, the amplitude of the oscillations increases from ΔM ≈ 30% to ΔM ≈ 40%. This is the numerical proof that the overstability mechanism as predicted by Vishniac (1983) takes place at the front of the BW in SNRs experiencing cooling. As a consequence, although a BW front is stable for a SNR with γ = 5∕3 (an initial perturbation vanishes with time and the VI does not develop), adding a cooling process in the SNR makes the BW front unstable against the VI.

Compared to the behavior obtained for γ = 1.1 (see Sect. 4), the instability is stronger for γ = 5∕3 with cooling. This is evidenced in Fig. 9d where the BW front deformation is more pronounced at t ≈ 93 kyr than at t ≈ 58 kyr (see Fig. 9c). This is also clearly seen from the comparison between Figs. 10a and b. Indeed, numerical simulations show that from t ≈ 58 kyr to t ≈ 93 kyr, the perturbation of the outer radius of theSNR grows from δR ≈ 0.6 × 10¹⁶ m to δR ≈ 1.2 × 10¹⁶ m, and assuming a time variation given by t^s_r, the value of the growth rate is 1.4 < s_r < 1.5.

Fig. 7 Variation of the radius of the SNR with respect to the time (full curve), in a log–log plot highlighting the transition from the Sedov–Taylor phase (dashed line,n = 0.4) to the asymptotic radiative regime (dotted line, n = 0.3); is expressed in pc and t in 1011 s.

Fig. 8 Energy (in 1042 J) radiated by the SNR as a function of time t (in kyr).

Fig. 9 Evolution in time and in space of mass density of shock front region undergoing cooling; density is given by color table in 10−20 kg m−3; panel a: at t ≈ 38 kyr, panel b: at t ≈ 43 kyr, panel c: at t ≈ 58 kyr, panel d: at t ≈ 93 kyr.

Fig. 10 Transverse mass flux (ρvT) map; panel a; at t ≈ 58 kyr, panel b: at t ≈ 93 kyr. We clearly see transverse fluxes in the dense shell, from the hills to the valleys, characteristic of the VI mechanism. They are still active at a very advanced age.

7 Conclusion

We have shown that the evolution of a spatial perturbation in a SNR is very different according to whether the SNR experiences cooling or not. First, for an adiabatic (i.e., not cooled) SNR, which is the situation studied theoretically by Vishniac (1983) and Ryu & Vishniac (1987), we have shown that for an adiabatic index γ = 1.1, the oscillating perturbation (overstability mechanism) is observed throughout the simulation time (t_final ≈ 200 kyr), as predicted bytheory. As far as we know, this is the first time that the evolution of the VI in SNRs is computed over such a long period.

In a second step, cooling was taken into account in the simulations in order to demonstrate the formation of a thin and very dense shell (radiative regime of the SNR). In addition, we have considered a SNR with γ = 5∕3 which, in opposition to γ = 1.1, is predicted to be stable against VI, and the influence of the cooling on the instability was investigated. In the first part of the evolution, the radius R(t) of the blast wave (BW) follows the Sedov–Taylor (ST) law, R ∝ t^2∕5, but for t ≥ 35 kyr the expansion is more decelerated, with R ∝ t^3∕10 because the energy radiated by the SNR, E_cool, becomes comparable to the energy released by the supernova, E_SN. During the ST stage of the BW, the perturbation decreases with time and the shock front recovers a smooth spherical shape. However, after the transient regime from R ∝ t^2∕5 to R ∝ t^3∕10, the alterations produced in the interior of the SNR by the outer initial perturbation with mode l reappear as a new growing, and therefore unstable, disturbance at the surface of the SNR with a mode l′ satisfying l′ = 2 × l. Cooling provokes therefore the VI and makes that the VI grows again at a later stage.

As a result, the evolution of a perturbed cooling SNR of adiabatic index γ = 5∕3 can be divided into threestages:

• In the first stage, radiative losses do not deeply modify the SNR structure, but the VI with the overstability mechanism is evidenced. However, as γ is large, the perturbation vanishes with time on the outer edge of the SNR.

• In the second stage, the material in the interior of the SNR, perturbed during the first step, is accelerated outwardly (see Sect. 5). When it reaches the outer dense thin shell it generates a new perturbation at the shock front with a doubled mode number.

• In the third stage, the new perturbation is seen to grow. The shock front undergoes large deformations and the SNR experiences again the VI.

The development of a doubled mode (l′ with l′ = 2 × l) is a new phenomenon which, according to us, is evidenced here for the first time. Furthermore, and strictly speaking, a BW is unstable against VI provided it evolves, first, adiabatically, and, second, γ satisfies γ ≤ 1.2 (Ryu & Vishniac 1987; Sanz et al. 2016). Nevertheless, the present work shows that cooling enhances this instability and, as a consequence, the SNRs we have studied undergo VI, although we have taken γ = 5∕3 and although they radiate energy (departure from adiabaticity). This result is also new and it suggests that any cooling SNR with γ < 5∕3 will be unstable against VI We note that for γ = 5∕3, the strong shock approximation provides an “effective” value of gamma, γ_eff, from the equality C = (γ_eff + 1)∕(γ_eff − 1), and with C ≈ 65, we get γ_eff ≈ 1.03, which is very low. A limitation to the comparison between analytical theoretical results and the simulations presented in this paper is dueto the form of the perturbation handled numerically. Theoretically, the spatial perturbation of the radius δR involves the spherical harmonics . As we have a 2D code, we cannot simulate perturbations depending on ϕ, and we have chosen to imprint a sinusoidal initial perturbation. However, it would be possible to work at least with spherical harmonics with m = 0. This particular case is thought to be representative of any Y _l,m as the perturbation growth rate does not depend upon the parameter m. Furthermore, the cooling function Λ we have used is highly idealized. It would be instructive to take a dependence upon mass density ρ and temperature T in Λ that matches an actual radiative cooling process. In our opinion, this would not change qualitatively our results but such more relevant simulations could be interesting in order to confirm the apparition of the VI in more realistic situations.

Fig. 11 Massvariation ΔMi (%) by region for the mode l = 24 depending on time, as explained in Sect. 4. The evolution shows oscillations typical of the VI mechanism from t ≈ 3 kyr to t ≈ 50 kyr.

Fig. 12 Massvariation ΔMi (%) by region for the mode l′ = 2 × l = 48 depending on time. From t ≈ 50 kyr to t ≈ 140 kyr, the regions corresponding to the mode l′ = 48 present the distinctive oscillations of ΔMi, confirming the development of the VI with a mode l′ twice the initial mode l.

Acknowledgements

This work was granted access to the HPC resources of MesoPSL financed by the Région île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the program Investissements d’Avenir supervised by the Agence Nationale pour la Recherche (ANR). J.S. is grateful to Commissariat à l’Énergie Atomique et aux Énergies Alternatives (CEA) and to Paris Observatory for their support for his stays in France, and S.E.B. thanks the Universidad Politecnica de Madrid for his visits to Spain.

Appendix A: Mesh resolution convergence

In Fig. A.1, we have plotted the mass variation Δ M_i as a function of time for five numerical resolutions for the angular sectors 8 and 9 shown in Fig. 5. The full blue lines corresponding to the conditions used for Fig. 5 are compared with lower and higher resolutions. Although they are symmetric, the dotted curves (at very low resolution) are obviously not correct because they oscillate about Δ M = 10% instead of ΔM = 0.

Fig. A.1 Mass variation ΔMi (%) of the two angular sectors i = 8 and i = 9 already shown in Fig. 5, for five different resolutions. The dotted curves correspond to a low resolution simulation with 2700 × 2700 mesh cells (dx = dy = 3.2 × 1014 m). The full cyan lines are obtained for the resolution used in Fig. 5 (high resolution with 54002 mesh cells and dx = dy = 1.6 × 1014 m). The other curves (full purple, dashed red, and dashed black) have been plotted for both intermediate (respectively for 36002, 43202 cells) and higher resolutions (72002 cells).

With such a large-cell mesh, the computing accuracy is not good enough. Consequently, the volume of the angular sectors is not accurately defined and their mass cannot be precisely evaluated. In contrast, the departures of the dashed lines compared to the full lines are about 5% and each of the pairs (full curves on the one hand, and dashed curves on the other hand) exhibits symmetry with respect to the horizontal line ΔM = 0. The calculation accuracy is high enough to satisfy mass preservation. In addition, the eight curves intersect at t ≈ 80 kyr as obtained in Fig. 5. It is therefore clear that the meshes we have used for the simulations presented in this work are converged.

References

All Figures

Fig. 1 Sketch of the mechanism of the VI done in parallel-plane geometry illustrating a small part of the shell. The structure of the SNR is modeled by a hot bubble (left part) and a distorted thin dense shell (middle) that expands in the ISM (right) at the horizontal velocity V s. In this perturbed configuration, the thermal force, Fth, and the ram force, Fram, are not colinear, so that matter flows (white arrows) from the extended front zone (A), called “hill”, to the receded shell zone (B), called “valley”.

In the text

	Fig. 2 Geometry of the configuration: the numerical domain 1, a quarter of disk (), undergoes an equatorial symmetry with respect to y () defining domain 2, and a rotational symmetry around the x-axis () giving domain 3. It is representative of a 3D axi-symmetric object.
In the text

	Fig. 3 Evolution in time and in space of mass density of shock front region along a radius; density is given by color table in 10−20 kg m−3; panel a: at t = t0 = 3 kyr, panel b: at t ≈ 4 kyr and highlighting transverse velocity with arrows, panel c: at t ≈ 18 kyr, panel d: at t ≈ 53 kyr.
In the text

	Fig. 4 Illustration of regions in which mass variation is calculated during the “overstability” mechanism; panel a: corresponding regions in the unperturbed case, panel b: regions in the perturbed case. Each region, numbered by , should be understood as a volume resulting from a 2π angular rotation around the x-axis.
In the text

	Fig. 5 Mass variation ΔMi (%) by region depending on time. The evolution shows oscillations typical of the VI The regions are numbered by increasing index, i, with respect to the polar angle θ (see the main text). The amplitude of ΔMi is approximately constant with time for each region. We notice that the amplitude of Δ Mi grows with i.
In the text

	Fig. 6 Comparison of profiles with γ = 5∕3 in the SNR at t ≈ 33 kyr for various values of Λ0 with the same color code, W kg−1; panel a: density, panel b: temperature, panel c: pressure, and panel d: velocity profiles. The case Λ 0 = 0 corresponds to the Sedov–Taylor case.
In the text

	Fig. 7 Variation of the radius of the SNR with respect to the time (full curve), in a log–log plot highlighting the transition from the Sedov–Taylor phase (dashed line,n = 0.4) to the asymptotic radiative regime (dotted line, n = 0.3); is expressed in pc and t in 1011 s.
In the text

	Fig. 8 Energy (in 1042 J) radiated by the SNR as a function of time t (in kyr).
In the text

	Fig. 9 Evolution in time and in space of mass density of shock front region undergoing cooling; density is given by color table in 10−20 kg m−3; panel a: at t ≈ 38 kyr, panel b: at t ≈ 43 kyr, panel c: at t ≈ 58 kyr, panel d: at t ≈ 93 kyr.
In the text

	Fig. 10 Transverse mass flux (ρvT) map; panel a; at t ≈ 58 kyr, panel b: at t ≈ 93 kyr. We clearly see transverse fluxes in the dense shell, from the hills to the valleys, characteristic of the VI mechanism. They are still active at a very advanced age.
In the text

	Fig. 11 Massvariation ΔMi (%) by region for the mode l = 24 depending on time, as explained in Sect. 4. The evolution shows oscillations typical of the VI mechanism from t ≈ 3 kyr to t ≈ 50 kyr.
In the text

	Fig. 12 Massvariation ΔMi (%) by region for the mode l′ = 2 × l = 48 depending on time. From t ≈ 50 kyr to t ≈ 140 kyr, the regions corresponding to the mode l′ = 48 present the distinctive oscillations of ΔMi, confirming the development of the VI with a mode l′ twice the initial mode l.
In the text

Fig. A.1 Mass variation ΔMi (%) of the two angular sectors i = 8 and i = 9 already shown in Fig. 5, for five different resolutions. The dotted curves correspond to a low resolution simulation with 2700 × 2700 mesh cells (dx = dy = 3.2 × 1014 m). The full cyan lines are obtained for the resolution used in Fig. 5 (high resolution with 54002 mesh cells and dx = dy = 1.6 × 1014 m). The other curves (full purple, dashed red, and dashed black) have been plotted for both intermediate (respectively for 36002, 43202 cells) and higher resolutions (72002 cells).

In the text