Transport of magnetic turbulence in supernova remnants
^{1} DESY, 15738 Zeuthen, Germany
email: robert.brose@desy.de
^{2} Institute of Physics and Astronomy,
University of Potsdam, 14476
Potsdam,
Germany
email: igor.telezhinsky@desy.de
Received:
11
September
2015
Accepted:
4
June
2016
Context. Supernova remnants are known as sources of Galactic cosmic rays for their nonthermal emission of radio waves, Xrays, and gamma rays. However, the observed soft broken powerlaw spectra are hard to reproduce within standard acceleration theory based on the assumption of Bohm diffusion and steadystate calculations.
Aims. We point out that a timedependent treatment of the acceleration process together with a selfconsistent treatment of the scattering turbulence amplification is necessary.
Methods. We numerically solve the coupled system of transport equations for cosmic rays and isotropic Alfvénic turbulence. The equations are coupled through the growth rate of turbulence determined by the cosmicray gradient and the spatial diffusion coefficient of cosmic rays determined by the energy density of the turbulence. The system is solved on a comoving expanding grid extending upstream for dozens of shock radii, allowing for the selfconsistent study of cosmicray diffusion in the vicinity of their acceleration site. The transport equation for cosmic rays is solved in a testparticle approach.
Results. We demonstrate that the system is typically not in a steady state. In fact, even after several thousand years of evolution, no equilibrium situation is reached. The resulting timedependent particle spectra strongly differ from those derived assuming a steady state and Bohm diffusion. Our results indicate that proper accounting for the evolution of the scattering turbulence and hence the particle diffusion coefficient is crucial for the formation of the observed soft spectra. In any case, the need to continuously develop magnetic turbulence upstream of the shock introduces nonlinearity in addition to that imposed by cosmicray feedback.
Key words: ISM: supernova remnants / acceleration of particles / turbulence
© ESO, 2016
1. Introduction
Diffusive shock acceleration (DSA) at the forward shock of a supernova remnant (SNR) is an efficient process that relies on selfgenerated turbulence (Blandford & Eichler 1987). Streaming cosmic rays (CRs) in the upstream region of the shock generate magnetic turbulence that enhances the acceleration process, which in turn leads to further turbulence growth. This process is terminated by escape of CRs or generally when the growth time of turbulence becomes longer than the evolutionary timescale of the system. Conventionally, it is considered that turbulence growth is the fastest process in the system, which is then followed by particle acceleration and, finally, by global magnetohydrodynamical (MHD) evolution of the SNR. It is then often assumed that for the first two processes a quasiequilibrium develops that slowly changes on account of the MHD evolution. Under the assumption of steady state for both turbulence and particle transport, the cosmicray distribution at the shock can be derived accounting for their feedback (Blasi 2002; Caprioli et al. 2009). It can then be imposed on global SNR models (Ellison et al. 2012), or, alternatively, one can solve the entire coupled system of turbulence, CRs, and SNR fluid under steadystate conditions (Bykov et al. 2014). It has been realized, however, that limited time is available both for turbulence growth and particle acceleration for SNR (Lagage & Cesarsky 1983; Bell et al. 2013; Schure & Bell 2013), and so the steadystate assumption for turbulence and CRs up to the highest energies is questionable.
Here, we introduce a fully timedependent calculation of the cosmicray acceleration coupled to the evolution of isotropic Alfvénic turbulence using an analytical selfsimilar description of the SNR magnetohydrodynamics. The advantage of our method lies in the full account of the time evolution of the hydrodynamical flow, the cosmicray distribution, and magnetic turbulence. Among the simplifications are the neglect of cosmicray feedback and the treatment of the shock as infinitely thin with parametrized cosmicray injection. The shock structure and particle preacceleration can be well studied with kinetic simulations, but even very large hybrid simulations (e.g., Caprioli & Spitkovsky 2014) cover at most about one hour of real time and a region not larger than one astronomical unit. Steadystate Monte Carlo studies including full feedback on the other hand may overestimate the available growthtime. In Bykov et al. (2014) about 180 yr of optimal wave growth would be needed to amplify the intensity of turbulence at the longest wavelengths (k/k_{0} = 10^{6}) by five orders of magnitude to its steadystate value, considering the intensity growth time of 16 yr. The simulation setup with a precursor size of 0.5 pc, which is reasonable for PeVscale particles in a 40 μG field assuming Bohmdiffusion, provides only 100 yr of time before the plasma has flown to the shock. Our study is complementary to simulations of the smallscale physics operating at shocks and to steadystate Monte Carlo studies including full feedback (Bykov et al. 2014) that neglect the evolution of the remnant. With our method we can examine the system over very long time, ten of thousands of years, and cover a considerable fraction of the lifetime of a SNR. Our model accounts for competing plasma processes without the assumption of a steadystate situation. This renders it possible to study additional effects that arise from the limited time that is available for turbulence growth and particle acceleration.
As a result we obtain the energy density of magnetic turbulence and particle spectra at any position on the grid, which allows us to study the propagation of the CRs and their escape from the source.
2. Particle acceleration
To model cosmicray acceleration we use a kinetic approach in the testparticle approximation (Telezhinsky et al. 2012a,b, 2013). The feedback of accelerated CRs on the shock structure is negligible, as long as the CR pressure stays below 10% of the shock ram pressure (Kang & Ryu 2010). Thus the acceleration process can be treated independent of the SNR evolution as long as the amount of energy contained in CRs is limited. In any case, the purpose of the present paper is to isolate and discuss the impact of timedependent wave growth. We numerically solve the timedependent transport equation for the differential number density of cosmic rays in spherically symmetric geometry (Skilling 1975): (1)where N is the differential number density of cosmic rays, D_{r} is the spatial diffusion coefficient, u is the advective velocity, ṗ are the energy losses, and Q is the source of thermal particles that is treated according to the thermalleakage injection model (Gabici et al. 2005). The transport equation is transformed to a frame comoving with the shock in which the radial coordinate is normalized to the shock radius, x = r/R_{sh}. To resolve the precursor of lowest energy CRs, we make another transformation, (x − 1) = (x^{∗} − 1)^{3}. For equidistant binning in x^{∗} this transformation guarantees a very good resolution close to the shock while simultaneously extending out to several tens of shock radii in the upstream region, allowing for all injected particles to remain in the simulation domain.
2.1. Diffusion coefficient
One of the crucial but still poorly known parameters for the acceleration process and subsequent propagation is the spatial diffusion coefficient (Yan et al. 2012; Telezhinsky et al. 2012b). This coefficient governs the efficiency of cosmicray acceleration and thus the maximum energy reached by the cosmic rays. It is also responsible for the spatial distribution of accelerated particles both upstream and downstream of the shock, which in turn impacts the subsequent emission from the source and its vicinity. The diffusion coefficient is usually assumed to be Bohmlike, i.e., (2)where v is the particle velocity, r_{g} is its gyroradius, and η is the ratio of background magnetic energy density to energy density in magnetic fluctuations and usually assumed to be order of unity for particles of all energies. The following arguments suggest that this approach is oversimplified. As was noted, the diffusion coefficient is directly connected to the magneticfield fluctuations. We assume that CRs are being scattered by Alfvén waves that satisfy the resonance condition (3)where k_{res} is the wavenumber, q is the particle charge, and B_{0} is the background magnetic field. Then, the diffusion coefficient reads (Bell 1978; Blandford & Eichler 1987) (4)where ℰ_{w} denotes the energy density per unit logarithmic bandwidth of Alfvén waves resonant with particles of momentum p according to resonance condition (3) and U_{m} is the energy density of the background magnetic field B_{0}. Bohm diffusion (2), for which η is a constant, is then equivalent to a featureless and flat magnetic turbulence spectrum. As pristine plasma is continuously advected toward the shock, however, magnetic turbulence must be constantly replenished in the upstream region. If the turbulence is a result of the growth of Alfvén waves (for instance owning to resonant amplification by streaming CRs) as well as their spatial transport, compression at the shock, damping by various mechanisms, and spectral energy transfer due to turbulence cascading, it is obviously very unlikely that the turbulence spectrum in the SNR and its vicinity is flat and featureless. Besides, from γray observations of SNRs and their surroundings we now understand that (i) it is hard to accommodate Bohm diffusion for particles in an energy band as wide as we see in SNRs; and (ii) the diffusion around SNRs is much slower than the average Galactic diffusion, but much faster than Bohm. This is in fact the expected behavior as turbulence must continuously be generated in the upstream region of the shock (Fujita et al. 2011; Yan et al. 2012). A timedependent calculation of the spectrum of magnetic turbulence is clearly needed to derive the selfconsistent diffusion coefficient in the precursor region and should result in a more realistic and selfconsistent picture of cosmicray acceleration in SNRs.
3. Magnetic turbulence
3.1. Spectral energy density
We consider Alfvén waves as scattering centers for CRs. Alfvén waves can be considered a small contribution to the magnetic field at some position, so that (5)where δb is the combined amplitude of all waves present at the given position. Averaging the energy density over sufficiently large times gives (6)The total energy density in the waves can be represented as (7)where W_{w} and E_{w} are the spectra of magnetic turbulence energy density per unit interval in k and lnk, respectively.
The energy density per unit logarithmic bandwidth of resonant Alfvén waves is then (8)and the diffusion coefficient of a particle with momentum p moving in the background field B_{0} can be calculated using expression (4).
3.2. Transport equation
We consider a 1D spherically symmetric geometry and treat the turbulence as isotropic. To investigate the growth, damping, and cascading of the Alfvén waves, as well as their propagation along the background magnetic field, we solve an equation for the transport of the magnetic turbulence along with the transport Eq. (1) for cosmic rays. The transport of magnetic turbulence can be described by a continuity equation for the spectral energy density, E_{w} = E_{w}(r,k,t), as follows: (9)where C_{w} = 1.5 denotes the prefactor for wave compression at the shock (McKenzie & Voelk 1982), D_{k} is the diffusion coefficient in wavenumber space representing cascading, and Γ_{g} and Γ_{d} are the growth and damping rates, respectively. Therefore, we consider growth, damping, advection, compression of turbulence at the shock, as well as spectral energy transfer through cascading. The transport equation for magnetic turbulence (9) is transformed to a frame comoving with the shock in the same manner as the cosmicray transport Eq. (1); this leads to a system of two coupled equations, which we solve numerically using implicit finitedifference methods (Guyer et al. 2009).
3.3. Wave growth
Particles streaming faster than the Alfvén speed should generate Alfvén waves at wavelengths similar to the gyroradii of the particles (Wentzel 1974, and references therein). In the diffusion limit the growth rate of waves can be related to the pressure gradient of the CRs. This mechanism is known as resonant amplification of Alfvén waves and the only wavedriving process considered in this paper. The growth rate due to resonant amplification is given as (Skilling 1975; Bell 1978) (10)Magnetic turbulence can also be produced on very small scales through Bell’s nonresonant instability (Lucek & Bell 2000; Bell 2004). The interaction of cosmic rays with this mode is also nonresonant and does not involve pure pitchangle scattering (Winske & Leroy 1984; Niemiec et al. 2010). Presumably, the mean free path for scattering is small only for very lowenergy particles whose Larmor radius is commensurate the wavelength of Bell’s mode (Bykov et al. 2014). Those particles are typically not present far out in the cosmicray precursor and so the impact of Bell’s mode is likely to be moderate at later times. Generally Bell’s mode might be operating only for a few efolding times. Following Niemiec et al. (2008), assuming Bohm diffusion and comparing the growthtime to the shockcapture time, one finds for the available number of efolding times N, (11)where U_{CR} are the energy density of the cosmic rays and U_{bulk} the energy density of the bulk plasma respectively. In our fully timedependent treatment another nonlinearity would occur as initially only lowenergetic particles are present in the simulation, whose current is quickly reduced by a decreasing diffusion coefficient on account of wave growth. The fastest growth rate would initially occur on scales where Ioncyclotron damping plays a role and cascading is fast. The resulting total growth rate Γ_{g} − Γ_{d} would thus be smaller.
Recently, fastmode waves were found to be efficient scatterers of cosmic rays through both resonant and nonresonant interactions (Yan & Lazarian 2004). Compressive modes such as fastmode waves are thermally damped with a rate that depends on the orientation of the wave vector, and so a 3D treatment of the wave spectrum would be needed which we defer to a future publication. Likewise, we ignore largescale modes that are driven by the cosmicray pressure or are produced as response to smallscale turbulence (Beresnyak et al. 2009; Bykov et al. 2011; Schure & Bell 2011).
3.4. Wave damping
For wave damping we consider neutralcharged collisions and ioncyclotron damping. The damping rate due to neutralcharged collisions is given as (Kulsrud & Cesarsky 1971; Bell 1978) (12)where n_{H} is the number density of neutral hydrogen and ⟨vσ⟩ is the velocityweighted cross section, averaged over the random velocity of ions. Normally, neutralcharged damping is relatively weak and independent of the wavenumber of Alfvén waves. This mechanism is mostly important in regions of low temperatures and high densities such as molecular clouds. As we noted before, if cosmic rays penetrate molecular clouds, their spectrum may be strongly modified because of the evanescence of magnetic turbulence that is responsible for particle scattering (Malkov et al. 2011). In this work the effect of the neutralcharged damping is negligible since no molecular clouds or other regions of sufficient low temperature and ionization fraction are considered.
Ioncyclotron damping is due to interaction of Alfvén waves with the thermal particles of the plasma and is strongest at small scales (13)where ω_{P} is the ion plasma frequency (Threlfall et al. 2011). This damping should transfer energy to the plasma via heating, which is not yet considered in this work, although we are aware that it might modify the spectrum around the scale of particle injection.
3.5. Wave cascading
The process of energy transfer through turbulence cascading from larger scales to smaller scales is not yet fully understood and is the subject of active research. Empirically it can be described as a diffusion process in wavenumber space (Zhou & Matthaeus 1990; Schlickeiser 2002). Given the assumption of isotropic turbulence, the corresponding diffusion coefficient is (14)If cascading is the dominant process, this phenomenological treatment will result in a Kolmogorovlike turbulence spectrum, E_{w} ∝ k^{− 2 / 3}. Because cascading is treated as a diffusion process, a small fraction of energy is also transferred to scales that are larger than the turbulence injection scale. This permits scattering of particles whose energy is higher than those currently driving the turbulence and is in effect similar to resonance broadening. The acceleration time for these particles is therefore reduced and the acceleration process more efficient.
3.6. Initial conditions
In the absence of initial turbulence there is nothing to grow, and so Eq. (9) requires some seed turbulence. We therefore take as an initial condition an interstellar medium (ISM) turbulence derived with Eq. (4) from the diffusion coefficient, (15)The value of D_{0} is a factor 100 lower than the ISM diffusion coefficient found in studies of Galactic CR propagation (e.g. Trotta et al. 2011), corresponding to a somewhat enhanced intensity of turbulence at large distance from the forward shock. This choice is computationally expedient and increases numerical stability during the initial stages of our simulations. As long as the energy density of the turbulence attained through selfconsistent amplification is orders of magnitudes higher than the initial intensity, the diffusion of particles near the shock is entirely governed by the selfgenerated turbulence. The intensity of turbulence is either determined by balance of growth and damping or cascading, in which case the intensity is insensitive to its initial value, or by the available time. In the latter situation our elevated initial intensity makes it easier to reach a very high level of turbulence. Any failure to generate the turbulence intensity needed for very efficient cosmicray acceleration is therefore intrinsic.
Moreover, although not particularly important for the scope of the current paper, this diffusion coefficient allows all injected particles to stay within the numerical grid.
4. Evolution of supernova remnant
The evolution of an SNR can be subdivided into three major stages: free expansion, adiabatic phase, and radiative stage. The initial free expansion is characterized by the fastest shocks, and cosmicray acceleration should be most efficient, but this stage is brief for a typical typeIa SNR. After several hundreds years the shock already sweeps up enough ISM to slow down considerably. The following adiabatic phase, given its long duration, should be most frequent among Galactic typeIa SNRs. Because the shock is still rather fast, the bulk of the CRs produced by the remnant should be accelerated during this stage. In the subsequent radiative stage, the shock becomes very slow, allowing for efficient recombination in the downstream region. The layer of cooling material just behind the radiative shock contains a large portion of neutrals and, therefore, acceleration should level off.
Fig. 1 Spectral evolution of turbulence energy density, E_{w} (top), and differential proton number density, N (bottom), for an SNR in the adiabatic stage at the age of 400 (red), 1000 (green), 3000 (blue) and 12 000 (gray) yr. Both the selfconsistent treatment (thick lines) and Bohmlike diffusion (thin lines) are presented. Spectra for the location at r = 1.15 R_{shock} are averages over a shell with thickness 0.5 pc. Particles with kinetic energy of 1 GeV are resonant with waves of wavenumber k_{0}. 

Open with DEXTER 
For corecollapse SNRs evolving in the winds of progenitor stars the picture is completely different. On account of strong variation in the wind parameters, the sequence of stages may be mixed up. Most of the corecollapse SNRs require a thorough numerical modeling, while typeIa SNRs can be well described analytically. Besides, it was shown (Acharya et al. 2015) that typeIa SNRs should be easiest to observe with the future CTA facility owing to a continuous increase of γray flux and size. Detectability and resolvability is a strong function of the age for many corecollapse SNRs.
As we aim to introduce our method, we consider an SNR in the adiabatic (SedovTaylor) stage evolving in a typical ISM with density of 0.4 cm^{3} and a magneticfield strength (B_{0}) of 5μG. We use analytic approximations for SedovTaylor solutions (Cox & Franco 1981) to describe the evolution of the plasmaflow parameters. The corresponding selfsimilar solutions for the background magnetic field profiles inside the SNR are taken from Korobeinikov (1964).
For the sake of comparison, we also performed simulations for two cases of freely expanding SNR in a uniform and powerlaw density medium. The simulations were carried out just for the first 1000 yr. For these simulations we used analytic expressions for the shock evolution (Truelove & McKee 1999), and the background magnetic field was calculated by solving the induction equation as described in Telezhinsky et al. (2013).
5. Results
We explore the evolution of magneticturbulence spectra and the corresponding particle spectra in adiabatic SNRs up to an age of 12 000 yr. We also present and discuss the results for freeexpansion solutions. We provide a comparison with Sedov SNRs of the same ages. The results obtained from selfconsistent calculations are then compared to the standard approach assuming Bohm diffusion. The specific model features Bohm diffusion everywhere downstream of the shock, whereas in the upstream region an exponential transition to the Galactic diffusion is assumed to occur between the forward shock and a location one SNR radius ahead (Telezhinsky et al. 2012b, 2013). In all calculations we use the same injection parameters. In order not to violate test particle approximation, we set the injection of thermal particles so that the cosmicray pressure always remains below 10% of the shock ram pressure.
5.1. SedovTaylor stage
5.1.1. Turbulence spectra
The evolution of magnetic turbulence is given in Fig. 1. The turbulence spectrum corresponding to Bohm diffusion in these coordinates, as can be seen from Eq. (4), would be a constant line with a value of 4/π. In contrast, our calculations show that the turbulence spectra at the shock exhibit a complicated shape; there is a very extended region of efficient growth that spans over several orders of magnitude in k space. This is not surprising because at the position of the shock the turbulence is driven by particles of all energies, and, since lowenergy particles dominate the cosmicray spectrum, the growth of turbulence is fastest at large k. However, the larger k is, the more important cascading becomes, and so at some point it starts playing a crucial role and dominates over growth. This is a region where a transition to the inertial range happens, and a break in the spectrum appears. Finally, at high k, a classical Kolmogorov powerlaw turbulence spectrum is observed. An interesting result of such an interplay between growth and cascading is that in a rather wide range of k the turbulence spectrum appears to be plateaulike and similar to that for Bohm diffusion.
This interpretation is well supported by the shape of the spectrum at some distance upstream. The turbulence spectrum at the shock of a young SNR suggests that particles in a wide energy band should be well confined to the shock. However, at low k there is a sharp cutoff in the turbulence spectrum and, thus, particles beyond some energy can freely escape from the shock region. Therefore, at some distance from the shock in the upstream region one expects to see a narrow particle distribution that is peaked at the energy of escaping particles. This particle distribution should drive turbulence only in a very narrow k region. For instance, at a location 15% of the SNR radius ahead of the shock, there are too few lowenergy particles to substantially impact the turbulence growth, whereas highenergy cosmic rays can easily diffuse further away from the shock. Hence the turbulence spectrum looks like a classical spectrum (Fig. 1, right). We observe just two regions: the first is the injection region peaked at the wavenumber of maximum growth corresponding to the peak in the particle energy spectrum and the second is a cascadingdominated region at higher wavenumbers. The damping range, which would be the third region of the classical turbulence spectrum, is not shown in Fig. 1 in the upstream region because it is not important for this work.
With passage of time, the sharp cutoffs in the turbulence spectra at the shock become mild and the peaks in the upstream spectra broader. This can be explained by the initially high growth rates and rather strong magneticfield amplification, so that only particles at the highest energy were escaping the shock. Later on, the cosmicray gradients are not as sharp, and the growth rates become low. Consequently, the turbulence level drops, allowing the escape of particles over a wide range in energy. Additionally, the advection of turbulence from far upstream broadens the observed spectral distribution at later times.
To be noted from both plots is that there is no steady state. The spectral shape and the maximum level of turbulence in both the shock and the upstream regions are continuously changing.
5.1.2. Particle spectra
To understand the impact of the selfconsistent, fully timedependent coupled treatment of magnetic turbulence and cosmicray acceleration, we compare the particle spectra obtained here to those derived earlier assuming Bohm diffusion. The downstream volumeintegrated spectra for that case are shown at the bottom panel of Fig. 1.
Both calculations show similar powerlaw indices, but the cutoff regions differ sharply. In the case of Bohm diffusion, a classical exponential cutoff is seen at very high energy. In contrast, the selfconsistent calculation systematically yields lower maximum energies and a smoother cutoff that is softer than exponential. Moreover, the spectral softening in the cutoff region shows an evolutionary trend, that arises because the diminishing turbulence amplitude provides for more effective particle escape as time passes. Escape becomes possible for particles in a broad range of energies, and the maximum particle energy decreases faster than for Bohm diffusion. In the Bohmian case of a constant turbulence amplitude, cosmicray escape is possible only at the highenergy tail of the distribution.
The evolution of the maximum energy of particles and the history of their escape from the shock is best illustrated in the distribution of the particles upstream of the shock. We consider a particle population in a shell upstream at some distance from the shock. Ahead of the shock, it is expected that the particle distribution is dominated by escaping particles. For Bohm diffusion this is a logparabola (Ellison & Bykov 2011) centered around the maximum energy, which is exactly what we observe for the Bohm case in Fig. 1 (bottom right). The contrast with the selfconsistently calculated spectra is huge, which resemble a log parabola only at the very beginning when the turbulence amplitude is sufficiently high. (Fig. 1, right). Later on, when the turbulence is weaker, the distribution of escaped particle shifts to lower energy and becomes substantially broader. One consequence is a softening and broadening of the cutoff region of the downstream particle spectrum, which under certain conditions can be fitted with a power law and a cutoff. Steadystate simulations that incorporate the CR feedback on the shock do not show this softening. Instead a hardening at the highest energies is observed (Bykov et al. 2014). The nonlinear impact of particles at the highest energy is the slowest process in the system, however, and it is possible that the limited growth time for waves scattering particles of the highest energies might attenuate the development of convex spectra and spectral hardening at the highest energies typical for NDSA. The hybridsimulations of Caprioli & Spitkovsky (2014), which include the nonlinear contributions, do not show any deviation from the DSApredicted s = −2 spectra, despite a transfer of 10−20% of the total energy to accelerated particles. In the end, the need to continuously develop magnetic turbulence upstream of the shock introduces nonlinearity in addition to that imposed by cosmicray feedback.
Fig. 2 Spectral evolution of turbulence energy density, E_{w} (top), and differential proton number density, N (bottom), for a SNR freely expanding into a uniform ISM (thick lines) at the age of 100 (black), 400 (red), and 1000 (green) years. For comparison, E_{w} and N derived with Bohmlike diffusion are given with thin lines. The shell and k_{0} are the same as in Fig. 1. 

Open with DEXTER 
Fig. 3 Same as Fig. 2 but for an SNR expanding into a wind zone. 

Open with DEXTER 
In general, the timedependent treatment introduces a connection between the maximum energy and the injection parameter. To reach higher energies more particles have to be injected because the growth rate (10) used here is proportional to the gradient of the CR distribution. This gradient determines the maximum level of turbulence that defines the cutoff energy. Moreover, the injection parameter is no longer an arbitrary parameter that simply scales the CR number density and thus photon fluxes. As in calculations of nonlinear diffusive shock acceleration with cosmicray feedback (e.g., Amato & Blasi 2006), it can be possibly constrained by accommodating both flux intensity and cutoff energy in the observational data.
The results for the SedovTaylor stage are somewhat affected by our use of a testparticle approach. Keeping the cosmicray pressure below 10% of the ram pressure at late times requires a low injection efficiency, ψ, resulting in a maximum particle energy of only about 1 TeV. Although we did not calculate electron spectra and their emission here, the cutoff of the turbulence spectrum in our simulation does not allow for the high electron energies needed to explain observations of nonthermal Xray emission from remnants such as SN1006 (Koyama et al. 1995). Our model would still be applicable to young remnants, where higher ψ values can be used. The model also indicates that there is a need for further improvements in the injection model. Arguably, the combination of particle feedback and a fully timedependent treatment would be highly desirable.
5.2. Free expansion stage
We simulated two cases for the freeexpansion stage of the SNR: (i) a shock expanding in a uniform medium of constant density; and (ii) a shock expanding into a windblown bubble. We terminate our calculations when the age reaches 1000 yr.
5.2.1. Turbulence spectra
The turbulence spectra for the selfconsistent simulations of the freely expanding SNR are shown in Figs. 2 and 3 for a uniform ambient medium and a wind zone with density ρ ∝ r^{2}, respectively. For both freeexpansion runs, the maximum level of turbulence is higher, and spectra extend to lower k compared to SedovTaylor cases of the same age. This is a consequence of the slower shock deceleration than in SedovTaylor solutions. Therefore, there are more particles of higher energy that can drive turbulence at small k. Other than that, the spectra are similar to those for the SedovTaylor phase.
There are slight differences in spectral evolution between the two freeexpansion simulations. Again they mainly arise from the difference in the evolution of the shock velocity. In the wind zone, the shock speed does not significantly decrease during the simulation. Hence, spectra at the shock and in the upstream region tend to extend to lower k than they do in uniformmedium simulations and, with time, this difference becomes more pronounced.
Besides, in windzone simulations the turbulence spectra in the upstream region initially have a somewhat concave structure beyond the peak. This can be seen comparing the red lines (400 yr) in the plot for the upstreamturbulence spectra (upper right panels) for a uniform ambient medium (Fig. 2) and for expansion into a wind zone (Fig. 3). In the case without a wind zone, the turbulence spectra for our timedependent treatment and for the Bohm assumption are parallel, both displaying a Kolmogorovlike k^{− 2 / 3} scaling for log (k) > −3.25. In the case with a windzone, the selfconsistently calculated turbulence spectrum is clearly softer than that corresponding to Bohm diffusion, at least for log (k) < −1.5, on account of the limited time available for turbulence cascading in the presence of a very young, fast shock and the absence of lowenergetic particles that could amplify turbulence.
5.2.2. Particle spectra
The particle spectra shown in Figs. 2 and 3 for the uniform and windzone scenario, respectively, resemble those obtained with Bohmlike diffusion. The maximum energies are consistently lower, though, at least for our choice of a low injection efficiency. Also, the maximum energy does not increase as it would for Bohmlike diffusion, but stays fairly constant (see distributions of escaped particles in Figs. 2 and 3). The slowly changing velocity during the freeexpansion stage is compensated by the evolution of the turbulence spectra, and we do not observe substantial growth in maximum energy. There is no spectral softening observed as was seen in SedovTaylor solutions. The simulation time of 1000 yr is too short for a decreasing turbulence level to have a significant effect on particles escape from the remnant, and the distributions of escaped cosmic rays look rather narrow and do not strongly deviate from a log parabola.
6. Conclusions
We developed a model for particle acceleration in SNR by simultaneously solving timedependent transport equations for magnetic turbulence and cosmic rays. The equations are solved in spherically symmetric geometry but so far are limited to the testparticle regime. We consider the cosmic rays being scattered by isotropic Alfvénic turbulence that is subject to compression, advection, cascading, damping and growth due to resonant amplification of Alfvén waves. The calculated turbulence spectra, independent of the hydrodynamical model, reveal the same features at the shock:

a sharp cutoff at low k numbers;

a wide plateaulike region at intermediate k that resembles that leading to Bohmlike diffusion;

a cascadedominated tail at high k, where the spectrum is Kolmogorovlike;

a dampingdominated part at very high k.
A wide plateaulike region in the kspectrum is essential for particle acceleration.
We found that even for old remnants a steady state is not reached, either in turbulence or in particle spectra. Even after 12 000 yr of SNR evolution both spectra are continuously changing, which raises concerns as to the validity of cosmicray acceleration models that rely on the assumption of a steady state, either openly or implicitly. In particular, the spectral shape and intensity of turbulence are strongly time dependent. As feedback, the particle spectra acquire significant features, which do not appear for Bohmlike diffusion. The more efficient escape in the selfconsistent treatment gives rise to the formation of softer spectra at late stages of SNR evolution. We note with interest that soft spectra, without any sign of hardening at the highest energies as predicted by NDSA models, are indeed observed in the highenergy γray band (Ackermann et al. 2013). The need to continuously develop magnetic turbulence upstream of the shock introduces nonlinearity in addition to that imposed by cosmicray feedback. Although our choice of a low injection efficiency is partly responsible for the low maximum energy of cosmic rays compared to that in the Bohm case, the evolutionary trend in maximum energy suggests that its evolutionary decrease proceeds much faster, and therefore at the age of γray emitting Galactic SNR the maximum energy is indeed lower than is estimated with steadystate models.
Acknowledgments
We acknowledge support by the Helmholtz Alliance for Astroparticle Physics HAP funded by the Initiative and Networking Fund of the Helmholtz Association.
References
 Acharya, B. S., Aramo, C., Babic, A., et al. 2015, Astropart. Phys., 62, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Amato, E., & Blasi, P. 2006, MNRAS, 371, 1251 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 1978, MNRAS, 182, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R. 2004, MNRAS, 353, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Beresnyak, A., Jones, T. W., & Lazarian, A. 2009, ApJ, 707, 1541 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R., & Eichler, D. 1987, Phys. Rep., 154, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Blasi, P. 2002, Astropart. Phys., 16, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Bykov, A. M., Osipov, S. M., & Ellison, D. C. 2011, MNRAS, 410, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Bykov, A. M., Ellison, D. C., Osipov, S. M., & Vladimirov, A. E. 2014, ApJ, 789, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Caprioli, D., & Spitkovsky, A. 2014, ApJ, 794, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Caprioli, D., Blasi, P., Amato, E., & Vietri, M. 2009, MNRAS, 395, 895 [NASA ADS] [CrossRef] [Google Scholar]
 Cox, D. P., & Franco, J. 1981, ApJ, 251, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Ellison, D. C., & Bykov, A. M. 2011, ApJ, 731, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Ellison, D. C., Slane, P., Patnaude, D. J., & Bykov, A. M. 2012, ApJ, 744, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Fujita, Y., Takahara, F., Ohira, Y., & Iwasaki, K. 2011, MNRAS, 415, 3434 [NASA ADS] [CrossRef] [Google Scholar]
 Gabici, S., Blasi, P., & Vannoni, G. 2005, in Astrophysical Sources of High Energy Particles and Radiation, eds. T. Bulik, B. Rudak, & G. Madejski, AIP Conf. Ser., 801, 369 [Google Scholar]
 Guyer, J. E., Wheeler, D., & Warren, J. A. 2009, Comput. Sci. Eng., 11, 6 [CrossRef] [Google Scholar]
 Kang, H., & Ryu, D. 2010, ApJ, 721, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Korobeinikov, V. P. 1964, Zhurnal Prikladnoi Mekhaniki i Tekhnicheskoi Fiziki, 4, 113 [Google Scholar]
 Koyama, K., Petre, R., Gotthelf, E. V., et al. 1995, Nature, 378, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Kulsrud, R. M., & Cesarsky, C. J. 1971, ApJ, 8, L189 [Google Scholar]
 Lagage, P. O., & Cesarsky, C. J. 1983, A&A, 125, 249 [NASA ADS] [Google Scholar]
 Lucek, S. G., & Bell, A. R. 2000, MNRAS, 314, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Malkov, M. A., Diamond, P. H., & Sagdeev, R. Z. 2011, Nature Commun., 2, 194 [NASA ADS] [CrossRef] [Google Scholar]
 McKenzie, J. F., & Voelk, H. J. 1982, A&A, 116, 191 [NASA ADS] [Google Scholar]
 Niemiec, J., Pohl, M., Stroman, T., & Nishikawa, K.I. 2008, ApJ, 684, 1174 [NASA ADS] [CrossRef] [Google Scholar]
 Niemiec, J., Pohl, M., Bret, A., & Stroman, T. 2010, ApJ, 709, 1148 [NASA ADS] [CrossRef] [Google Scholar]
 Schlickeiser, R. 2002, Cosmic Ray Astrophysics (Berlin: SpringerVerlag) [Google Scholar]
 Schure, K. M., & Bell, A. R. 2011, MNRAS, 418, 782 [NASA ADS] [CrossRef] [Google Scholar]
 Schure, K. M., & Bell, A. R. 2013, MNRAS, 435, 1174 [NASA ADS] [CrossRef] [Google Scholar]
 Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2012a, Astropart. Phys., 35, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2012b, A&A, 541, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2013, A&A, 552, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Threlfall, J., McClements, K. G., & De Moortel, I. 2011, A&A, 525, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Trotta, R., Jóhannesson, G., Moskalenko, I. V., et al. 2011, ApJ, 729, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Truelove, J. K., & McKee, C. F. 1999, ApJS, 120, 299 [NASA ADS] [CrossRef] [Google Scholar]
 Wentzel, D. G. 1974, ARA&A, 12, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Winske, D., & Leroy, M. M. 1984, J. Geophys. Res., 89, 2673 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., & Lazarian, A. 2004, ApJ, 614, 757 [NASA ADS] [CrossRef] [Google Scholar]
 Yan, H., Lazarian, A., & Schlickeiser, R. 2012, ApJ, 745, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, Y., & Matthaeus, W. H. 1990, J. Geophys. Res., 95, 14881 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Spectral evolution of turbulence energy density, E_{w} (top), and differential proton number density, N (bottom), for an SNR in the adiabatic stage at the age of 400 (red), 1000 (green), 3000 (blue) and 12 000 (gray) yr. Both the selfconsistent treatment (thick lines) and Bohmlike diffusion (thin lines) are presented. Spectra for the location at r = 1.15 R_{shock} are averages over a shell with thickness 0.5 pc. Particles with kinetic energy of 1 GeV are resonant with waves of wavenumber k_{0}. 

Open with DEXTER  
In the text 
Fig. 2 Spectral evolution of turbulence energy density, E_{w} (top), and differential proton number density, N (bottom), for a SNR freely expanding into a uniform ISM (thick lines) at the age of 100 (black), 400 (red), and 1000 (green) years. For comparison, E_{w} and N derived with Bohmlike diffusion are given with thin lines. The shell and k_{0} are the same as in Fig. 1. 

Open with DEXTER  
In the text 
Fig. 3 Same as Fig. 2 but for an SNR expanding into a wind zone. 

Open with DEXTER  
In the text 