Numerical simulations of diffusive shock acceleration in SNRs
A new numerical model of the nonlinear diffusive shock acceleration is presented. It is used for modeling of particle acceleration in supernova remnants. The model contains coupled spherically symmetric hydrodynamic equations and the transport equations for energetic protons, ions and electrons. The forward and reverse shocks are included in the consideration. The spectra of cosmic rays released into interstellar medium from a supernova remnant are determined. The role of the reverse shock in the production of CR ions and positrons is discussed.
keywords:cosmic rays, acceleration, supernova remnants
The diffusive shock acceleration (DSA) process (1); (2); (3); (4) is considered as the principal mechanism for production of galactic cosmic rays (CR) in supernova remnants (SNRs). A significant theoretical progress in the investigation of this mechanism was achieved (see e.g. Malkov & Drury (5) for a review). However only during the last decade the excellent results of X-ray and gamma-ray astronomy supplied the observational evidence of the presence of multi-TeV energetic particles in these objects (see e.g. Aharonian et al. (6)).
Two shocks are produced by super-sonically moving supernova ejecta after an explosion. A forward shock propagates in the circumstellar medium while a reverse shock propagates in the gas of ejecta. Generally it is believed that some part of thermal particles is injected at the shock fronts into acceleration.
In this paper we present a new numerical model of nonlinear diffusive shock acceleration. This model is a natural development of the existing models (7); (8). The solution of spherically symmetric hydrodynamic equations is combined with the energetic particle transport and acceleration at the forward and reverse shocks of a supernova remnant. Nonlinear response of energetic particles via their pressure gradient results in the self-regulation of acceleration efficiency.
Our previous studies that used this model dealt with the CR spectra produced by SNRs (9); (10). The input of the reverse shock was not taken into account there while the acceleration by both shocks was considered for modeling of non-thermal electromagnetic emission from the SNR RX J1713.7-3946 (11).
The paper is organized as follows. The short description of the model is given in Sect. 2. The results of modeling of evolution of a supernova remnant in the interstellar medium are presented in Sect. 3. Sect. 4 contains the discussion of our results. Our conclusions are given in the last Section. The numerical code is described in Appendix.
2 Nonlinear model of diffusive shock acceleration
Hydrodynamical equations for the gas density , gas velocity , gas pressure , and the equation for isotropic part of the CR proton momentum distribution in the spherically symmetrical case are given by
Here is the CR pressure, is the advective velocity of CRs,
is the adiabatic index of the gas, and is the CR diffusion coefficient.
It was assumed that the diffusive streaming of CRs results in the generation of magnetohydrodynamic (MHD)
waves. CR particles are scattered by these waves. That is why the CR advective velocity
may differ from the gas velocity . Damping of these waves results in an additional gas heating. It is
described by the last term in Eq. (3). Two last terms in Eq. (4)
correspond to the injection of thermal protons with momenta
, and mass at the fronts of the forward and
reverse shocks at and
The equation for ions is similar to Eq. (4). For ions with the mass and the mass number it is convenient to use the momentum per nucleon and the normalization of the ion spectra to the baryonic number density. Then the number density of ions is . The ion pressure is also taken into account in the cosmic ray pressure .
We shall neglect the pressure of energetic electrons. In other words they are treated as the test particles. The evolution of electron distribution is described by equation similar to Eq. (4) with additional terms describing synchrotron and inverse Compton (IC) losses.
CR diffusion is determined by magnetic inhomogeneities. Strong streaming of accelerated particles changes medium properties in the shock vicinity. CR streaming instability results in the high level of MHD turbulence (2) and even in the amplification of magnetic field in young SNRs (12). Due to this effect the maximum energy of accelerated particles may be higher in comparison with previous estimates of Lagage and Cesarsky (14).
According to the numerical modeling of this non-resonant instability, the magnetic field is amplified by the flux of run-away highest energy particles in the relatively broad region upstream of the shock (13). Magnetic energy density is a small fraction () of the energy density of accelerated particles. This amplified almost isotropic magnetic field can be considered as a large-scale magnetic field for lower energy particles which are concentrated in the narrow region upstream of the shock. Resonant streaming instability of these particles produces MHD waves propagating in the direction opposite to the CR gradient. Strong nonlinear damping of these waves results in the gas heating (see the last term in Eq.(3)). CR gradient is negative upstream of the forward shock and MHD waves propagate in the positive direction. The situation changes downstream of the forward shock where CR gradient is as a rule positive and MHD waves propagate in the negative direction. This effect is mostly pronounced downstream of the forward shock of SNR because the magnetic field is additionally amplified by the shock compression and the Alfvén velocity can be comparable with the gas velocity in the shock frame . As for CR diffusion coefficient, it is probably close to the Bohm value , where is the electric charge of particles.
We apply a finite-difference method to solve Eqs (1-4) numerically upstream and downstream of the forward and reverse shock (see Appendix). A non-uniform numerical grid upstream of the shocks at and allows to resolve small scales of hydrodynamical quantities appearing due to the pressure gradient of low-energy CRs. The gases compressed at the forward and reverse shocks are separated by a contact discontinuity at that is situated between the shocks. An explicit conservative TVD scheme (15) for hydrodynamical equations (1-3) and uniform spatial grid were used between the shocks.
Because of synchrotron losses the spacial scale of the high-energy electrons can be rather small downstream of the shocks. That is why we use a non-uniform numerical spacial grid for accelerated electrons downstream of the shocks.
The magnetic field plays no dynamical role in the model. Since we have not perform the modeling of the amplification and transport of magnetic field here, its coordinate dependence should be specified for determination of cosmic ray diffusion and for calculation of synchrotron emission and losses. We shall assume below that the coordinate dependencies of the magnetic field and the gas density coincide upstream and downstream of the forward shock:
Here is the gas density and is the Alfvén velocity of the circumstellar medium. The parameter determines the value of the amplified magnetic field strength. For low shock velocities the magnetic field is not amplified.
The magnetic energy is about 3.5 percent of the dynamical pressure according to estimates from the width of X-ray filaments in young SNRs (16). This number and characteristic compression ratio of a modified SNR shock correspond to . Since the plasma density decreases towards the contact discontinuity downstream of the forward shock, the same is true for the magnetic field strength according to Eq. (5). This seems reasonable because of a possible magnetic dissipation in this region.
Situation is different downstream of the reverse shock at . The plasma flow is as a rule strongly influenced by the Rayleigh-Taylor instability that occurs in the vicinity of the contact discontinuity and results in the generation of MHD turbulence in this region. We shall assume that the magnetic field does not depend on radius downstream of the reverse shock while the dependence in the upstream region is described by the equation similar to Eq. (5):
Here is the radius where the ejecta density has a minimum and equals . This radius is generally close to the reverse shock radius and is equal to it if the reverse shock is not modified by the cosmic ray pressure.
CR advective velocity may differ from the gas velocity on the value of the radial component of the Alfvén velocity calculated in the isotropic random magnetic field: . Here the factor describes the possible deviation of the cosmic ray drift velocity from the gas velocity. The similar expression for the cosmic ray drift velocity is used upstream of the reverse shock at . We shall use values and upstream of the forward and reverse shocks respectively, where Alfvén waves are generated by the cosmic ray streaming instability and propagate in the corresponding directions. The damping of these waves heats the gas upstream of the shocks (17) and limits the total compression ratios by a number close to 6. This Alfvén heating produce a more efficient limitation of the shock modification in comparison with the dynamical effects of the magnetic field considered by Caprioli et al. (18) which are neglected in our study. In the downstream region of the forward and reverse shock at we set and therefore in the major part of this paper except the consideration of the Alfvén drift effects downstream of the shocks at the end of Sect.3 and in Fig.10.
We shall use the following diffusion coefficient
Here the parameter depends on the type of nonlinear wave damping which is relevant only for low velocity shocks when the magnetic field is not amplified. The parameter describes the possible deviations of diffusion coefficient from the Bohm value for high-velocity shocks. Since the highest energy particles are scattered by small-scale magnetic fields, their diffusion is faster than the Bohm diffusion (13). The same is true for smaller energy particles because they can be resonantly scattered only by a fraction of the magnetic spectrum. We shall use the value throughout the paper.
We shall use the value of parameter below. It corresponds to the nonlinear wave damping of the weak turbulence theory. Note that a stronger Kolmogorov-type nonlinear damping used by Ptuskin & Zirakashvili (19) for estimate of maximum energy in SNRs corresponds to .
In real situation the level of MHD turbulence drops with distance upstream of the shock and diffusion may be even faster there. This is qualitatively taken into account by the exponential factors in Eq. (7). The characteristic diffusive scale of highest energy particles is a small fraction of the shock radius (13) and is determined by the generation and transport of MHD turbulence in the upstream region (20); (21). The value is determined by ratio of diffusion coefficient in the circumstellar medium and diffusion coefficient in the vicinity of the shock. The MHD turbulence is amplified exponentially in time before the shock arrival from the background level by cosmic ray streaming instability. The characteristic range of is (13). We shall use the value below.
It is believed that the supernova ejecta has some velocity distribution just after the supernova explosion (22)
Here the index characterizes the steep power-low part of this distribution. The radial distribution of ejecta density is described by the same expression with . The characteristic ejecta velocity can be expressed in terms of energy of supernova explosion and ejecta mass as
3 Numerical results
Figures (1)-(9) illustrate the numerical results that are obtained for the SNR shock propagating in the medium with a hydrogen number density cm, magnetic field strength G and temperature K. The fraction of helium nuclei was assumed. The gas of ejecta does not contain hydrogen in the case considered. We use the ejecta mass , the energy of explosion erg and the parameter of ejecta velocity distribution . The value of the parameter was assumed.
The initial forward shock velocity is km s. The injection efficiency is taken to be independent on time , and the injection momenta are , . Protons with a mass are injected at the forward shock while ions with mass number and charge number are injected at the reverse shock. This high injection efficiency results in the significant shock modification already at early stages of SNR expansion while the thermal sub-shock compression ratio is close to 2.5 during the simulation. This is in agreement with the radio-observations of young extragalactic SNRs (23) and with the modeling of collisionless shocks (24). Similar values of the injection efficiency were found in hybrid modeling (see e.g. (25)) and at the Earth bow shock in the solar wind (26).
As for the electron injection we assume a rather high injection energy of electrons MeV. This qualitatively corresponds to some models of suprathermal electron injection. Partially ionized ions accelerated at the shocks up to relativistic energies can produce multi-MeV electrons in the upstream region in the course of photo-ionization by Galactic optical and infrared radiation (27). MeV electrons and positrons are also present in the radioactive supernova ejecta while gamma-rays from Co decay in ejecta produce energetic electrons via Compton scattering in the circumstellar medium (28). These energetic particles may be additionally pre-accelerated via stochastic acceleration in the turbulent upstream regions of the shocks.
Below we assume that electrons are injected at the forward shock with efficiency while positrons are injected at the reverse shock with efficiency . These numbers are expected for the injection mechanisms mentioned above (see Sect.4 for details). Since electrons are considered as test particles our results can be easily rescaled for any other injection efficiency. The chosen injection rate at forward shock maintains the electron to proton ratio of the order of throughout the simulation while time-independent positron injection at the reverse shock results in positron to ion ratio increase from in the very beginning of SNR evolution up to at several thousand years shortly before the disappearance of the reverse shock.
In the real SNR the ions are also injected at the forward shock. We do not consider this process here in order to find the spectra of ions produced by the reverse shock. Because of the same reason we assumed the absence of hydrogen in the ejecta. Although this is the case for Ia/b/c and IIb supernovae it is not true for IIP supernovae. It is expected that the spectra of ions injected at the forward and reverse shock are similar to the spectra of protons. The production of CR ions at the forward shock was recently considered by Caprioli et al. (29).
The dependencies on time of the shock radii and , the forward and reverse shock velocities and , CR energy are shown in Fig.1. The calculations were performed until yr, when the value of the forward shock velocity drops down to km s and the forward shock radius is pc.
At early times of SNR evolution the distance between reverse and forward shocks is only 10 of the remnant radius. This is less than 23 thickness for automodel Chevalier-Nadezhin solution with (30) and should be attributed to a strong modification of both shocks by CR pressure. The reverse shock is strongly decelerated only when the forward shock sweeps the gas mass comparable with the ejecta mass at yr and when the transition to the Sedov phase begins.
Radial dependencies of physical quantities at yr are shown in Fig.2. The contact discontinuity between the ejecta and the interstellar gas is at pc. The reverse shock in the ejecta is situated at pc. At the Sedov stage the reverse shock moves in the negative direction and reach the center after seven thousand years after the supernova explosion. We stop the calculations in the region when the reverse shock radius .
It should be noted that our one-dimensional calculations cannot adequately describe the development of the Rayleigh-Taylor instability of the contact discontinuity. In the real situation the supernova ejecta and the circumstellar gas are mixed by turbulent motions in this region (see e.g. MHD modeling of Jun & Norman (31)).
Spectra of accelerated protons and electrons at yr are shown in Fig.3. At this epoch maximum energy of protons accelerated in this SNR is about 100 TeV, while higher energy particles have already left the remnant.
Radial dependencies of physical quantities at later epoch yr are shown in Fig.4. The reverse shock has reached the center of the remnant earlier. A weak reflected shock is clearly visible at pc.
The spectra of particles ar yr are shown in Fig.5.
Radial dependencies of physical quantities at the end of simulation at yr are shown in Fig.6. At this epoch the remnant is deeply in Sedov stage. The contact discontinuity is at pc.
The spectra of particles at yr are shown in Fig.7. The shock modification is not strong because the Alfvénic Mach number is close to 6 and the corresponding Alfvén heating upstream of the forward shock results in the lower compression ratio and acceleration efficiency. That is why the spectra of particles are steeper in comparison with ones at earlier epochs.
The spectra of particles produced during the whole evolution of the remnant are shown in Fig.8. They are obtained as the sum of the spectra integrated throughout simulation domain and of the time-integrated diffusive flux at the simulation boundary at . At yr the maximum energy of currently accelerated particles drops down to 100 GeV because of nonlinear damping. Higher energy particles have already left the remnant. Note that stronger Kolmogorov-type damping with will result even in lower energies of the order of 1 GeV. However we found that the spectra do not change in this case.
We found that the maximum energy of CR protons is somewhat less than eV. It is almost an order of magnitude lower in comparison with the results of Ptuskin et al. (10) where more optimistic assumptions and the spatially-uniform CR diffusion coefficient upstream of the forward shock were used (cf. Eq.(7)).
Note that the synchrotron losses of run-away electrons and positrons were taken into account in our modeling. The cut-off energy of the leptonic spectra TeV shown in Fig.8 is determined by the magnetic field strength G in the circumstellar medium and by the remnant age yr.
Evolution of non-thermal emission produced in the SNR at distance 4 kpc is shown in Fig.9. It is worth noting that a significant part of the IC emission is produced in the central region of the remnant at late epochs when the reverse shock have disappeared. This is because the magnetic field is rather weak in these regions.
Although only about of supernova energy is transferred to the particles accelerated at the reverse shock, they cannot be neglected. First of all the ejecta has absolutely different composition in comparison with interstellar medium where the forward shock propagates. Since the solar abundance corresponds to in the mass of heavy elements while the ejecta can contain up to 50 of heavy elements it is clear that the reverse shock will dominate in the production of heavy high-energy CR nuclei.
The relative input of the reverse shock at high energies is determined by the relative energetics of the reverse and forward shocks that is of the order of 1/10 for SNRs expanding in the uniform circumstellar medium.
According to our results, approximately of supernova energy is transferred to particles accelerated by forward shock. It is significantly higher than the estimate of 10-20 needed to maintain CR density in the Galaxy if the supernova rate is 1/30 yr. One of possibilities to resolve this contradiction is the assumption that CRs are accelerated only at small part of the forward shock surface. This can be due to the dependence of the proton and ion injection on the shock obliqueness (32). This effect is observed in SN 1006 and in the interplanetary medium.
This effect does not influence strongly the ion injection at the reverse shock. It is expected that the random magnetic field is the main component of the field in the expanded ejecta. This is because the magnetic field of ejecta originates from the magnetic field of the exploded star. Thus the random magnetic field strength of the red super giant progenitor of IIP supernovae is of the order of G similar to the magnetic field strength in the Sun interior while the regular field is of the order of 1 G. After a homogenous expansion from the initial stellar radius cm up to the radius cm of a young SNR the frozen-in magnetic field drops down to G. Although this value is significantly lower than the magnetic field in the interstellar medium it is strong enough for acceleration of particles up to 100 GeV in SNRs. A more realistic non-homogenous expansion will result in the stretching of the field in the radial direction. This can increase the magnetic field strength and ion injection efficiency at the reverse shock. Magnetic fields can be also amplified by non-resonant CR streaming instability suggested by Bell (12). If so the relative contribution of the reverse shock to the over-all CR spectrum increases.
The effect also depends on the type of the supernova explosion. It is known that acceleration at the reverse shock occurs in Cas A SNR (33); (34). The progenitor of the core-collapse Cas A supernova had the radius of the order of cm. On the other hand white dwarfs that are progenitors of Ia supernova explosions have small radii of the order of cm. Interior magnetic field G of the white dwarf will drop down to G after the homogenous expansion of the young SNR. Such a weak magnetic field can result in the ineffective DSA at the reverse shock of Ia supernovae. This effect is probably observed in Tycho SNR (see Warren et al. (35) for details).
Another possibility to suppress the CR production is related with the Alfvén drift downstream of the forward shock (9). It results in the steeper spectrum of CRs accelerated at the forward shock. On the other hand the Alfvén drift downstream of the reverse shock may produce even an opposite effect because the CR gradient is positive in this region (see Fig.2). As a result the input of the reverse shock will be significant at TeV energies.
The over-all spectrum according to this model is shown in Fig.10. The value of instead of is used downstream of the shocks. Because of the Alfvén drift the positron spectrum at reverse shock is significantly harder than the electron spectrum at the forward shock.
It should be noted that the spectra of ions accelerated at the reverse shock are harder than the proton spectra at the forward shock (see Figs 3,8) in spite of the same level of the shock modification for both shocks. This is because the shocks propagate in the media with different properties. When the reverse shock reach the flat part of the ejecta density distribution it propagates in the medium with decreasing in time density. That is why the number of freshly injected ions is low in comparison with higher energy particles accelerated earlier. It results in the spectral hardening. This effect is absent at the forward shock propagating in the medium with a constant density. The Alfvén drift strengthens this effect (see Fig.10).
We adjust the electron (positron) injection to produce a sufficient number of electrons and positrons. The spectra shown in Figs 8 and 10 can explain CR Galactic electrons and positrons. The expected positron injection efficiency from the radioactive decay of Ti is estimated as (28) for Ti mass of the order of as observed in SNRs (36). We used this number in our simulation. As for the electron injection at the forward shock the relative number of energetic electrons from photo-ionization of accelerated single charged He ions is of the order of . Here is the gamma-factor of single-charged He ion photo-ionized by galactic ultra-violet photons with energy eV and eV is the ionization potential of helium. This slightly overestimates the electron injection in young SNRs where the acceleration is fast enough, is closer to and the ionization is provided by eV optical photons. However we used this crude estimate that is justified in old SNRs for electron injection in our simulations. This injection mechanism suggested by Morlino (27) produces one order of magnitude higher number of energetic electrons in comparison with the number of Compton scattered electrons energized by gamma-photons from Co radioactive decay in supernova ejecta (28). Here is the forward shock radius expressed in parsecs.
The efficiency of electron (positron) acceleration at the reverse shock in SNRs can be observationally checked in the gamma-ray band. Young and middle-age ( yr) SNRs can be bright in TeV energies (see Fig.9). The electrons (positrons) accelerated before at the reverse shock do not strongly influenced by synchrotron losses in the weak magnetic field of the central part of the SNR and emit IC gamma-rays.
The main purpose of the present paper was the presentation of a new numerical code for the modeling of hydrodynamics and nonlinear shock acceleration in SNR. The code is described in details in five Appendixes. It develops the modeling of particle acceleration by spherical shocks fulfilled earlier by other authors (7); (8). Besides some important technical details, the main novel features of our code include the particle acceleration by two shocks - forward and reverse, the evolving with a SNR age magnetic field, which determines the cosmic ray transport, the account for the Alfvén drift effects both upstream and downstream of the forward and reverse shocks. The first version of the model of acceleration (9) and a number of important applications including the explanation of the overall spectrum of galactic cosmic rays (10), and the modeling of particle acceleration and nonthermal radiation in SNR RX J1713.7-3946 (11) were developed upon our code refinements.
The special attention in the astrophysical applications of the modeling discussed in the present work was focused on the particle acceleration by the reverse SNR shock. It was shown that the reverse shock can give a non-negligible contribution to the production of CR ions and positrons as compared to the contribution of the forward shock. The spectra of particles accelerated at the reverse shock can be harder than the spectra at the forward shock, see Figures 8 and 10. It may offer a new interpretation of cosmic ray data (37) that suggests the presence of high energy primary positrons. The acceleration of positrons by the reverse shock moving through the ejecta material, which contains the positrons from Ti radioactive decays, was proposed in (28). It is also possible (38) that the hard spectrum of accelerated nuclei and the low abundance of hydrogen in the material of supernova ejecta accelerated by the reverse shock may explain the difference in the energy spectra of hydrogen and helium in cosmic rays. We plan detailed study of all these effects in a future work.
In the light of the problem of electron injection it is of interest that the models of suprathermal electron injection (28); (27) reproduce the required amount of Galactic CR electrons and positrons if the leptons are pre-accelerated up to MeV in the upstream regions of supernova shocks.
The work was supported by the Russian Foundation for Basic Research grant 10-02-00110a. We thank the anonymous referee for a number of valuable comments.
Appendix A Details of the numerical method
We used the following change of variables upstream of the shocks:
and downstream of the shocks:
We shall use the dimensionless parameters instead of time and instead of momentum . It is convenient to use the new variable instead of momentum distribution . For the relativistic momenta this variable corresponds to a partial energy density of cosmic rays.
In these new variables the equations (1)-(4) have the following form upstream of the shocks:
Here indexes and are referred to the forward () and reverse shock () respectively.
The equations (1)-(4) have the following form downstream of the shocks:
Here the speeds and are determined as and respectively. The corresponding CR advective speeds and are determined as and respectively. The quantities are the distances between the forward shock and the contact discontinuity and between the reverse shock and the the contact discontinuity . We introduced energy density of the gas . Then the equations (A7)-(A9) are written in the conservative form that is convenient for the numerical method. The radius in the last equations should be expressed using Eq. (A2):
Appendix B Solution of hydrodynamic equations in the upstream regions
Since the flow upstream of the shocks is supersonic the hydrodynamic equations can be solved using the following implicit numerical scheme in these regions. Let introduce the grid . Integration of equations (A3)-(A5) on from to results in the following expressions for density , velocity and the gas pressure upstream of the forward shock () at the grid knot with a number in terms of and :
Here the superscript ’’ is referred to quantities at the previous time step and is the step in the dimensionless time . The quantity is given by the expression
In spite of implicity Eqs (B1-B3) are readily solved recurrently from down to corresponding to the position of the forward shock. The value of the forward shock radius is related with the forward shock radius at the previous time step as in Eqs (B1-B4).
Performing the integration of equations (A3)-(A5) on from to upstream of the reverse shock we obtain the expressions for for density, velocity , and and the gas pressure upstream of the reverse shock () in terms of and :
The quantities and are given by expressions
Eqs (B5-B7) are solved recurrently from to corresponding to the position of the reverse shock. The value of the new reverse shock radius in Eqs (B5-B8) was extrapolated using the reverse shock radius and velocity as .
Appendix C Solution of hydrodynamic equations in the downstream regions
Hydrodynamical quantities are determined at the centers of the cells with the indexes in the downstream regions. We used an explicit numerical scheme of Trac & Pen (15) for solution of Eqs. (A7)-(A9). In accordance with these equations we use the following variables
The variables are introduced as
Here is the so called freezing speed.
The values at the next instant of time are given by
Here is the grid step downstream of the forward and reverse shock respectively. Note that a uniform grid in the downstream region was used.
The fluxes at the grid knots between discontinuities are given by
where the fluxes at the grid knots are given by
Here the function is the nonlinear flux limiter.
The numerical scheme is reduced to the first order Godunov’s one for in Eqs. (C8), (C9). Such a scheme is very dissipative. The simplest second order scheme corresponds to . However it is unstable. One should use more complex flux limiters for stability. The different types of corresponding nonlinear limiters are available. Trac & Pen (15) used a so-called Van-Leer limiter:
We used a more dissipative ”minmod” limiter:
For calculations of fluxes near discontinuities we also use or .
The fluxes and just downstream of the forward and reverse shock at the knots with numbers and are calculated using hydrodynamical quantities just upstream of the forward and reverse shocks:
The fluxes at the position of the contact discontinuity at the knot with number are given by