LETTER TO THE EDITOR
3D simulations of supernova remnants evolution including non-linear particle acceleration
G. Ferrand1 - A. Decourchelle1 - J. Ballet1 - R. Teyssier1,2 - F. Fraschetti3,4
1 - Laboratoire AIM (CEA/Irfu, CNRS/INSU, Université Paris VII), CEA Saclay, Bât. 709, 91191 Gif-sur-Yvette, France
2 - Institute for Theoretical Physics, University of Zürich, 8057 Zürich, Switzerland
3 - LUTh, Observatoire de Paris, CNRS-UMR8102 et Université Paris VII, 92195 Meudon Cedex, France
4 - Lunar and Planetary Lab & Dept. of Physics, University of Arizona, Tucson, AZ, 85721, USA
Received 13 November 2009 / Accepted 19 December 2009
If a sizeable fraction of the energy of supernova remnant shocks is channeled into energetic particles (commonly identified with Galactic cosmic rays), then the morphological evolution of the remnants must be distinctly modified. Evidence of such modifications has been recently obtained with the Chandra and XMM-Newton X-ray satellites. To investigate these effects, we coupled a semi-analytical kinetic model of shock acceleration with a 3D hydrodynamic code (by means of an effective adiabatic index). This enables us to study the time-dependent compression of the region between the forward and reverse shocks due to the back reaction of accelerated particles, concomitantly with the development of the Rayleigh-Taylor hydrodynamic instability at the contact discontinuity. Density profiles depend critically on the injection level of particles: for modifications are weak and progressive, for modifications are strong and immediate. Nevertheless, the extension of the Rayleigh-Taylor unstable region does not depend on the injection rate. A first comparison of our simulations with observations of Tycho's remnant strengthens the case for efficient acceleration of protons at the forward shock.
Key words: ISM: supernova remnants - instabilities - cosmic rays - acceleration of particles - methods: numerical
Supernova remnants (SNRs) are believed to be the major contributors to Galactic cosmic rays. But although acceleration of electrons in these objects is well established, direct evidence of the acceleration of protons remains elusive: recent detections of non-thermal gamma rays do not yet allow unambiguously distinguishing leptonic and hadronic contributions (see e.g. Gabici 2008, and references therein).
To assess particle acceleration in SNRs, a promising alternative approach consists in diagnosing the impact of energetic particles on the SNR dynamics. Indeed, if SNRs are efficient accelerators, as claimed, then energetic nuclei must make a sizable impact on their evolution (see e.g. Jones & Ellison 1991; Malkov & Drury 2001). Such a study is now possible thanks to the performance of modern X-ray observatories such as Chandra and XMM-Newton, which allow the structure of young Galactic SNRs to be probed in great detail. Warren et al. (2005) have been able to determine the positions of the three waves (forward shock, contact discontinuity, reverse shock) in Tycho's SNR with rather good accuracy, and have found that it does not match any pure hydrodynamic model. The same effect has been reported in SN 1006 (Cassam-Chenaï et al. 2008; Miceli et al. 2009) and in Cas A (Patnaude & Fesen 2009). This is evidence that not all the kinetic energy from the explosion is converted into heat downstream of the shock, but that a sizeable part is channeled elsewhere - probably into energetic particles.
This effect has been studied by Decourchelle et al. (2000), using 1D self-similar simulations coupled with a simple model of non-linear acceleration (from Berezhko & Ellison 1999). They have shown how the shocked region shrinks in the case of efficient acceleration of particles at the shocks: as the injection fraction rises the shocks get closer to the contact discontinuity. These results have been confirmed and extended by 1D hydrodynamic simulations of radially symmetric SNRs (see Ellison et al. 2007, and references therein).
But the contact discontinuity between the shocked ISM and the shocked ejecta is known to be hydrodynamically unstable: this interface is subject to the Rayleigh-Taylor instability, as the ejecta are being decelerated by the ambient medium of lower density. Thus the morphological study of a SNR requires a 3D (or at least 2D) modelling (see Chevalier et al. 1992; Dwarkadas 2000; Wang & Chevalier 2001, and references therein). As the instability develops the contact discontinuity is profoundly modified: the ejecta protrude inside the shocked ISM, forming characteristic finger-like structures. These features are particularly apparent in Tycho's SNR, where the ejecta exhibit a fleecy aspect in both X-rays (Warren et al. 2005) and in radio (Velazquez et al. 1998). Thus, to diagnose the back reaction of energetic particles, one needs to take hydrodynamic instabilities into account. To assess their impact, Blondin & Ellison (2001) have made 2D and 3D hydrodynamic simulations of a slice of SNR, mimicking the presence of energetic particles by lowering the adiabatic index of the fluid (uniformly in space and time).
In this work, for the first time we combine all these previous approaches; that is, we make full 3D simulations of an SNR evolution including a space- and time-dependent model of acceleration and back reaction of particles, to be able to interpret the X-ray observations.
As we are interested in the time-dependent and non-linear interplay between the SNR and the energetic particles it accelerates, we resort to numerical simulations. To compute the remnant evolution, we used the code RAMSES (Teyssier 2002). This 3D Cartesian Eulerian code includes a second-order hydrodynamic solver, and implements a tree-based structure allowing for versatile adaptive mesh refinement (AMR).
Although RAMSES has been already extensively tested and used, we had to adapt it to the study of supernova remnants (Fraschetti et al. 2009). The main point is that we use a grid that is comoving with the contact discontinuity; that is, we work in an expanding frame. Because this frame is non-inertial, we have to modify the Euler equations. Although it is computationally very interesting to factor out the global expansion of the remnant in this way, we have to face the numerical instabilities associated to quasi-stationary shock waves. Still, we can accurately follow the SNR evolution as shown in Fig. 1. In the unmodified case, the position and speed of the shocks exactly coincide with analytical predictions by Truelove & McKee (1999). In the simulations presented here, we actually simulate one eighth of a sphere, assuming symmetry.
Evolution of the shocks.
|Open with DEXTER|
To compute acceleration of particles by shocks, we used the semi-analytical model of Blasi (2002,2004), a non-linear model of diffusive shock acceleration (DSA) that takes the back reaction of accelerated particles on the fluid structure into account. This model solves the particle spectrum f(p) and the fluid velocity profile U(p) jointly as functions of the momentum p of particles. It includes the escape of particles with the highest energy upstream of the shock, which carry energy away from the accelerator. We also include the effect of Alfvén wave heating in the precursor, which limits shock modifications, but we do not include magnetic field amplification (Amato & Blasi 2006).
As inputs, the model requires (i) information on the shock (its speed uS and Mach number MS are provided by the hydrodynamic code, averaged over the remnant surface); (ii) an injection recipe (we assume that a fraction of the particles crossing the shock enter the accelerator, at injection momentum equal to times the mean downstream thermal momentum); and (iii) a cutoff recipe (we limit the maximum momentum according to the age and the size of the remnant, assuming Bohm diffusion, i.e. a diffusion coefficient with ). We consider particle acceleration only at the forward shock, as there is less theoretical and observational evidence of efficient acceleration at the reverse shock (see Ellison et al. 2005 for a discussion of this issue).
Amongst the outputs, the acceleration model provides the total shock compression ratio . To couple the hydrodynamic evolution of the remnant with particle acceleration, we vary the adiabatic index of the fluid as done in Ellison et al. (2004): at each time-step, we compute the index , which would have produced the same ratio , and affect it in RAMSES to the cells located just upstream of the shock front. Then is advected inside the shocked region, so that it remains constant in each fluid element, which thus remembers modifications induced by particle acceleration at the time it was shocked. Ellison et al. (2004) have shown that there is good agreement between such a pseudo-fluid approach and two-fluid calculations in 1D.
We initialise our simulations at a young age (here 10 years), using self-similar profiles from Chevalier (1983), including the pressure of accelerated particles (as computed from the acceleration model presented in the previous section). Assuming that both the ejecta (but for a central uniform core) and the ambient medium have power-law density profiles (of indices respectively n and s), hydrodynamic profiles are obtained by integration of ordinary differential equations.
Here we are interested in a Tycho-like SNR, that is a supernova of type Ia, referring to 1051 erg of kinetic energy in 1.4 solar masses, with an n=7 power-law distribution. We assume a uniform (s=0) and tenuous ( ) ambient medium.
Finally, we mention that Rayleigh-Taylor instabilities are not explicitly seeded in the simulation, but are spontaneously triggered at the contact discontinuity by numerical perturbations seeded by the grid.
Evolution of the remnant structure.
|Open with DEXTER|
3.1 Remnant evolution The temporal evolution of some key parameters is shown in Fig. 2. The acceleration model depicted in Sect. 2.2 provides the compression ratios plotted on top of the figure. For ease of interpretation, we ran multiple simulations with different fixed injection rates . If (in dark blue), the system is almost in the linear (test-particle) regime, and there is a single strong shock of compression ratio r=4. As we raise to 10-3 (as colour gets warmer), it progressively enters the non-linear (modified) regime, and the shock discontinuity is reduced to a subshock of compression between 3 and 4, whereas the overall compression ratio increases to more than 10. (Contrary to ordinary hydrodynamics, always significantly depends on the shock Mach number, even when the latter is very high, see e.g. Blasi 2002). The corresponding effective adiabatic index is plotted in the middle of the figure. As expected, it decreases from 5/3 (thermal fluid) to 4/3 (relativistic fluid) or even below (as particles escape). In initially poorly efficient accelerators (low ), shock modifications are small, but steadily increase over time. On the other hand, in very efficient accelerators (high ), back reaction effects are very strong even at a very early stage, but then decrease over time. In all cases, we see that parameters evolve rather slowly after the first hundred years. These results are similar to those obtained by Decourchelle et al. (2000) in 1D, assuming self-similarity of the hydrodynamic profiles and using a different acceleration scheme, which cross-validates the models.
Maps of the density.
|Open with DEXTER|
Finally, the evolution of the relative positions of the waves is shown at the bottom of the figure.
Two main comments are in order.
First, the region affected by the Rayleigh-Taylor instability
does not seem to be affected by the acceleration of particles.
Second, the forward shock can get very close to this turbulent region
if the injection level is high enough (above
- the reverse shock also reaches the contact discontinuity,
but this is caused by the development of the Rayleigh-Taylor instability,
independent of acceleration efficiency.
These observations agree with the results of Blondin & Ellison (2001) and of Fraschetti et al. (2009),
obtained with no space- and time-dependent model of acceleration.
Our simulations also show that the average distance between the forward
shock and the contact discontinuity does not depend much on time in the
case of efficient acceleration. (Rayleigh-Taylor fingers grow steadily
in time, but the forward shock moves away as back reaction effects are
On the left of the figure, we show slices of the density, which highlights the remnant's structure. The main effect of particle back reaction is apparent, namely that the shocked region shrinks. For this level of injection, this causes density to rise by a factor typically of two downstream of the forward shock. Nevertheless, the size of the region perturbed by the Rayleigh-Taylor instability, the development of which depends on the density contrast, is basically unchanged.
On the right of the figure we show projected maps of the square of the density of shocked ejecta, which roughly approximates the intensity of their thermal X-ray emission. The interior of the remnant is filled with a filamentary texture, structured over scales compatible with the picture of Warren et al. (2005). Even in projection, the unstable region is still clearly visible as a bright annulus. But here again, cases with and without efficient acceleration cannot be distinguished from the structure of the shocked ejecta alone.
Although the detailed modelling of a specific object is beyond the scope of this Letter, we can show the interest of our approach by comparing our numerical results with the observational results of Warren et al. (2005) regarding Tycho's remnant (see also Völk et al. 2008, for a study of this object using a 1D non-linear kinetic model). Their main finding is that the forward shock is very close to the contact discontinuity, with a mean FS:CD ratio of 1:0.93. Assuming a uniform medium (s=0, a reasonable hypothesis for Tycho), the theoretical ratio is 1:0.85 for an ejecta indice n=7 (assuming a self-similar evolution without particle acceleration). As n increases, the predicted ratio decreases, but is still only 1:0.89 for n=14. As pointed out by Dwarkadas & Chevalier (1998), an exponential profile may be more adapted to the early stages of an SN Ia; however, this would produce an even broader shocked region at the time considered (Cassam-Chenaï et al. 2008). On the other hand, such a high FS:CD ratio can be obtained readily with our code through the back reaction of accelerated particles. The procedure used by Warren et al. (2005)seems to extract the envelope of the ejecta, which can be obtained in our simulations by setting some threshold on the ejecta fraction - we found the value of 10% used here to give comparable results. Then the simulations can reproduce the observations provided that is of the order of (Völk et al. 2008 obtained ). Although we have not conducted a complete parametric study yet, we believe that uncertainties in SNR parameters or in observed ratios could be accounted for by varying the injection level - and reciprocally, a good knowledge of the object will allow this crucial parameter to be constrained.
Warren et al. (2005) also found that that the reverse shock is deep inside the ejecta, with a mean FS:RS ratio of 1:0.71, which is quite puzzling. For n=7 and s=0 the predicted ratio(without particle acceleration) is 1:0.79. Projection effects might cause underestimates of the radius of the reverse shock, but according to our maps this effect is too small to explain the discrepancy. Considering more realistic profiles from explosion models might help understanding the situation. In any case, this observation does not favour efficient acceleration at the reverse shock, which would shrink the whole remnant's structure even more. However, other potentially important effects (like the composition and temperature of the ejecta) have to be included in our model before drawing firm conclusions.
We have presented a new code that couples a 3D hydrodynamic description of a supernova remnant, allowing consideration of hydrodynamic instabilities such as the Rayleigh-Taylor instability, with a realistic kinetic model of non-linear acceleration at the blast wave, allowing evaluation of the efficiency of particle acceleration. We are thus now able for the first time to simulate the morphology of SNRs undergoing efficient DSA without limiting assumptions on its spatial structure or temporal evolution.
Our first results confirm the most notable previous findings regarding particle back reaction on the SNR morphology: (i) the shock structure is all the more compact since acceleration is efficient, which provides a clear observational diagnostic; and (ii) the development of the Rayleigh-Taylor instability is not significantly affected by acceleration at the forward shock, but it has to be taken into account when interpreting observations.
Regarding the case of Tycho's remnant, comparison of our simulations with X-ray observations strengthens the case for efficient acceleration of protons at the forward shock.Acknowledgements
This work has been partially funded by the ACCELRSN project ANR-07-JCJC-0008. The numerical simulations were performed with the DAPHPC cluster at CEA/Irfu. The authors thank the anonymous referee for his/her useful comments.
- Amato, E., & Blasi, P. 2006, MNRAS, 371, 1251 [NASA ADS] [CrossRef] [Google Scholar]
- Berezhko, E. G., & Ellison, D. C. 1999, ApJ, 526, 385 [NASA ADS] [CrossRef] [Google Scholar]
- Blasi, P. 2002, Astropart. Phys., 16, 429 [NASA ADS] [CrossRef] [Google Scholar]
- Blasi, P. 2004, Astropart. Phys., 21, 45 [NASA ADS] [CrossRef] [Google Scholar]
- Blasi, P., Gabici, S., & Vannoni, G. 2005, MNRAS, 361, 907 [NASA ADS] [CrossRef] [Google Scholar]
- Blondin, J. M., & Ellison, D. C. 2001, ApJ, 560, 244 [NASA ADS] [CrossRef] [Google Scholar]
- Cassam-Chenaï, G., Hughes, J. P., Reynoso, E. M., Badenes, C., & Moffett, D. 2008, ApJ, 680, 1180 [NASA ADS] [CrossRef] [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. C., & Ballet, J. 2000, ApJ, 543, L57 [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. 2004, A&A, 413, 189 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ellison, D. C., Decourchelle, A., & Ballet, J. 2005, A&A, 429, 569 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Ellison, D. C., Patnaude, D. J., Slane, P., Blasi, P., & Gabici, S. 2007, ApJ, 661, 879 [NASA ADS] [CrossRef] [Google Scholar]
- Fraschetti, F., Teyssier, R., Ballet, J., et al. 2009, A&A, submitted [Google Scholar]
- Gabici, S. 2008, in Proceedings of the XXI European Cosmic Ray Symposium [Google Scholar]
- Jones, F. C., & Ellison, D. C. 1991, Space Sci. Rev., 58, 259 [NASA ADS] [CrossRef] [Google Scholar]
- Malkov, M. A., & Drury, L. O. 2001, Rep. Progr. Phys., 64, 429 [NASA ADS] [CrossRef] [Google Scholar]
- Miceli, M., Bocchino, F., Lakubovskyi, D., et al. 2009, A&A, 501, 239 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Patnaude, D. J., & Fesen, R. A. 2009, ApJ, 697, 535 [NASA ADS] [CrossRef] [Google Scholar]
- Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Truelove, J. K., & McKee, C. F. 1999, ApJS, 120, 299 [NASA ADS] [CrossRef] [Google Scholar]
- Velazquez, P. F., Gomez, D. O., Dubner, G. M., de Castro, G. G., & Costa, A. 1998, A&A, 334, 1060 [NASA ADS] [Google Scholar]
- Völk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2008, A&A, 483, 529 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Wang, C.-Y., & Chevalier, R. 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]
Evolution of the shocks.
|Open with DEXTER|
|In the text|
Evolution of the remnant structure.
|Open with DEXTER|
|In the text|
Maps of the density.
|Open with DEXTER|
|In the text|
Copyright ESO 2010