Simulation of the growth of the 3D RayleighTaylor instability in supernova remnants using an expanding reference frame
F. Fraschetti^{1,2,3}  R. Teyssier^{1,4}  J. Ballet^{1}  A. Decourchelle^{1}
1  Laboratoire AIM, CEA/DSM  CNRS  Université Paris Diderot, IRFU/SAp, 91191 Gif sur Yvette, France
2  LUTh, Observatoire de Paris, CNRSUMR8102 and Université Paris VII, 5 Place Jules Janssen, 92195 Meudon Cedex, France
3  Lunar and Planetary Lab & Dept. of Physics, University of Arizona, Tucson, AZ, 85721, USA
4  Institute for Theoretical Physics, University of Zurich, 8057 Zurich, Switzerland
Received 14 June 2009 / Accepted 26 February 2010
Abstract
Context. The RayleighTaylor instabilities that are
generated by the deceleration of a supernova remnant during the
ejectadominated phase are known to produce fingerlike structures in
the matter distribution that modify the geometry of the remnant. The
morphology of supernova remnants is also expected to be modified when
efficient particle acceleration occurs at their shocks.
Aims. The impact of the RayleighTaylor instabilities from the
ejectadominated to the SedovTaylor phase is investigated over one
octant of the supernova remnant. We also study the effect of efficient
particle acceleration at the forward shock on the growth of the
RayleighTaylor instabilities.
Methods. We modified the Adaptive Mesh Refinement code RAMSES to
study with hydrodynamic numerical simulations the evolution of
supernova remnants in the framework of an expanding reference frame.
The adiabatic index of a relativistic gas between the forward shock and
the contact discontinuity mimics the presence of accelerated particles.
Results. The great advantage of the supercomoving coordinate
system adopted here is that it minimizes numerical diffusion at the
contact discontinuity, since it is stationary with respect to the grid.
We propose an accurate expression for the growth of the RayleighTaylor
structures that smoothly connects the early growth to the asymptotic
selfsimilar behaviour.
Conclusions. The development of the RayleighTaylor structures
is affected, although not drastically, if the blast wave is
dominated by cosmic rays. The amount of ejecta that reaches the
shocked interstellar medium is smaller in this case.
If acceleration were to occur at both shocks, the extent of
the RayleighTaylor structures would be similar but the reverse shock
would be strongly perturbed.
Key words: ISM: supernova remnants  acceleration of particles  hydrodynamics
1 Introduction
In young supernova remnants (SNR), the dense shell of material ejected by the explosion and decelerating in a rarefied interstellar medium (ISM) is expected to be subject to hydrodynamic instabilities (Shirkey 1978; Gull 1975,1973) of RayleighTaylor (RT) type. These instabilities modify the morphology of the SNR causing a departure of the ejecta from spherical symmetry. They manifest themselves as fingerlike structures of material protruding from the contact discontinuity between the two media into the ISM heated by the forward shock, as shown by numerical simulations (e.g. Dwarkadas & Chevalier 1998; Chevalier et al. 1992; Wang & Chevalier 2001; Dwarkadas 2000). During this process, the shocked ejecta and the shocked ISM remain two distinct fluids. The Xray observations of Tycho's SNR (Warren et al. 2005) exhibit structures that are consistent with these effects. Deviations from spherical symmetry may also be caused by initial asymmetries in the explosion of the progenitor or local inhomogeneities in the circumstellar medium. In Cassiopeia A, the spatial inversion observed by the Chandra observatory of the iron and silicon layers (Hughes et al. 2000) provides a strong indication of an asymmetric explosion. An even more radical example of deviation from spherical symmetry is SN 1993J, a stellar wind case, where the optical and radio observations can be reconciled by assuming a strong asphericity for the preexisting progenitor activity (Bietenholz et al. 2001).
Twodimensional hydrodynamic simulations of SNRs can take into account the RT instability. Threedimensional simulations have shown an enhancement of smallscale structures and more severe deformation of the reverse shock surface (Jun & Norman 1996); the perturbation seems to grow faster by than in the twodimensional case (Kane et al. 2000). We chose to pursue threedimensional hydrodynamic numerical simulations, focusing on the deviations from spherical symmetry of the SNR in the ejectadominated phase induced by the RT instabilities.
Several distinct physical processes occur in young SNRs, for instance, the dispersion of synthesized materials in the circumstellar medium and the propagation of collisionless shocks. Galactic SNRs are also considered to be a strong candidate source of galactic cosmic rays up to the knee, i.e. 10^{15} eV (Lagage & Cesarsky 1983), since the rate and energy budget can account for the galactic energy density of cosmic rays; however, a direct identification of SNRs as sources of cosmicray nuclei is still lacking.
In the present paper, we aim to investigate with hydrodynamic equations the growth of RayleighTaylor instabilities in the context of efficient particle acceleration. We adapt to SNR physics the hydrodynamic version of the AMR (Adaptive Mesh Refinement) code RAMSES (Teyssier 2002), designed originally to study the largescale structure formation in the universe with high spatial resolution. The first application, considered here, is to study the effect of the growth of RT instabilities on the profile of the hydrodynamic variables by considering the evolution of a full octant of the SNR, i.e. a larger angular region than previously considered (Blondin & Ellison 2001). We do not introduce any seed perturbation into the 3D radial velocity field, so that any departure from spherical symmetry is naturally produced by the numerical fluctuations; this is supported by the nonlinear phase of the instability being insensitive to initial perturbations (Chevalier et al. 1992).
Blondin & Ellison (2001) described the efficient acceleration of particles by changing the effective adiabatic index throughout. A higher compression ratio at the shocks was found to only slightly affect the growth of RT instabilities. However, the suggestion that the reverse shock can efficiently accelerate particles continues to be debated (Ellison et al. 2005). In this paper, we study the impact of cosmicray acceleration on the instabilities using a simplified description that assumes that the adiabatic index is 4/3 in the shocked ISM region, but remains 5/3 inside the contact discontinuity (surface separating ejecta and ISM). This simulates the case of a cosmic raydominated blast wave and gasdominated reverse shock.
To follow the expansion of SNRs, we use synergically two numerical approaches: AMR and the Moving Computational Grid. As a result, the largescale turbulence in SNR can be simulated, allowing for an accurate description of the instabilities.
The hydrodynamic equations can be formulated with respect to two distinct classes of coordinate systems: Eulerian, i.e., fixed spacecoordinates in time, whose major concern is the production of numerical diffusion due to the advective terms; and Lagrangian, i.e., coordinates comoving with the bulk fluid, which is free in principle from numerical diffusion, but possible grid distortions require a rezoning which produces new numerical diffusion. Therefore, Eulerian coordinate systems are generally selected for multidimensional flow simulations. The numerical approach of Moving Computational Grid consists of using a computational grid comoving, or quasicomoving, with the hydrodynamic flow to minimize the local fluid velocity. The idea of a computational grid adapted to follow the global motion of the fluid was first proposed in cosmological numerical simulations by Gnedin (1995). However, further analysis (Gnedin & Bertschinger 1996) highlighted problems in the coupling of the hydrodynamics solver with the gravity solver. Strong mesh deformations were also found in the gravitational clustering. In the approach of Moving Computational Grid, the mesh moves continuously and in addition to a full transformation of coordinates (position and time), a transformation of hydrodynamic variables (density, velocity and pressure) is performed. In contrast, in a purely Lagrangian approach the spacecoordinates are comoving with the bulk flow, leaving the hydrodynamic quantities unchanged.
The main disadvantage of simulating a SNR in a Eulerian fixed computational grid is that in the bulk flow of the SNR the total energy is roughly equal to the kinetic energy. Therefore, the thermal energy, computed as the difference between total and kinetic energies, can be very small leading possibly to local negative thermodynamic pressure. An algorithm is required to control the sign of pressure.
We apply a combination of the AMR with the Moving Computational Grid from the young phase, starting soon after the selfsimilar profile of Chevalier (1983,1982) has been established. The modification introduced enabled us to follow the evolution of one eighth of the volume of a young SNR, a large volume with respect to previous 3D simulations (e.g. Jun & Norman 1996), and for a long interval of time, namely until the transition to the SedovTaylor phase. We do not provide all the details of the RAMSES code, which can be found in Teyssier (2002), but a description of the modifications introduced here is outlined.
The plan of the paper is the following. In Sect. 2, we discuss the initial conditions adopted for a young SNR, namely the mapping of the unidimensional solution of the hydrodynamic equation over a cartesian grid. In Sect. 3, we present the hydrodynamic equations for the SNR flow, written in both the laboratory frame and the reference frame which is comoving with the contact discontinuity. In Sect. 4, we sketch the main steps of the implementation of the Godunov method in RAMSES, we discuss the modifications and the new variables introduced: the active variable to change locally the equation of state; the passive scalar f, which traces the surface of the contact discontinuity; and the ionization age . In Sect. 5, we present the results of SNR simulations: the growth of the RT structures and the elongation in time; the influence of a cosmicraydominated blast wave on the RT instabilities. In Sect. 6, we discuss the main findings and present additional numerical tests. In Sect. 7, we conclude and present forthcoming applications.
2 Initial conditions for selfsimilar expansion of young SNRs
Supernovae observations have confirmed (for the case of SN 1987A, see Arnett 1988) that at early times, typically a few hours after the supernova explosion, the matter density of the outer ejecta reaches a powerlaw profile as a function of radius given by . A kinematical study of the preSedovTaylor SNR of Kepler (Vink 2008a) confirmed the powerlaw density profile. The difficulties in attaining a measurement of n makes the comparison of theory with observations difficult. The value of n depends on the properties of the progenitor star and the acceleration of the shock during the explosion. The value n=7, or alternatively an exponential density profile for ejecta (Dwarkadas & Chevalier 1998), is commonly used to describe type Ia supernovae, which are believed to result from the explosion of a C/O white dwarf accreted in a binary system (Hoyle & Fowler 1960). The value n=9 is believed to describe type II supernovae, which result from the corecollapse of a massive star. If the density of the ambient medium has a powerlaw behaviour, namely , selfsimilar solutions can be found for the density, velocity, and pressure of the interaction region in young SNRs (Chevalier 1982); from dimensional analysis, the time evolution of the radius of the contact discontinuity has the form (Chevalier 1982) , where .
During the expansion, a colder and denser fluid, the shocked ejecta, slows down in a hotter and less dense fluid, the shocked ISM. Seen in the reference frame of the contact discontinuity, this deceleration is equivalent to a gravity field pushing the shocked ejecta toward the shocked ISM. This instability is at the basis of the familiar RT structures (Chevalier et al. 1992). The structure of the interaction region between the ballistically expanding ejecta and a uniform ISM in young SNR was first explored numerically by Gull (1975,1973), who described how the formation of optical filaments is brought back to the RT instability and provided a qualitative explanation of the amplification of the magnetic field in the interaction region. Given spherical initial conditions, we investigate the growth of RT instabilities in different situations.
In this paper, a steep powerlaw density ejecta is assumed to expand into a uniform ISM (s=0).
In the innermost part, the common cutoff in radius
is adopted to introduce the density plateau for the ejecta core
where , and g^{n} is the normalization constant in Chevalier (1982). As expected, dilutes from the initial high values in the inner core down to a negligible value during the SedovTaylor phase.
At a given time ,
the physical parameters that fix unequivocally the profiles of , ,
and P are the total kinetic energy
of the explosion, the total ejecta mass
,
the circumstellar medium density
and the indices of the powerlaw density profile of ejecta and ISM, namely (n,s). The outer radius of the inner plateau is given by
(2) 
A SNR can be defined to be young when the mass of interstellar matter swept up is a small fraction of the mass of ejecta given by . A more quantitative criterion is that the young SNR phase finishes when the reverse shock radius has reached the radius of the core of the remnant.
The selfsimilar profiles of mass density , fluid velocity and pressure P for a young SNR (Chevalier 1983,1982) are used as initial conditions (see Fig. 1). By assuming spherical symmetry, the selfsimilar solution is mapped in a 3D computational grid. At the initial time t_{0} = 10 yr, the selfsimilar profiles are assumed to be wellestablished. It is reasonable to assume that such a spherical symmetry has not yet been significantly modified by the convective motions grown between the explosion of the supernova and the time t_{0}. Given the value of t_{0}, we can fix and therefore the density profile is determined down to r=0, where the origin of the explosion is located. We assume the reasonable values of SNR parameters E_{0} = 1.6 10^{51} erg and , where 10^{33} g is the solar mass and the circumstellar medium density amu cm^{3}. The only information concerning the progenitor of the SNR is therefore contained in the index n and the mass . For the previous set of parameters, we find the initial radii of the contact discontinuity, forward shock and reverse shock to be, respectively, pc, pc, and pc.
3 Equation of hydrodynamics in accelerated frame
3.1 Euler equations in laboratory frame
In an inertial cartesian (laboratory) reference frame
,
the nonrelativistic Euler equations of fluid dynamics, in the
absence of viscosity and dissipative processes, are written in the
quasiconservative form
where is the mass density, u_{i} is the ith component of fluid velocity, P is the pressure, is the total energy density, and is the internal energy density. The solution of these equations requires the use of an equation of state to find the quantity given the quantity (Synge 1957). We use the common shortcut for the adiabatic case, consisting of defining the parameter by means of the equation . In the nonrelativistic case, can be identified as the ratio of the specific heats of the gas ( ). If the pressure is dominated by the relativistic particles, . The active scalar is introduced to treat in a purely hydrodynamic way the effect of efficient particle acceleration at the forward shock (see Sect. 5.2).
Equation (3) is written in conservative form, except the socalled ``quasiconservative'' form of equation for . If the conservative variable is used, higher order schemes of shock capturing applied to multicomponent flows can produce oscillations at the interface of different fluids (Shyue 1998; Johnsen & Colonius 2006).
For a shock propagating into an ideal gas, the compression factor at the shock can be expressed as (Zel'dovich & Raizer 2002), where () is the downstream (upstream) density of the shock and is the Mach number in the unshocked medium. In the limit of a very strong shock, as in young SNRs, and , depending only on the thermodynamic properties of the medium.
The passive scalar f, representing the mass fraction of ejecta, is initially defined as
where is the radius of the contact discontinuity. The gradient of the function f traces the surface of contact discontinuity, namely the fingerlike structures arising because of the RT instabilities, and f is propagated by means of
The Xray spectra of SNRs usually exhibit significant emission lines of heavy elements, such as Fe in Kepler (Smith et al. 1989). Therefore, the computation of the Xray spectrum of SNRs requires an evaluation of the ionization produced by electrons after the passage of the reverse shock. In this paper, we do not compute the ionization state by coupling the hydrodynamics equations with the equations for the evolution of ionization and recombination (for a 1D example of this see Rothenflug et al. 1994). We follow the ionization age, which allows us to compute a posteriori the ionization structure assuming a constant temperature in the fluid element.
An additional passive scalar representing the ionization age, ,
was introduced into the code, defined as
where is the time since crossing the forward or reverse shock at point and is the number density of postshock electrons in the heated plasma. As a passive scalar, does not cause any variations in the dynamics of the SNR. The ionization age satisfies the equation
where is the Heaviside function, is the pressure in the unshocked ISM, and we can assume for the pressure in the free expansion region, i.e., , that , therefore Eq. (7) is homogeneous outside the interaction region. The constant B is the threshold ratio for the detection algorithm of the interaction region. Since the pressure in the interaction region satisfies , we can choose B=10^{2}. No failures in the detection algorithm have been found for lower values of B down to B=10.
3.2 Euler equations in comoving reference frame
Now we consider a noninertial reference frame . We deal with the simplified case of a purely radially contracting accelerated motion, which was introduced in Martel & Shapiro (1998) as supercomoving coordinate system; the more general case in the presence of rotation is treated in Poludnenko & Khokhlov (2007). The frame is defined by the transformation
where and are two scaling parameters. The scale factor a(t) must be a nonvanishing twicedifferentiable function of time only. Once the transformation law of the three fundamental quantities, namely length, time and mass, is fixed, the tranformation law of all the others can be found. From Eqs. (8) and (9), we obtain the transformation law of the velocity vector
where . The velocity is decomposed into the velocity in the inertial frame and the relative velocity of the bulk motion. To establish a computational grid comoving with the hydrodynamic flow, we choose the expansion law
where we consider the simple case of being independent of time. In the cosmological framework, the transformation to supercomoving coordinates maps the uniform and isotropic Hubble flow in a matterdominated universe to a uniform matter distribution at rest with constant thermodynamic properties, in the adiabatic case, namely . In that case, the aim is to focus on deviations from the Hubble flow that drive the structure formation in the universe. In the case of selfsimilar expansion of SNR, the transformation of Martel & Shapiro (1998) has the effect of converting the mesh from Eulerian to quasiLagrangian, which is exactly comoving with the contact discontinuity while the selfsimilar expansion applies, but does not experience grid distortions as in the case of purely Lagrangian mesh.
From Eqs. (8) and (9), we also have
In 3D, the choice , for any value of , preserves the invariance of the mass conservation equation. However, as shown explicitly in Poludnenko & Khokhlov (2007), only the choice
and for a perfect monoatomic gas with preserves the invariance of the Euler equation (Eq. (3)). In this paper, we report results obtained by fixing and in Eq. (15) and , fixing the value of according to the physical parameters of the specific problem. Hereafter, the choice a_{0} = 1 is applied without reducing the generality of the treatment. In the reference frame , the Euler equations (Eq. (3)) become
where , and
Here, for a generic function we have
(18)  
(19) 
In most physical cases, the bulk motion is accelerated; therefore in the momentum equation written in the noninertial reference frame a noninertial force term will appear on the righthand side. In the righthand side of Eq. (16), both the second line and the first term of the third line represent the noninertial force; the second term of the third line, treated as a source term since it acts as the gravitational potential term in the cosmological framework, may create an energy sink, losing the numerical advantage of the conservative form (see Sect. 4 for details).
Since it is an advection equation, Eq. (5) is unchanged by the choice of Eq. (15)
Equation (7) changes as follows
In this paper, we consider the early phase of SNR expansion and the transition to the SedovTaylor phase, thus the radiative cooling is negligible with respect to the total energy content and the approximation of adiabatic expansion of a perfect fluid is justified. We assume that the fluid mainly consists of fully ionized hydrogen, therefore . We notice that for the last term in the third relation of Eq. (16) disappears. If , as in Sect. 5.2, this term can be regarded as a cooling or heating term, representing a microscopic form in which the internal energy of a nonmonoatomic perfect gas can be deposited in the fluid as vibrorotational degrees of freedom.
The physical time t is related to the time
for
and by
In the cases considered, ( and ), thus the factor for every . Therefore, the evolution in ranges from , for t = t_{0}, to for . This squeezing of the total time interval depends on the parameters of the progenitor of the SNR. We also note that, because of Eq. (9), since , the time in is compressed relative to the physical time in by a factor , reducing the time step of the numerical integration.
Figure 1: Angleaveraged radial density profile at t=50 yr for , or (n,s) = (7,0), and with 13 levels of refinement (128^{3} in yellow, (128256)^{3} in green, (128512)^{3} in black). For comparison, the selfsimilar profile is shown in red. The density profiles are normalized to their values at the forward shock. From this diagram, it can be seen that the localization of the forward shock and its compression ratio are not strongly dependent on the maximal refinement level. 

Open with DEXTER 
4 Numerical code
4.1 Setup
The simulations were performed with the hydrodynamics version of the code RAMSES (Teyssier 2002). RAMSES was originally developed to study, with the AMR technique, structure formation in the universe at high spatial resolution. The solver is based on a second order Godunov scheme. We modified the constitutive equations of RAMSES as shown in Eqs. (16) and (17) to solve the hydrodynamics flow in an accelerated reference frame expanding according to Eqs. (8)(10). The modifications concern the hydrodynamic solver routines and the external gravity routine, because the change of reference frame introduces a noninertial force in the Euler equations, which can be treated as a source term in the momentum and energy equation. We also introduced an additional variable , treated as a passive scalar in Eqs. (3) and (16), to change locally the equation of state and two more passive scalars (see Sect. 3.1 for details).
The simulations were performed over three levels of refinement, of value l=79, corresponding to a cartesian grid with a number of cells 128^{3}, 256^{3}, and 512^{3}, respectively. The effective maximal resolution attained in the interaction region at the setup of the simulation is 0.4/512 10^{4} . Our 3D numerical simulations were carried out across one octant of a sphere.
The refinement algorithm is applied to the surfaces of the contact discontinuity and the two shocks. The surface of the contact discontinuity is detected by means of the gradient in the tracing function f (see Sect. 3.1), rather than the density gradient; in this way, the effectiveness of the detection is independent of the parameters (n,s). The forward and reverse shocks are detected by measuring the pressure discontinuity. The advantage is that the pressure discontinuity is greater than the discontinuities in density or radial velocity and independent of the assumptions about the density powerlaw index in the free expansion medium and the equation of state of the gas in the interaction region. Figure 1 shows the angleaveraged radial density profile at t=50 yr for , or (n,s) = (7,0), and with 13 levels of refinement. The resolution of the RT structures is crucially improved by applying the AMR.
During the expansion of a young SNR, the only dynamically interesting region, represented by the interaction region, occupies a fraction of about 40% (and constant during the selfsimilar regime) of the whole remnant volume. Moreover the KelvinHelmoltz instabilities triggered by the growth of the RT fingerlike structures during the ``selfsimilar'' phase lead to a strongly turbulent phase. Therefore, AMR hydrodynamic simulation that attempts to spatially resolve the growth of the RT instabilities must refine over an increasingly large relative volume. The contact discontinuity turns out to be the most computationally expensive region of the flow. In contrast, the core region within the powerlaw ejecta, where the value of the uniform density is fixed by mass conservation, can be easily handled at the lowest possible refinement level. In this context, the AMR technique, combined with the approach of Moving Computational Grid, is found to be suitable.
We define the boundary conditions of the conservative variables in the corresponding six ghost zones. Since we treat 1/8 of the volume of the remnant sphere, in some ghost zones (three) the boundary condition is inflow, produced by the advection of the interstellar material, and in some other ghost zones (three) the boundary condition is reflexive, according to the location of the interface of the ghost cell. Changes in the order of the assignment of the ghost cells significantly change neither the resulting intercell flow at the shock nor the physical content of the result. In the laboratory rest frame, the ISM is cold and stationary: . In the shock frame, this implies that the ISM is advected onto the SNR shock surface with a velocity field given by (see Eq. (11)). Therefore, a linear interpolation (second order) for the reconstruction of the velocity field at the boundaries more clearly describes the inflowing material.
The time step is controlled using a standard Courant factor stability constraint. Based on our expanding coordinate system, we also included a time step limitation, where a(t) is not allowed to vary by more than over the time step.
The cartesian coordinate system that we use here is not welladapted to spherically symmetric objects such as SNRs and spurious geometrical effects are expected. It has long been known that an Eulerian code should take care of its violation of Galileian invariance: a turbulent phenomenon observed in a fixed computational grid can disappear if it is observed from a computational grid moving at constant velocity with respect to the fixed grid. The goal of our comoving coordinate system is precisely to eliminate this effect. This improvement does, however, have the problem, which is common to all shock capturing schemes dealing with stationary shocks, of the socalled ``oddeven decoupling'' and the associated carbuncle phenomenon (Peery & Imlay 1988; Quirk 1994; Pandolfi & D'Ambrosio 2001), namely a disturbance in the computation of the flow emerging when the shock direction is aligned with one of the grid coordinate directions. This effect may be amplified when a flow converges in a cell from different directions, as in the case of spherical flow on a cartesian mesh. The result can strongly deform the shock surface and form ripples or bubblelike regions of local subdensity. Local changes to the solver were proposed (Kifonidis et al. 2003), based on an algorithm for the detection of the shock. Reconstruction of the interface values of the state variables is performed, but in the cells detected by that algorithm a more diffusive solver is used, such as HLL, HLLC or LaxFriedrichs (Toro 1999), while in all the remaining regions a conventional Riemann exact solver or an approximate Roe solver can be used. Solving this longstanding issue is beyond the scope of this paper, so we did not implement these possible fixes here.
4.2 Integration method
We briefly recall the main steps of the integration method of RAMSES (Toro 1999).
In the absence of a gravitational potential,
the threedimensional hydrodynamic Euler equations can be written
in the conservative form as (see Eq. (3))
where has components , has components , and both i and j have values from 1 to the space dimension d=3. The spacetime integration is performed by defining a cell centered on (x_{i}, y_{j}, z_{k}) of volume defined by the coordinates and a time interval by , which is not constant through the simulation. In the present section, the quantity n indicates the nth time step of the integration algorithm. The averaged, cellcentered state is defined by
(24) 
Apart from a smooth flow of the conservative variables , the discretization of space imposes the solution to a unidimensional Riemann problem along all the space directions. Inside the cell V_{i,j,k}, the solution can be reconstructed with distinct algorithms, for instance constant (first order), linear (second order), to obtain the values at the interfaces , , , , , and . The intercell flux is then computed by using the values at the interface. The timeaveraged intercell flux is defined by integrating over the planes separating neighboring cells
(25) 
(26) 
(27) 
The conservative system can then be written
(28) 
This integral form remains exact for the corresponding Euler system since no approximation has been made. Since the intercell flux is computed at time , the method is second order in time. The main changes to RAMSES reported here concern the definition of the inflow of interstellar matter into the simulation box, as explained above, and the introduction of the source term . RAMSES allows us to define an analytical acceleration vector, which is treated as an effective gravity vector. We have used this option, the grid motion being defined analytically by Eq. (12). The three independent physical quantities, i.e., length, time, and mass, are rescaled according to Eqs. (8)(10). All the quantities adopted until the end of the present section are computed in the frame.
In presence of a source term, Eq. (23) becomes
where denotes a numerical approximation to the cellaveraged value of at time for the cell (i,j,k), and the space indices (i,j,k) and the time index (n) are kept for simplicity the same in the frame. The numerical discretization of the Euler equations with source terms S^{n+1/2}_{i,j,k} is given by:
The time centered fluxes , , across cell interfaces are computed using a secondorder Godunov method. According to Eq. (16), the gravitational source terms are given by
where is defined in Eq. (17) and is the internal energy density. The second term on the righthand side of Eq. (31) is equal to zero if or, equivalently, .
Although the integration of the source term S^{n+1/2}_{i,j,k} in Eq. (30) is of second order in time, higher order methods, such as Rosenbrock (Kaps & Rentrop 1979), can be chosen.
5 Application to SNR
5.1 Growth of RayleighTaylor structures
Figure 2: Angleaveraged radial density profile at various ages for , or (n,s) = (7,0), and with 1 and 3 levels of refinement (128^{3} in yellow, 128^{3}512^{3} in black). For comparison, the red line shows the selfsimilar profiles normalized to their forward shock values, Chevalier profiles for the first five panels and Sedov for the last four panels. The departure from the ejectadominated phase, occurring when the reverse shock encounters the plateau outer radius , and the transition to SedovTaylor phase are shown. The inward motion of the reverse shock, initially only in the comoving frame of the contact discontinuity and later on in the laboratory frame as well, is also depicted. 

Open with DEXTER 
Figure 3: Temporal evolution of the forward shock and reverse shock radii for , or (n,s) = (7,0), and . Our 3D numerical solution (solid) matches the 1D semianalytical solution (Truelove & McKee 1999), dashed, and, in the asymptotic regime for , matches the selfsimilar solution in the ejectadominated phase (Chevalier 1982), dotted. See Sect. 5.1 for the units and ( pc and yr for the chosen parameters). The spherical symmetry of the reverse shock is strongly modified by the RT instability; here the innermost value of the radius of the reverse shock has been chosen (cf. Sect. 5.1). The ejectadominated solution (dotted) is prolonged indefinitely to emphasize the inward motion of the reverse shock. 

Open with DEXTER 
Many developed codes have satisfactorily tackled the growth of RT instabilities as a mandatory validity test. In contrast to previous treatments, we consider the evolution of a full octant of the SNR. Thus, within a larger SNR volume, a large number of RTinstabilitiy lengthscales, spread over a wider range, sum together. In other words, the growth of RT instabilities is more accurately described not only with high spatial resolution, but also with volumes comparable to the global SNR volume. In particular, since the initial departures from the spherical symmetry are sampled over a larger surface, the initial phase of the RT growth is statistically more accurately described (see Sect. 6.3).
Even in the case of a stationary advective shock, a small volume has been dealt with, in the SNR literature, in comparison with the volume of the interaction region. The description here provides a more comprehensive overview of the dynamics of RT instabilities in the SNR.
The simulations presented here cover three phases of the adiabatic SNR expansion (see Fig. 2): the selfsimilar expansion (Chevalier 1982) with an exponent depending on the properties of the ejecta and the ISM, i.e. on (n,s); the second phase, not selfsimilar, representing the transition to the SedovTaylor regime ( ); and the third phase (SedovTaylor), with a blast wave radius that expands as , for s=0, and . The transition between the two selfsimilar phases is shown by means of various snapshots of the angleaveraged density in Fig. 2. In this case, we assume that , or (n,s) = (7,0), and in the whole simulation box. The inward motion of the reverse shock and the progressive establishment of the SedovTaylor profile is manifest. We extended our computation until values of time t when the inward motion of has reached the inner core, and show that during the transition to the SedovTaylor phase, the angleaveraged profile is not significantly altered by the development of the RT instability. The 1D SedovTaylor density profile is superimposed in the last three panels of Fig. 2, by assuming that the SedovTaylor phase is fully established at t = 5 , where (Truelove & McKee 1999). Compared with the original version of the code in the case of a spherical Sedov blast wave (Teyssier 2002), we note the improvement in the resolution of the compression factor; in that previous work relatively small values of final output times had been considered, the aim being to test the code.
In Fig. 3, the outermost value of the position of the forward shock and the innermost value of the position of the reverse shock in 3D are compared with the semianalytical 1D solution of Truelove & McKee (1999). The units used here correspond to (pc) = and (yr) = , where is the mean mass per hydrogen nucleus assuming cosmic abundances. In particular, the inward motion of the reverse shock is compared with the solution of Truelove & McKee (1999); in the case (n, s) = (7, 0), we find agreement within between our 3D simulation and the 1D solution of Truelove & McKee (1999). The radius of the reverse shock is strongly modified by the RT instability where the innermost value is chosen. In contrast, the radius of the forward shock is unequivocally defined, since it is not perturbed by the RT instabilities for the physically relevant values of (n,s) (see also Fig. 5).
In Fig. 4 (upper panel), the mass fraction , defined in Eq. (4), is shown in a planar slice parallel to the coordinate plane XY at time t=100 yr (left) and t=500 yr (right); in this case, , or (n,s) = (7,0) and . The mixture of the two media, shocked ejecta and shocked ISM, as well as the deformation of the contact discontinuity surface appear at intermediate values of . In the middle panel, the matter density is shown with the same spatial extension and times as in the upper panel. In the lower panel, an estimate of the frequency emissivity produced by bremsstrahlung in purely hydrogen gas is shown, integrated along the line of sight towards the reader. Only the interaction region where the temperature is high contributes to the emissivity. The emissivity is computed based on the assumption that the plasma is optically thin (bremsstrahlung). Therefore, the emissivity is given by , where is the electron number density and is the electron temperature given by the perfect gas equation of state such that , where is the mean molecular weight, 10^{24} g is the proton mass, and 10^{16} erg/K is the Boltzmann constant. For simplicity, in this case we chose , corresponding to the case of fully ionized hydrogen. We implicitly assumed that , where is the ion temperature, while at an early stage the SNR is not yet completely thermalized. A more detailed study should account for inclusion of the evolution of , but this indicative diagram provides nevertheless a starting point.
In Fig. 5, a 3D large view image of the density field in the full simulation box is shown. This picture shows that while the sphericity of the forward shock surface is maintained since the RT instabilities do not reach the shock, the reverse shock surface is strongly modified. The large relative volume of the refinement region can be clearly seen.
In Fig. 6, the streamlines of the velocity field in the reference frame of the contact discontinuity in distinct snapshots show the development of the chaotic flow in the interaction region due to the convective motions. As the time proceeds, the central density decreases and the relative size of the region of turbulence increases.
Figure 4: Snapshots at two distinct ages (t=100 yr at left; t=500 yr at right) for , or (n,s) = (7,0) and . In the upper panel, the mass fraction , defined in Eq. (4), is shown in a comoving slice plane parallel to the coordinate plane XY at height z=0.1 a(t). The mixture of the two media, namely ejecta shocked and ISM shocked, as well as the deformation of the contact discontinuity surface is shown. In the middle panel, matter density is shown in the same plane as in the upper panel. Density values are coded according to the color bars given for each frame. Note the change in the radial scale. In the bottom panel, frequency emissivity due to bremsstrahlung integrated along the line of sight towards the reader is shown, assuming a mean molecular weight , corresponding to the case of fully ionized hydrogen. 

Open with DEXTER 
For a general acceleration field g, the initial growth of
small amplitude perturbations at the accelerated contact discontinuity
between two incompressible fluids is exponential with a growth
rate .
The quantity
is related to the wavenumber k by the relation
,
where
and ,
are the densities of the two adjacent media at the contact
discontinuity. At later times, the growth of the ripples at the
contact discontinuity enters the socalled ``selfsimilar phase''; this
phase was first quantitatively analyzed by Fermi & von Neumann (1953). The selfsimilar growth of the RT instable structures is governed by the equation
where the constant is a dimensionless growth parameter depending on the dimensionality and the geometry of the system, and g is the acceleration. If , the asymptotic quadratic law in time is found to be . A second order polynomial fit was found to reproduce quite satisfactorily the mixing of ejecta and ISM in a young SNR (Dwarkadas 2000). However, in view of the selfsimilar deceleration of the contact discontinuity of an SNR, it is reasonable to assume a timedependent g. We propose that a more accurate result can be found by assuming that the acceleration g is also selfsimilar, namely . The integration of Eq. (32) gives readily
where is the initial radius of the contact discontinuity. The solution in Eq. (33) reproduces our simulation (see Fig. 7). The comparison in Fig. 7 shows a disagreement of a factor of more than 2 with the second order polynomial fit at early times, because the deceleration of the contact discontinuity has not been properly taken into account.
Figure 5: The threedimensional density in the full simulation box of side L = 4.2 pc is shown at t=500 yr for , or (n,s) = (7,0) and . The black line indicates the height of the planar slice at z = 0.1 L used in Fig. 4, upper and middle panels. 

Open with DEXTER 
Figure 6: Snapshots showing streamlines of velocity field superposed to matter density for , or (n,s) = (7,0) and in the frame of contact discontinuity. The planar slice is in all panels z=0. From left to right, in above panel (t=100 yr) the inflow motion of the ISM and the outflow motion of the freelyexpanding ejecta are clearly evident. In the other panels (t=500 yr, t=1000 yr, t=2000 yr), the turbulent flow motion due to RT instabilities is shown. Note that the scale in color bar is log in amu cm^{3} in the first panel and linear in all the other panels. 

Open with DEXTER 
The limited spatial resolution attainable at the contact discontinuity with a cartesian grid makes it hard to verify early time exponential growth (Fermi & von Neumann 1953), prior to the selfsimilar regime.
5.2 Nonthermal particles in the shocked ISM
Multiwavelength surveys have inferred that cosmicray acceleration can be efficient in young SNRs (see Vink 2008b, and references therein). The idea at the basis of particle acceleration at SNR shocks is the Fermi firstorder mechanism of acceleration, namely that particles accelerate by repeatedly crossing the shock and colliding at every passage with magnetic turbulences. Therefore, beside other parameters, the efficiency of acceleration depends on the intensity and the spatial configuration of the magnetic field. The Xray detections of narrow synchrotron emitting filaments at SNR shocks indicate magnetic field intensities up to 1 mGauss (Vink & Laming 2003; Lazendic et al. 2003; Bamba et al. 2003). Magnetic field amplification at shocks has been interpreted as the result of turbulence in the postshock material (Giacalone & Jokipii 2007) or cosmicray induced streaming instability (Bell 2004; Amato & Blasi 2009).
The SNR hydrodynamics can be modified by efficient particle acceleration at the forward shock (Decourchelle et al. 2000). Therefore, a complete numerical description of the SNR expansion should take into account the feedback of the nonthermal population of particles on the dynamics of the SNR, i.e., the system of Euler equations (Eqs. (3)) should be coupled to the equation for the cosmicray phasespace distribution function.
Blondin & Ellison (2001) investigated this scenario by changing the adiabatic index
globally. They found that the development of the RT instability
was relatively unaffected and the main effect was that, because the
compression ratio increases and the shocked region shrinks when
decreases, the RT fingers approached the forward shock and could even reach it for
.
We adopt a similar approach but investigate when particle acceleration
is efficient at the forward shock only. We call this the hybrid case in
the remainder of the paper. A local equation of state of the form
is assumed. The shocked ISM is initialized as a relativistic gas
(dominated by cosmic rays) and the shocked ejecta is initialized
as a nonrelativistic gas such that
the corresponding then being propagated according to Eq. (16).
The initial adiabatic distribution in Eq. (34) should provide a more realistic representation of forward shock accelerating particles than a uniform distribution across the whole interaction region. There has been relatively little observational or theoretical evidence that particle acceleration is very efficient at the reverse shock (Ellison et al. 2005). An interesting result of that simulation is that the density behind the reverse shock is lower than behind the forward shock in the selfsimilar solution, at least for n = 7 (Fig. 8 as opposed to Fig. 1). This may hinder the development of the RT fingers.
6 Discussion
6.1 Intrinsic issues in the stationary advective shock problem
In simulating stationary advective shocks, numerical instabilities related to the ``carbuncle phenomenon'' can be found by using the Roe solver. The carbuncle phenomenon was first observed by Peery & Imlay (1988) in blunt body simulations using Roe's method and numerically studied in detail by Quirk (1994), and has remained a crucial numerical issue (Pandolfi & D'Ambrosio 2001; Xu & Li 2001) until recently (Nishikawa & Kitamura 2008). The carbuncle phenomenon is considered to be a spurious numerical effect affecting hydrodynamic simulations when the intercell flow is aligned with one of the coordinate directions of the computational grid; it shows up when a supersonic shock, or more generally a discontinuity in the flow variables, is comoving with the computational grid. The Riemann solver produces unphysical ripples of the shock and in the postshock region, which at late times can modify the structure of the shock surface itself. The numerical methods that are likely to produce carbuncle phenomenon have not yet been unequivocally identified, an attempt to classify them being Dumbser et al. (2003).
The socalled ``Hcorrection'' (Stone et al. 2008), which consists of introducing an anisotropic dissipative flux in the direction perpendicular to the direction of propagation of the shock, cannot be applied here. We did not add a diffusive term in the intercell flux, which is usually adopted to avoid the growth of numerical instabilities. The addition of artificial viscosity could damp physical effects emerging during the development of the RT instabilities, which we are interested in studying in detail. For the same reason, we did not use more diffusive solvers such as HLL or HLLC or LaxFriedrichs (Toro 1999).
In the case of supernova corecollapse simulations, a possible solution has been to compute the intercell flux in the zones of strong shock with an Einfeldt solver (Quirk 1994,1997), while all the other grids in the computational box are treated with a less diffusive Riemann solver (Kifonidis et al. 2003). We preferred not to introduce a solverswitching algorithm to avoid introducing other spurious effects possibly arising from the matching of the intercell fluxes.
To avoid the numerical instability occurring at the stationary shock, we used the exact solver. The failure of the Roe solver manifests itself in local fluctuations of density in regions behind the shock and strong deformations of the shock surface, which grow as the nonlinear regime of RT instabilities is reached; consequent unphysical effects due to velocity shear, resembling KelvinHelmholtz instabilities, develop at the same time. Both effects are absent if the shock is not stationary with respect to the computational grid and depend on the chosen simulation parameters, for instance CFL number or interpolation scheme of conservative variables; thus they must be attributed to the choice of the solver. The choice of the exact solver did not significantly increase the computational time.
6.2 Quasicomoving grid
A different solution to the carbuncle problem could involve defining a computational grid that drifts with respect to the SNR, slowly enough to follow the expansion of the remnant for a long time and to avoid boundary reflections possibly occurring with a remnant overtaking the grid. However, this is an even less viable solution since, as shown in the following, the angleaveraged density profiles strongly deviate from the expected compression ratio at the two shocks. In contrast, as shown earlier, in the case of an exactly comoving computational grid, the postshock ripples average to zero.
Figure 7: The protruding RT structures measured in our simulations at different times are quite well reproduced by the exact solution H(t) with (see Eq. (33)). The edge of the RT structures is defined by . The earlytime exponential growth can hardly be verified with the available spatial resolution. A second order polynomial fit is also shown for comparison. 

Open with DEXTER 
Figure 8: Angleaveraged radial density profile shown at distinct ages for , or (n,s) = (7,0) in the hybrid case in black, is compared with the corresponding initial profile, in red. 

Open with DEXTER 
Independently of the initial density profile, i.e., of the value of , the selfsimilar regime slows down towards the SedovTaylor solution, which has a selfsimilar exponent equal to . Therefore, the contact discontinuity in the computational grid will progressively acquire a drift velocity depending on the physical parameters of the remnant. To test the numerical stability of the shock, we introduced a ``quasicomoving'' grid, whose small drift with respect to the contact discontinuity velocity is parametrized by : .
Figure 9: Matter density in a comoving planar slice parallel to the coordinate plane XY at height z=0.1 a(t) (t=100 yr at left; t=500 yr at right) for , or n=7, s=0, and . Density values are coded according to the color bars given for each frame. Note the change in the radial scale. In this case, . Since , the computational grid expands faster than the SNR and the volume fraction occupied in the simulation box by the SNR decreases. It is evident that even a small drift of the computational grid from the shock velocity erases the effects of the numerical instability possibly emerging behind the shocks. 

Open with DEXTER 
The quasicomoving grid, even with the Roe solver, smears out the ripples produced behind the forward shock surface by the stationarity of the computational grid, therefore locally is oddeven instabilityfree (see Fig. 9). However, the angleaveraged density profile, even at the first timesteps, is affected by an error in the compression factor at the two shocks (see Fig. 10, upper panel), the position of the shock being unchanged. The redistribution of the mass through the two shocks in the presence of a grid drift (overcompression at reverse shock and undercompression at forward shock) is to be attributed to the mapping of a spherically symmetric object on a cartesian grid. It manifests itself during the transition to the SedovTaylor phase in terms of an error of in the resolution of the compression at the shock, and should be considered an intrinsic limit of the approach used here. This error does not depend on the number of levels of the AMR (see Fig. 10, lower panel). As a consequence, even if a larger number of computational cells is considered, the angleaveraged value at the shock does not fulfill the value of the compression predicted by the RankineHugoniot conditions. In contrast, for an exactly comoving grid, even with the Roe solver, strong deformations of the shock surface are triggered by the numerical instability but the compression at the shocks is preserved. The same difference of results between the exactly comoving and quasicomoving has been found with the exact Riemann solver, shown in Fig. 10.
As a consequence, during the transition to the SedovTaylor solution, the generation of a numerical error must be expected in the compression at the shock. A possible solution is the introduction of a grid velocity given by .
6.3 Cosmicraydominated blast wave
Figure 10: Angleaveraged radial density profiles between the reverse shock and the forward shock at time t=12 yr = 0.0063 , with , or ( n,s)=(7,0), and . The AMR is performed over three levels of refinement. The radial coordinate is normalized to the forward shock position corresponding to ; the density is normalized to the value at the shock. Upper panel. The profile in the case of the exactly comoving grid is compared with the profiles having a quasicomoving grid . The initial profile, at t_{0} = 10 yr, is also depicted for comparison. It is evident that even small values of produce oscillations and an error of the order of in the compression factor for . Outside the interaction region, the profiles are superimposed. Lower panel. The curves without AMR and with AMR are compared with the initial profile at t_{0}=10 yr. In this case, . The independence of the error in the compression factor at the two shocks from the space resolution is clear. 

Open with DEXTER 
Figure 11: Angleaveraged radial profile of the mass fraction of the ejecta f shown at distinct ages for , or (n,s) = (7,0); the cases with uniform , and the hybrid case are compared. The radius is in units of the corresponding radius of the contact discontinuity. The vertical lines on the right of the panels indicate the corresponding positions of the forward shock, which at this age still meets the selfsimilar expansion (see Fig. 3 for the uniform case). The vertical line on the left of the panels indicates the position of the reverse shock. 

Open with DEXTER 
Figure 12: The protrusion of RT structures measured in our simulations in units of the corresponding selfsimilar position of the contact discontinuity, whose surface is directly modified by the instabilities; the uniform and the hybrid cases are compared ( ). The cosmic rays affect the elongation of the instabilities from the very beginning by a factor of the order of . 

Open with DEXTER 
For s = 0, a convenient case to consider when analyzing the growth of RT structures is n = 7, for which the width of the interaction region is large, thus the elongation of the instabilities can be studied for a longer time and the forward shock remains unaffected. The effect of changing the equation of state on the growth of the instabilities is represented in Fig. 11. It shows how far the RT fingers reach during their growth (100 yr) and when they are fully developed (500 yr). At 500 yr, the fraction of ejecta material transported well into the shocked ISM (at , say) is lower in the hybrid case. Nevertheless, the outermost ejecta reach nearly as far, and very close to the forward shock.
In the hybrid case, the position of the reverse shock is unaffected by a relativistic gas at the forward shock, and the departure from the selfsimilar position is less than . If the relativistic gas is uniformly distributed across the whole interaction region, the departure of the reverse shock from the selfsimilar position is of the order of . Similar deviations have been found by Blondin & Ellison (2001), who considered however a relatively small angular region in the cases (n,s) = (7,2) or (n,s) = (12,0). The bumps in f(r) found at originate in the RT structures themselves, independent of the distribution. In the two cases and , the width of the shocked ISM region ( ) shrinks from 0.18 to 0.10. Comparable shrinking has been found in the angleaverage analysis of Xray observations of Tycho's SNR (Warren et al. 2005).
As shown in Fig. 12, the elongation of the RT structures in the hybrid case decelerates relative to the nonrelativistic case. The deceleration is caused by the higher density gradient ahead of the contact discontinuity (cf. Figs. 8 and 1). This can be easily understood because the RT growth rate depends through A on the average ratio of the two densities in the interaction region, i.e., (cf. Fig. 8), thus lower density ratios produce lower growth rates (see Fig. 12). The higher compression in the ISM shocked region in the hybrid case produces the slowing down shown in Figs. 11 and 12. When on both sides as in Blondin & Ellison (2001), remains larger than 1 and the RT fingers reach further than in the hybrid case at .
The region which drives the instability at very early time is that part of the density profile where density decreases outwards (see Fig. 1). This is entirely on the ejecta side and therefore unaffected by changing in the shocked ISM (see Fig. 11). In other words, at a very early stage, when the instabilities grow exponentially, the process occurs locally in the outer ejecta and the cosmic rays cannot affect it. However, particle acceleration will be able to affect RT structures as soon as they reach the shocked ISM, since the density gradient of the shocked ISM is higher for higher compression ratios. For the reasons explained above, a phenomenological equation of the type of Eq. (32) can only qualitatively reproduce the global behaviour. We could not find a more accurate description in this case. A comparison of Figs. 12 with 7 of Blondin & Ellison (2001) might indicate that the large fluctuations they observed in the latetime radial extent of RT structures are not an intrinsic property of the instability but a result of the too small physical size of the region considered.
7 Conclusions
We have adapted the AMR code RAMSES, which is based on a second order Godunov method in an expanding frame called ``supercomoving coordinates'', to follow the evolution of SNRs. In this approach, not only the spacetime variables have been modified but also the hydrodynamics variables. The comoving coordinate system allows an eighth of the total volume of the SNR to be described, which is much larger than previously considered. A longer time interval can be investigated, of the order of thousands of years, until the transition to the asymptotic SedovTaylor phase. Such a large volume allows the convective instabilities to be modeled more accurately since it considers not only the shortest wavelengths, which have the greatest growth rate, but also the longer wavelengths, which grow more slowly but make a significant contribution to the morphology of the SNR over timescales of thousands of years. Our larger spatial sampling allows a statistically more accurate description of the instabilities.
The great advantage of the method adopted here is that it minimizes the velocity of the fluid relative to the grid, and the numerical diffusion at the contact discontinuity, where the instability should develop. In the comoving reference frame, the contact discontinuity, in the absence of distortions due to convective instabilities and before the transition to the SedovTaylor phase, would be exactly stationary in time. The analytical noninertial term resulting from this coordinate transformation is strictly equivalent to a gravitational acceleration: the Euler equations have therefore been integrated with the corresponding additional source terms.
The elongation of the RayleighTaylor structures slowly reaches the asymptotic behaviour , by direct solution of the selfsimilar theory according to the assumption of a selfsimilar acceleration instead of a constant acceleration.
We have presented a simple way to numerically investigate the effect of efficient particle acceleration at the forward shock, approximating the shocked ISM as a relativistic gas and the shocked ejecta as a nonrelativistic gas, with a different adiabatic exponent in each fluid. The density behind the forward shock may be higher than at the reverse shock in that case. A deceleration in the protruding of RT structures is caused by the higher compression of the shocked ISM. The conclusion of Blondin & Ellison (2001) that the RT fingers can travel very close to the shock in the presence of accelerated particles remains valid. The elongation of the instabilities has here been more precisely quantified. We propose that this will allow us to understand why ejecta are found very close to the blast wave in the remnant of SN 1006 (CassamChenaï et al. 2008). The setup of the code presented here will allow future studies of the backreaction of particle acceleration on the SNR evolution.
AcknowledgementsThe numerical simulations have been performed with the DAPHPC cluster at CEA, DSM, for which the technical assistance is gratefully acknowledged. F.F. wishes to thank for useful discussions K. Kifonidis, E. Muller, A. Poludnenko. The work of F.F. was supported by CNES (the French Space Agency) and was carried out at Service d'Astrophysique, CEA/Saclay and partially at the Center for Particle Astrophysics and Cosmology (APC) in Paris. The authors are also grateful to the anonymous referee for the useful suggestions.
References
 Amato, E, & Blasi, P. 2009, MNRAS, 392, 1591 [NASA ADS] [CrossRef] [Google Scholar]
 Arnett, W. D. 1988, ApJ, 331, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Bamba, A., Yamazaki, R., Yoshida, T., et al. 2005, ApJ, 621, 793 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 2004, MNRAS, 353, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Bietenholz, M. F., Bartel, N., & Rupen, M. P. 2001, ApJ, 557, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., & Ellison, D. C. 2001 ApJ, 560, 244 [Google Scholar]
 CassamChenaï, G. Hughes, J. P., Reynoso, E. M., et al. 2008, ApJ, 680, 1180 [NASA ADS] [CrossRef] [Google Scholar]
 Chevalier, R. A. 1982, ApJ, 258, 790 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Chevalier, R. A. 1983, ApJ, 272, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Chevalier, R. A., Blondin, J. M., & Emmering, R. T. 1992, ApJ, 392, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Decourchelle, A., Ellison, D., & Ballet, J. 2000, A&A, 543, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Dumbser, M., Moschetta, J.M., & Gressier, J. 2004, J. Comput. Phys., 197, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Dwarkadas, V. V. 2000, ApJ, 541, 418 [NASA ADS] [CrossRef] [Google Scholar]
 Dwarkadas, V. V., & Chevalier, R. A. 1998, ApJ, 497, 807 [NASA ADS] [CrossRef] [Google Scholar]
 Ellison, D. C., Decourchelle, A., & Ballet, J. 2005, A&A, 429, 569 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fesen, R. A., & Gunderson, K. S. 1996, ApJ, 470, 967 [NASA ADS] [CrossRef] [Google Scholar]
 Fermi, E., & von Neumann, J. 1953, Technical Report No. AECU2979, Los Alamos Scientific Laboratory [Google Scholar]
 Giacalone, J., & Jokipii, J. R. 2007, ApJ, 663, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Gnedin, N. Y. 1995, ApJ, 97, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Gnedin, N. Y., & Bertschinger, E. 1996, ApJ, 470, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Gull, S. F. 1973, MNRAS, 161, 47 [NASA ADS] [Google Scholar]
 Gull, S. F. 1975, MNRAS, 171, 263 [NASA ADS] [Google Scholar]
 Hoyle, F., & Fowler, W. A. 1960, ApJ, 132, 565 [NASA ADS] [CrossRef] [Google Scholar]
 Hughes, J. P., Rakowski, C. E., Burrows, D. N., et al. 2000, ApJ, 528, L109 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Johnsen, E., & Colonius, T. 2006, J. Comput. Phys., 219, 715 [NASA ADS] [CrossRef] [Google Scholar]
 Jun, B.I., & Norman, M. L. 1996, ApJ, 472, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Kane, J., Arnett, D., Remington, B. A., et al. 2000, ApJ, 528, 989 [NASA ADS] [CrossRef] [Google Scholar]
 Kaps, P., & Rentrop, P. 1979, Numer. Math., 33, 55 [CrossRef] [Google Scholar]
 Kifonidis, K., Plewa, T., Janka, H.Th., et al. 2003, A&A, 408, 621 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lagage, P. O., & Cesarsky, C. J. 1983, A&A, 125, 249 [NASA ADS] [Google Scholar]
 Laming, J. M., & Hwang, U. 2003, ApJ, 597, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Lazendic, J. S., Slane, P. O., Gaensler, B. M., et al. 2003, ApJ, 593, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Martel, H., & Shapiro, P. R. 1998, MNRAS, 297, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Nishikawa, H., & Kitamura, K. 2008, J. Comput. Phys., 227, 2560 [NASA ADS] [CrossRef] [Google Scholar]
 Pandolfi, M., & D'Ambrosio, D. 2001, J. Comput. Phys., 166, 271 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Peery, K. M., & Imlay, S. T. 1988, Blunt Body Flow Simulations, AIAA Paper 882924 [Google Scholar]
 Poludnenko, A. Y., & Khokhlov, A. M. 2007, J. Comput. Phys., 220, 678 [NASA ADS] [CrossRef] [Google Scholar]
 Quirk, J. J. 1994, Int. J. Num. Meth. Fluids, 18, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Quirk, J. J. 1997, in Upwind and HighResolution Schemes, ed. M. Yousuff Hussaini, B. van Leer, & J. van Rosendale (Berlin: Springer), 550 [Google Scholar]
 Rothenflug, R., Magne, B., Chieze, J. P., et al. 1984, A&A, 291, 271 [NASA ADS] [Google Scholar]
 Shirkey, R. C. 1978, ApJ, 224, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Shyue, K. M. 1998, J. Comput. Phys., 142, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, A., Peacock, A., Arnaud, M., et al. 1989, ApJ, 347, 925 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Gardiner, T. A., Teuben, P., et al. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Synge, J. L. 1957 The Relativistic Gas (Amsterdam: North Holland Publishing Co.) [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Toro, E. F. 1999, Riemann Solvers and Numerical Methods for Fluid Dynamics (SpringerVerlag) [Google Scholar]
 Truelove, J. K., & McKee, C. F. 1999, ApJS, 120, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. 2008a, ApJ, 689, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Vink, J. 2008b, Proc. of the 4th International Meeting on High Energy GammaRay Astronomy, AIP Conf. Proc., 1085, 169 [NASA ADS] [Google Scholar]
 Vink, J., & Laming, J. M. 2003, ApJ, 584, 758 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, C.Y., & Chevalier, JR. A. 2001, ApJ, 549, 1119 [NASA ADS] [CrossRef] [Google Scholar]
 Warren, J. S., Hughes, J. P., Badenes, C., et al. 2005, ApJ, 634, 376 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, K., & Li, Z. 2001, Int. J. Num. Meth. Fluids, 37, 1 [CrossRef] [Google Scholar]
 Young, Y.N., Tufo, H., Dubey, A., et al. 2001, J. Fluid Mech., 447, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Zel'dovich, Ya. B., & Raizer, Yu. P. 2002, Physics of shock waves and hightemperature hydrodynamic phenomena (New York: Dover Publications) [Google Scholar]
All Figures
Figure 1: Angleaveraged radial density profile at t=50 yr for , or (n,s) = (7,0), and with 13 levels of refinement (128^{3} in yellow, (128256)^{3} in green, (128512)^{3} in black). For comparison, the selfsimilar profile is shown in red. The density profiles are normalized to their values at the forward shock. From this diagram, it can be seen that the localization of the forward shock and its compression ratio are not strongly dependent on the maximal refinement level. 

Open with DEXTER  
In the text 
Figure 2: Angleaveraged radial density profile at various ages for , or (n,s) = (7,0), and with 1 and 3 levels of refinement (128^{3} in yellow, 128^{3}512^{3} in black). For comparison, the red line shows the selfsimilar profiles normalized to their forward shock values, Chevalier profiles for the first five panels and Sedov for the last four panels. The departure from the ejectadominated phase, occurring when the reverse shock encounters the plateau outer radius , and the transition to SedovTaylor phase are shown. The inward motion of the reverse shock, initially only in the comoving frame of the contact discontinuity and later on in the laboratory frame as well, is also depicted. 

Open with DEXTER  
In the text 
Figure 3: Temporal evolution of the forward shock and reverse shock radii for , or (n,s) = (7,0), and . Our 3D numerical solution (solid) matches the 1D semianalytical solution (Truelove & McKee 1999), dashed, and, in the asymptotic regime for , matches the selfsimilar solution in the ejectadominated phase (Chevalier 1982), dotted. See Sect. 5.1 for the units and ( pc and yr for the chosen parameters). The spherical symmetry of the reverse shock is strongly modified by the RT instability; here the innermost value of the radius of the reverse shock has been chosen (cf. Sect. 5.1). The ejectadominated solution (dotted) is prolonged indefinitely to emphasize the inward motion of the reverse shock. 

Open with DEXTER  
In the text 
Figure 4: Snapshots at two distinct ages (t=100 yr at left; t=500 yr at right) for , or (n,s) = (7,0) and . In the upper panel, the mass fraction , defined in Eq. (4), is shown in a comoving slice plane parallel to the coordinate plane XY at height z=0.1 a(t). The mixture of the two media, namely ejecta shocked and ISM shocked, as well as the deformation of the contact discontinuity surface is shown. In the middle panel, matter density is shown in the same plane as in the upper panel. Density values are coded according to the color bars given for each frame. Note the change in the radial scale. In the bottom panel, frequency emissivity due to bremsstrahlung integrated along the line of sight towards the reader is shown, assuming a mean molecular weight , corresponding to the case of fully ionized hydrogen. 

Open with DEXTER  
In the text 
Figure 5: The threedimensional density in the full simulation box of side L = 4.2 pc is shown at t=500 yr for , or (n,s) = (7,0) and . The black line indicates the height of the planar slice at z = 0.1 L used in Fig. 4, upper and middle panels. 

Open with DEXTER  
In the text 
Figure 6: Snapshots showing streamlines of velocity field superposed to matter density for , or (n,s) = (7,0) and in the frame of contact discontinuity. The planar slice is in all panels z=0. From left to right, in above panel (t=100 yr) the inflow motion of the ISM and the outflow motion of the freelyexpanding ejecta are clearly evident. In the other panels (t=500 yr, t=1000 yr, t=2000 yr), the turbulent flow motion due to RT instabilities is shown. Note that the scale in color bar is log in amu cm^{3} in the first panel and linear in all the other panels. 

Open with DEXTER  
In the text 
Figure 7: The protruding RT structures measured in our simulations at different times are quite well reproduced by the exact solution H(t) with (see Eq. (33)). The edge of the RT structures is defined by . The earlytime exponential growth can hardly be verified with the available spatial resolution. A second order polynomial fit is also shown for comparison. 

Open with DEXTER  
In the text 
Figure 8: Angleaveraged radial density profile shown at distinct ages for , or (n,s) = (7,0) in the hybrid case in black, is compared with the corresponding initial profile, in red. 

Open with DEXTER  
In the text 
Figure 9: Matter density in a comoving planar slice parallel to the coordinate plane XY at height z=0.1 a(t) (t=100 yr at left; t=500 yr at right) for , or n=7, s=0, and . Density values are coded according to the color bars given for each frame. Note the change in the radial scale. In this case, . Since , the computational grid expands faster than the SNR and the volume fraction occupied in the simulation box by the SNR decreases. It is evident that even a small drift of the computational grid from the shock velocity erases the effects of the numerical instability possibly emerging behind the shocks. 

Open with DEXTER  
In the text 
Figure 10: Angleaveraged radial density profiles between the reverse shock and the forward shock at time t=12 yr = 0.0063 , with , or ( n,s)=(7,0), and . The AMR is performed over three levels of refinement. The radial coordinate is normalized to the forward shock position corresponding to ; the density is normalized to the value at the shock. Upper panel. The profile in the case of the exactly comoving grid is compared with the profiles having a quasicomoving grid . The initial profile, at t_{0} = 10 yr, is also depicted for comparison. It is evident that even small values of produce oscillations and an error of the order of in the compression factor for . Outside the interaction region, the profiles are superimposed. Lower panel. The curves without AMR and with AMR are compared with the initial profile at t_{0}=10 yr. In this case, . The independence of the error in the compression factor at the two shocks from the space resolution is clear. 

Open with DEXTER  
In the text 
Figure 11: Angleaveraged radial profile of the mass fraction of the ejecta f shown at distinct ages for , or (n,s) = (7,0); the cases with uniform , and the hybrid case are compared. The radius is in units of the corresponding radius of the contact discontinuity. The vertical lines on the right of the panels indicate the corresponding positions of the forward shock, which at this age still meets the selfsimilar expansion (see Fig. 3 for the uniform case). The vertical line on the left of the panels indicates the position of the reverse shock. 

Open with DEXTER  
In the text 
Figure 12: The protrusion of RT structures measured in our simulations in units of the corresponding selfsimilar position of the contact discontinuity, whose surface is directly modified by the instabilities; the uniform and the hybrid cases are compared ( ). The cosmic rays affect the elongation of the instabilities from the very beginning by a factor of the order of . 

Open with DEXTER  
In the text 
Copyright ESO 2010