© ESO, 2016

1. Introduction

Diffusive shock acceleration (DSA) at the forward shock of a supernova remnant (SNR) is an efficient process that relies on self-generated 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 quasi-equilibrium develops that slowly changes on account of the MHD evolution. Under the assumption of steady state for both turbulence and particle transport, the cosmic-ray 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 steady-state 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 steady-state assumption for turbulence and CRs up to the highest energies is questionable.

Here, we introduce a fully time-dependent calculation of the cosmic-ray acceleration coupled to the evolution of isotropic Alfvénic turbulence using an analytical self-similar description of the SNR magnetohydrodynamics. The advantage of our method lies in the full account of the time evolution of the hydrodynamical flow, the cosmic-ray distribution, and magnetic turbulence. Among the simplifications are the neglect of cosmic-ray feedback and the treatment of the shock as infinitely thin with parametrized cosmic-ray injection. The shock structure and particle pre-acceleration 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. Steady-state Monte Carlo studies including full feedback on the other hand may overestimate the available growth-time. 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₀ = 10^-6) by five orders of magnitude to its steady-state value, considering the intensity growth time of 16 yr. The simulation setup with a precursor size of 0.5 pc, which is reasonable for PeV-scale particles in a 40 μG field assuming Bohm-diffusion, provides only 100 yr of time before the plasma has flown to the shock. Our study is complementary to simulations of the small-scale physics operating at shocks and to steady-state 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 steady-state 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 cosmic-ray acceleration we use a kinetic approach in the test-particle 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 time-dependent wave growth. We numerically solve the time-dependent transport equation for the differential number density of cosmic rays in spherically symmetric geometry (Skilling 1975): $\begin{array}{ccc}\frac{\mathit{\partial N}}{\mathit{\partial t}}& \mathrm{=}& \end{array}$(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 thermal-leakage 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)³. 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 cosmic-ray 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 Bohm-like, i.e., $\begin{array}{ccc}{\mathit{D}}_{\mathit{r}}\mathrm{=}\frac{\mathit{v}}{\mathrm{3}}{\mathit{r}}_{\mathrm{g}}\mathit{\eta ,}& & \end{array}$(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 magnetic-field fluctuations. We assume that CRs are being scattered by Alfvén waves that satisfy the resonance condition $\begin{array}{ccc}{\mathit{k}}_{\mathrm{res}}\mathrm{=}\frac{\mathit{q}{\mathit{B}}_{\mathrm{0}}}{\mathit{pc}}\mathit{,}& & \end{array}$(3)where k_res is the wavenumber, q is the particle charge, and B₀ is the background magnetic field. Then, the diffusion coefficient reads (Bell 1978; Blandford & Eichler 1987) $\begin{array}{ccc}{\mathit{D}}_{\mathit{r}}& \mathrm{=}& \end{array}$(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₀. 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 time-dependent calculation of the spectrum of magnetic turbulence is clearly needed to derive the self-consistent diffusion coefficient in the precursor region and should result in a more realistic and self-consistent picture of cosmic-ray 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 $\begin{array}{ccc}{\mathit{B}}_{\mathrm{tot}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}\mathrm{+}\mathit{\delta b,}& & \end{array}$(5)where δb is the combined amplitude of all waves present at the given position. Averaging the energy density over sufficiently large times gives $\begin{array}{ccc}{\mathit{B}}_{\mathrm{tot}}^{\mathrm{2}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}\mathrm{⟨}\mathit{\delta }{\mathit{b}}^{\mathrm{2}}\mathrm{⟩}\mathit{.}& & \end{array}$(6)The total energy density in the waves can be represented as $\begin{array}{ccc}\mathrm{⟨}\mathit{\delta }{\mathit{b}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{\int }{\mathit{W}}_{\mathit{w}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{k}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{\int }{\mathit{E}}_{\mathit{w}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\ln \mathit{k,}& & \end{array}$(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 $\begin{array}{ccc}{\mathrm{ℰ}}_{\mathit{w}}\mathrm{=}\mathrm{\int }{\mathit{E}}_{\mathit{w}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{\delta }\mathrm{(}\mathit{k}\mathrm{-}{\mathit{k}}_{\mathrm{res}}\mathrm{)}\hspace{0.17em}\mathrm{d}\ln \mathit{k,}& & \end{array}$(8)and the diffusion coefficient of a particle with momentum p moving in the background field B₀ 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: $\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{E}}_{\mathit{w}}}{\mathit{\partial t}}\mathrm{+}{u}\mathrm{·}\mathrm{\left(}\mathrm{\nabla }{\mathit{E}}_{\mathit{w}}\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathit{w}}\mathrm{(}\mathrm{\nabla }\mathrm{·}{u}\mathrm{)}{\mathit{E}}_{\mathit{w}}\mathrm{+}\mathit{k}\frac{\mathit{\partial }}{\mathit{\partial k}}\left({\mathit{k}}^{\mathrm{2}}{\mathit{D}}_{\mathit{k}}\frac{\mathit{\partial }}{\mathit{\partial k}}\frac{{\mathit{E}}_{\mathit{w}}}{{\mathit{k}}^{\mathrm{3}}}\right)\mathrm{=}\\ & & \end{array}$(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 cosmic-ray transport Eq. (1); this leads to a system of two coupled equations, which we solve numerically using implicit finite-difference 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 wave-driving process considered in this paper. The growth rate due to resonant amplification is given as (Skilling 1975; Bell 1978) $\begin{array}{ccc}{\mathrm{\Gamma }}_{\mathrm{g}}& \mathrm{=}& \end{array}$(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 pitch-angle scattering (Winske & Leroy 1984; Niemiec et al. 2010). Presumably, the mean free path for scattering is small only for very low-energy 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 cosmic-ray 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 e-folding times. Following Niemiec et al. (2008), assuming Bohm diffusion and comparing the growth-time to the shock-capture time, one finds for the available number of e-folding times N, $\begin{array}{ccc}\mathit{N}& \mathrm{=}& \mathrm{3.3}\left(\frac{{\mathit{v}}_{\mathrm{sh}}}{\mathrm{4000}{\mathrm{km}\hspace{0.17em}\mathrm{s}}^{-1}}\right){\left(\frac{{\mathit{v}}_{\mathrm{A}}}{\mathrm{10}{\mathrm{km}\hspace{0.17em}\mathrm{s}}^{-1}}\right)}^{-1}\left(\frac{{\mathit{U}}_{\mathrm{CR}}}{\mathrm{0.1}{\mathit{U}}_{\mathrm{bulk}}}\right)\mathit{,}\end{array}$(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 time-dependent treatment another nonlinearity would occur as initially only low-energetic 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 Ion-cyclotron damping plays a role and cascading is fast. The resulting total growth rate Γ_g − Γ_d would thus be smaller.

Recently, fast-mode waves were found to be efficient scatterers of cosmic rays through both resonant and nonresonant interactions (Yan & Lazarian 2004). Compressive modes such as fast-mode 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 large-scale modes that are driven by the cosmic-ray pressure or are produced as response to small-scale turbulence (Beresnyak et al. 2009; Bykov et al. 2011; Schure & Bell 2011).

3.4. Wave damping

For wave damping we consider neutral-charged collisions and ion-cyclotron damping. The damping rate due to neutral-charged collisions is given as (Kulsrud & Cesarsky 1971; Bell 1978) $\begin{array}{ccc}{\mathrm{\Gamma }}_{\mathrm{d}\mathit{,}\mathrm{nc}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{n}}_{\mathrm{H}}\mathrm{⟨}\mathit{v\sigma }\mathrm{⟩}& & \end{array}$(12)where n_H is the number density of neutral hydrogen and ⟨vσ⟩ is the velocity-weighted cross section, averaged over the random velocity of ions. Normally, neutral-charged 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 neutral-charged damping is negligible since no molecular clouds or other regions of sufficient low temperature and ionization fraction are considered.

Ion-cyclotron damping is due to interaction of Alfvén waves with the thermal particles of the plasma and is strongest at small scales $\begin{array}{ccc}{\mathrm{\Gamma }}_{\mathrm{d}\mathit{,}\mathrm{IC}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{A}}\mathit{c}{\mathit{k}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\omega }}_{\mathrm{P}}}\mathit{,}& & \end{array}$(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 $\begin{array}{ccc}{\mathit{D}}_{\mathit{k}}\mathrm{=}{\mathit{k}}^{\mathrm{3}}{\mathit{v}}_{\mathrm{A}}\sqrt{\frac{{\mathit{E}}_{\mathit{w}}}{\mathrm{2}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}}\mathrm{·}& & \end{array}$(14)If cascading is the dominant process, this phenomenological treatment will result in a Kolmogorov-like 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, $\begin{array}{ccc}{\mathit{D}}_{\mathrm{0}}\mathrm{=}{\mathrm{10}}^{\mathrm{27}}{\left(\frac{\mathit{pc}}{\mathrm{10}\mathrm{GeV}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\left(\frac{{\mathit{B}}_{\mathrm{0}}}{\mathrm{3}\mathit{\mu }\mathrm{G}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}\mathit{.}& & \end{array}$(15)The value of D₀ 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 self-consistent amplification is orders of magnitudes higher than the initial intensity, the diffusion of particles near the shock is entirely governed by the self-generated 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 cosmic-ray 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 cosmic-ray acceleration should be most efficient, but this stage is brief for a typical type-Ia 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 type-Ia 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, Ew (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 self-consistent treatment (thick lines) and Bohm-like diffusion (thin lines) are presented. Spectra for the location at r = 1.15 Rshock are averages over a shell with thickness 0.5 pc. Particles with kinetic energy of 1 GeV are resonant with waves of wavenumber k0.

For core-collapse 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 core-collapse SNRs require a thorough numerical modeling, while type-Ia SNRs can be well described analytically. Besides, it was shown (Acharya et al. 2015) that type-Ia 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 core-collapse SNRs.

As we aim to introduce our method, we consider an SNR in the adiabatic (Sedov-Taylor) stage evolving in a typical ISM with density of 0.4 cm^-3 and a magnetic-field strength (B₀) of 5μG. We use analytic approximations for Sedov-Taylor solutions (Cox & Franco 1981) to describe the evolution of the plasma-flow parameters. The corresponding self-similar 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 power-law 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 magnetic-turbulence 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 free-expansion solutions. We provide a comparison with Sedov SNRs of the same ages. The results obtained from self-consistent 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 cosmic-ray pressure always remains below 10% of the shock ram pressure.

5.1. Sedov-Taylor 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 low-energy particles dominate the cosmic-ray 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 power-law 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 plateau-like 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 low-energy particles to substantially impact the turbulence growth, whereas high-energy 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 cascading-dominated 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 magnetic-field amplification, so that only particles at the highest energy were escaping the shock. Later on, the cosmic-ray 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 self-consistent, fully time-dependent coupled treatment of magnetic turbulence and cosmic-ray acceleration, we compare the particle spectra obtained here to those derived earlier assuming Bohm diffusion. The downstream volume-integrated spectra for that case are shown at the bottom panel of Fig. 1.

Both calculations show similar power-law 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 self-consistent 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, cosmic-ray escape is possible only at the high-energy 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 log-parabola (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 self-consistently 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. Steady-state 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 hybrid-simulations of Caprioli & Spitkovsky (2014), which include the nonlinear contributions, do not show any deviation from the DSA-predicted 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 cosmic-ray feedback.

Fig. 2 Spectral evolution of turbulence energy density, Ew (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, Ew and N derived with Bohm-like diffusion are given with thin lines. The shell and k0 are the same as in Fig. 1.

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

In general, the time-dependent 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 cosmic-ray 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 Sedov-Taylor stage are somewhat affected by our use of a test-particle approach. Keeping the cosmic-ray 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 X-ray 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 time-dependent treatment would be highly desirable.

5.2. Free expansion stage

We simulated two cases for the free-expansion stage of the SNR: (i) a shock expanding in a uniform medium of constant density; and (ii) a shock expanding into a wind-blown bubble. We terminate our calculations when the age reaches 1000 yr.

5.2.1. Turbulence spectra

The turbulence spectra for the self-consistent 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 free-expansion runs, the maximum level of turbulence is higher, and spectra extend to lower k compared to Sedov-Taylor cases of the same age. This is a consequence of the slower shock deceleration than in Sedov-Taylor 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 Sedov-Taylor phase.

There are slight differences in spectral evolution between the two free-expansion 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 uniform-medium simulations and, with time, this difference becomes more pronounced.

Besides, in wind-zone 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 upstream-turbulence 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 time-dependent treatment and for the Bohm assumption are parallel, both displaying a Kolmogorov-like k^{− 2 / 3} scaling for log (k) > −3.25. In the case with a windzone, the self-consistently 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 low-energetic particles that could amplify turbulence.

5.2.2. Particle spectra

The particle spectra shown in Figs. 2 and 3 for the uniform and wind-zone scenario, respectively, resemble those obtained with Bohm-like 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 Bohm-like diffusion, but stays fairly constant (see distributions of escaped particles in Figs. 2 and 3). The slowly changing velocity during the free-expansion 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 Sedov-Taylor 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 time-dependent transport equations for magnetic turbulence and cosmic rays. The equations are solved in spherically symmetric geometry but so far are limited to the test-particle 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 plateau-like region at intermediate k that resembles that leading to Bohm-like diffusion;

• a cascade-dominated tail at high k, where the spectrum is Kolmogorov-like;

• a damping-dominated part at very high k.

A wide plateau-like region in the k-spectrum 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 cosmic-ray 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 Bohm-like diffusion. The more efficient escape in the self-consistent 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 high-energy γ-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 cosmic-ray 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 steady-state 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

All Figures

Fig. 1 Spectral evolution of turbulence energy density, Ew (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 self-consistent treatment (thick lines) and Bohm-like diffusion (thin lines) are presented. Spectra for the location at r = 1.15 Rshock are averages over a shell with thickness 0.5 pc. Particles with kinetic energy of 1 GeV are resonant with waves of wavenumber k0.

In the text

	Fig. 2 Spectral evolution of turbulence energy density, Ew (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, Ew and N derived with Bohm-like diffusion are given with thin lines. The shell and k0 are the same as in Fig. 1.
In the text

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