Free Access
Issue
A&A
Volume 634, February 2020
Article Number A59
Number of page(s) 11
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201936567
Published online 07 February 2020

© ESO 2020

1. Introduction

Supernova remnants (SNRs) are known to accelerate cosmic-rays (CRs) to relativistic energies (Ackermann et al. 2013), and the highest energies are likely reached during the earliest phases of SNR evolution (Bell et al. 2013) and before the transition to the Sedov-phase (Ptuskin & Zirakashvili 2003). As soon as the peak maximum energy of a SNR is reached, the highest energetic CRs start to leak from the remnant. Whatever the CR spectrum inside the SNR at some point in time, the CR contribution of an SNR to the sea of Galactic CRs is given by the time integral of the CR leakage into the interstellar medium (ISM). Analytic calculations showed that the release spectra can be significantly softer than the spectra inside the SNRs (Ptuskin & Zirakashvili 2005; Ohira et al. 2010). Recently, Celli et al. (2019) developed an analytic description for both the released CR spectrum and the spectrum of CRs remaining inside the remnant. They showed that the spectrum inside the remnants can reassemble broken power-laws similar to the spectra observed in middle-aged remnants (Acciari et al. 2009; Ackermann et al. 2013). At the same time, the release spectrum can be a s = 2 power-law, if the acceleration spectrum has an index of s <  2, or can have the same spectral index as the accelerated spectrum if s >  2.

Models for the propagation of Galactic CRs indeed require injection spectra with a break at a few GV/c in rigidity that are hard at low energies and assume spectral indices around s ≈ 2.4 above the break energy (Putze et al. 2009; Trotta et al. 2011), which are softer than those produced by linear diffusive shock acceleration (DSA). The electron injection spectra appear to be even softer than those above 30 GeV (Moskalenko & Strong 1998), which most likely reflects electron energy losses prior to their release from the remnant (Diesing & Caprioli 2019). Moreover, the gamma-ray emission of the middle-aged SNRs W44 and IC443 indicates CR spectra with an spectral index of s ≈ 2.7 at the highest energies (Acciari et al. 2009; Ackermann et al. 2013) and thus even softer than the injection spectra required by the propagation models. Furthermore, most remnants seem to show spectral breaks in their high-energy gamma-ray spectra with a break energy that typically decreases with age (Zeng et al. 2019).

Recent models for CR acceleration that include the CR feedback on the plasma flow, as well as the amplification of magnetic turbulence, are either steady-state calculations (e.g Bykov et al. 2014) or based on kinetic simulations that only cover very short time-scales (Caprioli et al. 2009). The late phases of CR acceleration have been studied by Zirakashvili & Ptuskin (2012) with a nonlinear DSA model not considering the self-consistent amplification of turbulence. Reacceleration of CRs at slow shocks may arise (Cardillo et al. 2016), as does a spectral modification, due to an enhanced escape-rate of CRs (Malkov et al. 2011). Furthermore, energy transfer from CRs to magnetic turbulence (Bell et al. 2019) and deviations from spherical symmetry (Malkov & Aharonian 2019) might lead to softer CR-spectra.

We present a time-dependent test particle calculation of CR acceleration over most of the lifetime of an SNR, including the beginning of the post-adiabatic phase. We show that the maximum energy of the accelerated CRs evolves only by one order of magnitude for Bohm-like diffusion, but much more strongly than that if the self-consistent amplification of Alfvén waves is taken into account. The spectra inside the remnant, as well as the total production-spectra, feature spectral breaks and softer spectral indeces than those predicted by standard DSA.

2. Basic equations and assumptions

In this section, we present the basic assumptions we made in our numeric approach to the DSA problem. We combined a kinetic treatment of the CRs with a thermal leakage injection model, a fully time-dependent treatment of the magnetic turbulence, and a PLUTO-based hydro calculation.

2.1. Cosmic rays

We model the acceleration of cosmic rays using a kinetic approach in the test-particle approximation (Telezhinsky et al. 2012a,b, 2013), and we choose parameters for which the CR pressure always stays below 10% of the shock ram pressure. The time-dependent transport equation for the differential number density of cosmic rays N (Skilling 1975) is given by:

N t = ( D r N u N ) p ( ( N p ˙ ) · u 3 N p ) + Q , $$ \begin{aligned} \frac{\partial N}{\partial t} = \nabla (D_r\nabla N-\boldsymbol{u} N) -\frac{\partial }{\partial p}\left( (N\dot{p})-\frac{\nabla \cdot \boldsymbol{u}}{3}Np\right)+Q , \end{aligned} $$(1)

where Dr denotes the spatial diffusion coefficient, u the advective velocity, energy losses and Q the source of thermal particles.

We solve this transport equation in a frame co-moving with the shock. The radial coordinate is transformed according to (x − 1) = (x* − 1)3, where x = r/Rsh. For a equidistant binning of x* this transformation guarantees a very fine resolution close to the shock, and an outer grid boundary that extends to several tens of shock-radii upstream for x* ≫ 1. Thus, all accelerated particles can be kept in the simulation domain.

The background of Galactic cosmic rays is introduced as initial condition outside of the remnant and as boundary condition for the differential cosmic-ray density very far upstream, which is equivalent to assuming an infinite supply of cosmic rays to the boundary. We describe the spectrum of hadronic CRs as a power law in total energy, modified at low energy by the particle speed, β. The electron spectrum is a log-parabola at low energies,

N PR ( E ) = N 1 β ( E ) ( E + m c 2 ) s 1 $$ \begin{aligned}&N_\mathrm{PR} (E) = N_1\beta (E) (E+mc^2)^{s_1}\end{aligned} $$(2)

N EL ( E ) = { N 2 E exp ( log 2 ( E / E c ) σ ) for E E B N 3 E s 2 for E > E B . $$ \begin{aligned}&N_\mathrm{EL} (E) = {\left\{ \begin{array}{ll} \frac{N_2}{E}\exp \left(-\frac{\log ^2({E}/{E_c})}{\sigma }\right)\,\mathrm{for}\,E\le E_B\\ N_3E^{s_2} \,\mathrm{for}\,E>E_B. \end{array}\right.} \end{aligned} $$(3)

Both electron and proton background spectra can be directly measured above a few GeV, where solar modulation is unimportant. Whereas the galactic electron spectrum can be constrained by measurements of diffuse radio emission, the spectral slope of the proton spectrum at low energies remains unclear. At high energies the local CR spectrum has an index s1 = −2.75, which gives a spectrum harder than E−2 at low energies; we chose the normalization in accordance with Moskalenko et al. (2002). We also investigated an alternative, softer background spectrum and discuss the differences in Appendix B.

For the electron spectrum, we fit the spectra given in Jaffe et al. (2011) for a galactocentric radius of 6.5 kpc with expression (3) as the spectrum shows a continuous change in the index below 4 GeV. The spectral index at high energies is found to be s2 = −3.04 and EB = 5 GeV. These spectra are compatible with direct observations of the electron spectra in the local ISM by Voyager 1, which also show spectra harder than s = 2.0 at low energies (Cummings et al. 2016).

We inject particles according to the thermal leakage injection model (Blasi et al. 2005; Malkov 1998). Here, the efficiency of injection ηi is given by

η i = 4 3 π ( σ 1 ) ψ 3 e ψ 2 , $$ \begin{aligned} \eta _i = \frac{4}{3\sqrt{\pi }}(\sigma -1)\psi ^3\mathrm{e}^{-\psi ^2}, \end{aligned} $$(4)

where σ is the shock compression ratio, and ψ is the multiple of the thermal momentum, at which we inject particles. Hanusch et al. (2019) suggests an additional dependence of ηi on the shock Mach number. However, this behaviour is only observed for low-Mach-number shocks and will be ignored here. Reville & Bell (2013) used a spherical harmonics expansion of the cosmic-ray Fokker–Planck equation to find a quasi-universal behaviour of shocks irrespective of the magnetic field orientation very far upstream of the shocks, which suggests that injection is only weakly dependent on the shock orientation1. Recent MHD-PIC simulations seem to support this notion (van Marle et al. 2018), and so we do not differentiate between quasi-parallel and quasi-perpendicular shocks. This injection scenario should not be taken literally, in particular not for electrons, for which pre-acceleration to a few tens of MeV is required and established at the shock (Matsumoto et al. 2017; Li et al. 2018; Bohdan et al. 2019). We are interested in particles at energies well above 100 MeV, and so the particulars of that pre-acceleration can be ignored.

It should be noted that ηi is not the only factor controlling injection. The shock velocity and the upstream density determine the rate with which particles pass through the shock. As the shock speed decreases, the injection rate decreases as well. While in the wind zone of a core-collapse supernova a large part of the particles is injected early in the evolution, for a type-Ia SNR expanding into a uniform ISM, the area-integrated injection rate increases with time, on account of the enlarging shock surface. As a result, the most recently injected particles dominate the volume-averaged particle spectra.

Usually, ηi is assumed to be constant, but it is not inconceivable that ηi might change with time, as it is a simple parametrization of nonlinear microphysical processes operating at the shock front (Völk et al. 2003; Petruk & Kopytko 2016). To investigate the effects of a time-dependent injection efficiency, we also consider a variable ηi of the form

η i , t ( t ) = η i ( t t 0 ) a . $$ \begin{aligned} \eta _{i,\mathrm{t}}(t) = \eta _i\left(\frac{t}{t_0}\right)^{a}. \end{aligned} $$(5)

2.2. Magnetic field and diffusion coefficients

To obtain the large-scale magnetic field, we assume it is dynamically irrelevant, and hence solve the induction equation following (Telezhinsky et al. 2013). The magnetic field is assumed to be constant in the upstream of the shock at a value of 5 μG. The field strength in the immediate downstream of the shock is 11 · 5 μ G 16 μ $ \sqrt{11}\cdot5\,\mu\mathrm{G}\approx16\,\mu $G2.

Observations indicate an amplification of the magnetic field in the downstream to several 100 μG for very young SNRs (Berezhko et al. 2003), at least part of which likely arises from MHD processes at and behind the shock (Giacalone & Jokipii 2007). The theory also suggests efficient magnetic-field amplification in the precursor of the remnant (Lucek & Bell 2000; Bell & Lucek 2001). However, the nonlinearity imposed by this magnetic-field amplification does not allow simple scaling (or averaging) of the results to the entire population of Galactic cosmic rays; modeling it must, therefore, be beyond the scope of this paper, and we checked magnetic turbulence to only reach levels of δB ≈ B0. The amplification of magnetic field well beyond its initial value imposes energy loss of CRs that can steepen the spectrum of the accelerated CRs. Bell et al. (2019) find a spectral steepening that scales with vshock/c and is hence only relevant during the early phases of SNR evolution. Here, this effect can be neglected, as the amplification factor of the field only reaches a value of ≈1, considerably less than assumed by Bell et al. (2019).

We choose two models for the diffusion that can be considered as limiting cases. In one case, we apply Bohm-like diffusion at the shock and have an exponential transition to the Galactic diffusion coefficient further upstream at 2Rsh (following Telezhinsky et al. 2012b).

In the second case, we solve, in parallel, a transport equation for the magnetic turbulence spectrum, assuming Alfvén waves only, and thus calculate the diffusion coefficient self consistently (Brose et al. 2016). In this case, the diffusion coefficient is time-dependent. However, we inject only a small number of particles, and thus the magnetic field is not amplified above the background magnetic field value in this case. The former case can be considered as a lower limit, and the latter ansatz as an upper limit for the diffusion coefficient.

In both cases, the diffusion process is assumed to be isotropic. In reality, D ≪ D, and so the escape of particles along the field lines is enhanced. This situation can be neglected as long as the coherence length of the magnetic field is smaller than the size of the remnant (López-Coto & Giacinti 2018). In this study, we are mainly interested in middle-aged and old remnants whose size is sufficient to justify the isotropy assumption.

2.3. Hydrodynamics

The evolution of an SNR without CR feedback can be described with the standard gas-dynamical equations:

t ( ρ m E ) + ( ρ v mv + P I ( E + p ) v ) T = ( 0 0 L ) $$ \begin{aligned}&\frac{\partial }{\partial t}\left( \begin{array}{c} \rho \\ \boldsymbol{m}\\ E \end{array} \right) + \nabla \left( \begin{array}{c} \rho \boldsymbol{v}\\ \boldsymbol{mv} + P\boldsymbol{I}\\ (E+p)\boldsymbol{v} \end{array} \right)^T = \left( \begin{array}{c} 0\\ 0\\ L \end{array}\right)\end{aligned} $$(6)

ρ v 2 2 + P γ 1 = E , $$ \begin{aligned}&\frac{\rho \boldsymbol{v}^2}{2}+\frac{P}{\gamma -1} = E, \end{aligned} $$(7)

where ρ is the density of the thermal gas, v the plasma velocity, m = vρ the momentum density, P the thermal pressure of the gas, L the energy losses due to cooling, and E the total energy density of the ideal gas with γ = 5/3. This system of equations is solved under the assumption of spherical symmetry in 1D using the PLUTO code (Mignone et al. 2007). The nonequilibrium cooling function, L, is taken from Sutherland & Dopita (1993).

In this work, we display results for type-Ia supernova explosions and discuss the basic differences to core-collapse explosions in Appendix A. Therefore, we initiate the simulations with ejecta profiles:

ρ SN = A exp ( v / v e ) t i 3 , and v = r / t i , $$ \begin{aligned}&\rho _\mathrm{SN} = A\exp (-v/v_{\mathrm{e}})t_i^{-3,}\,\mathrm{and}\, v = r/t_i, \end{aligned} $$(8)

with v e = ( E ex 6 M ej ) 1 / 2 , and A = 6 3 / 2 8 π M ej 5 / 2 E ex 3 / 2 $$ \begin{aligned}&\,\mathrm{with}\, v_{\mathrm{e}} = \left( \frac{E_\mathrm{ex} }{6M_\mathrm{ej} } \right)^{1/2,} \,\mathrm{and}\, A =\frac{6^{3/2}}{8\pi }\frac{M_\mathrm{ej} ^{5/2}}{E_\mathrm{ex} ^{3/2}} \end{aligned} $$(9)

as initial conditions (Dwarkadas & Chevalier 1998). Here ti = 20 yr is the start time of our simulation, Mej = 1.4 Msol the ejecta mass, Eex = 1051 erg the explosion energy, and r the spatial coordinate.

The initial age of about 20 years is rather large, but the solution quickly converges against solutions with a lower initial age. In any case, we are mainly interested in the later stages of the evolution. The density of the ambient medium was chosen to be 0.4 cm−3 unless stated differently.

The results of the hydro simulation for the density, velocity, pressure, and temperature distributions are then mapped onto the spatial coordinate of the CR and turbulence grid, respectively. The shock, that is typically a few bins wide in the hydro-solution, needs to be resharpened in order to guarantee a realistically high acceleration rate from GeV to TeV energies. This procedure is repeated for each time step of the CR and turbulence grid. One of these time steps typically requires many time steps of the hydro solver.

3. Results

Using the framework described above, we followed the evolution of the remnant for 100 000 years. The shock speed was 9200 km s−1 after 100 yrs. The transition to the Sedov phase happened after 1300 years when the swept-up mass was approximately 10 Msol and the shock speed vshock = 2350 km s−1. The remnant reached the post-adiabatic phase after 35 000 years3 with vshock = 300 km s−1. After roughly 85 000 years (vshock = 130 km s−1), the shock compression ratio started to fall significantly below 4, and, at the end of the simulation, it had a value of 3.15. Figure 1 displays the time evolution of the speed and position of the forward shock (cf. Fig. 4a by Petruk et al. 2016).

thumbnail Fig. 1.

Curve shows the radial position of the forward shock as a function of time. The transitions between free expansion and the Sedov phase, and between Sedov and the post-adiabatic phase, appear roughly at 1.3 kyr and 35 kyr, respectively.

Already after 10 000 years, the shock speed is down to 640 km s−1, and so the difference in speed between the end of the Sedov phase and the post-adiabatic phase is only a factor of ten. The maximum energy, that can be reached during the post-adiabatic phase should then be only one order of magnitude lower than that during the Sedov phase.

In this section, we first present results for a post-adiabatic remnant under different assumptions about CR diffusion. We describe the escape of CRs from the remnant, the reacceleration of pre-existing CRs, and the inverse Compton and Pion-decay gamma-ray spectra.

3.1. Escape

To distinguish between the escape of CRs produced in the remnant and reacceleration of Galactic CRs, in this section, we set the density of Galactic CRs to zero.

We first evaluate the evolution of the CR-spectra for the two diffusion scenarios that we introduced in Sect. 2.2 – Bohm-like diffusion in a 5 μG field around the remnant and a diffusion coefficient obtained from the self-consistent amplification of Alfvénic turbulence.

Figure 2 shows the evolution of the volume-averaged cosmic-ray spectrum in the downstream region at three points in time that represent the three stages of SNR evolution. These downstream spectra represent the particle population that is mainly responsible for the detectable emission from the remnant.

thumbnail Fig. 2.

Volume-averaged downstream proton spectra for Bohm-like (red) and self-consistent (green) diffusion. The background of Galactic CRs is neglected.

What should be noted from the figure is that for Bohm-like diffusion (red lines), the spectra are power laws with the typical test-particle index of s = 2.0 at all stages of the evolution. The normalization decreases as the remnant decelerates and particles get injected at lower energies, where most particles reside. The maximum energy decreases approximately from 20 TeV to 5 TeV. The decrease in maximum energy is small, because for a constant upstream magnetic-field strength, the maximum energy decreases very slowly in the Sedov phase; in fact the shock speed only drops by a factor 20 between the beginning of the Sedov phase and the sharp decrease of the shock-velocity at an age of 85 000 yr, and the time available for particle acceleration at the late, slow shock is roughly one order of magnitude longer than the lifetime of the fast shock.

With explicit treatment of turbulence transport and Alfvénic diffusion, the decreasing normalization of the CR density and hence the CR pressure gradient reduces the driving of turbulence. As a consequence, the diffusion coefficient and the timescale for acceleration increase, and the maximum energy falls from 5 TeV to 10 GeV. As it takes time for particles to escape from the remnant, high-energy particles are still present in small numbers, and the cosmic-ray spectra display a break at the momentary maximum energy and are soft at higher energies. Figure 4 shows the total number of CRs for the Alfvénic scenario at different times. It is clear that all CRs with an energy above 1 TeV are produced within the first 10 kyr of the remnants evolution. At the later stages, CRs from the downstream escape to the upstream, which forms softer spectra inside the remnant. The break in the spectrum occurs at the energy the SNR is currently able to accelerate CRs to. This behaviour was also obtained in the analytic calculations of Celli et al. (2019), who found that high-energy particles that already escaped the acceleration process can be trapped inside and close to the remnant, and they introduce a spectral break in the spectrum of particles present inside the SNR. The particle spectrum above the break energy is not a simple power law, but can be reasonably well approximated with a s = 2.7 power-law index. Interestingly, this spectral index is in rough agreement with that measured in IC443 and W44 (Ackermann et al. 2013). In both remnants, the CR spectra also feature a spectral break around a few hundred GeV. We acknowledge that both IC443 and W44 interact with dense material, and locally they have a wide range of evolutionary ages. Our model assumes spherical symmetry and an external medium with a constant density. In first-order approximation, a composite model based on spectra calculated for different ages should permit a rough comparison with interacting SNRs, though.

Our findings are compatible with analytic calculations by Malkov et al. (2012). There, the authors suggested that an evanescence of Alfvén waves due to strong neutral-charged collisions in SNRs interacting with dense environments leads to a spectral steepening by exactly one power (Δs = 1.0). In this case, particles above a break energy around 10 GeV can escape the shock precursor, and the spectra above the break energy feature softer spectra. The authors noted that their model explains the gamma-ray observations of W44 (Ackermann et al. 2013; Giuliani et al. 2011) well, but it overpredicts the spectral steepening observed in IC443 and W28. Here, the spectral index changes only by Δs = 0.6 − 0.7 (Ackermann et al. 2013; Acciari et al. 2009; Abdo et al. 2010), which is in agreement with our findings.

In our model, the reason for the softening of our spectra is also an evanescence of Alfvén waves, but it is caused by reduced driving of turbulence at later stages of the remnant’s evolution, which introduces additional nonlinearity and time-dependence. A second difference to the scenario of Malkov et al. (2012) is that escaping particles are still scattered and confined to the shock region in our simulations.

We can also evaluate the total production spectrum of CRs, including those that left the remnant. Our spatial grid extends to 65 Rshock, and all accelerated particles stay in the simulation domain. Then, the CR spectrum integrated over the whole simulation domain represents the total CR yield of the SNR. These total production spectra then represent the spectrum of CRs that the remnant injects into our galaxy4. Still, spectra above the current maximum energy are only weakly modified during later times of the SNR’s evolution, where they are subject to further modification during galactic propagation. These spectra are typically different from the downstream spectra that represent the particle population responsible for the emission detected from the remnant. Figure 3 shows the total production spectra for the two diffusion regimes discussed in this section. Again, for Bohm-like diffusion, the spectrum is well-represented by a power law with index s = 2.0. In the Alfvénic regime, a spectral break occurs at ≈10 GeV after 50 000 yr, where the spectral index changes approximately from s = 2.0 to s = 2.4.

thumbnail Fig. 3.

Total volume-averaged proton spectrum in the case of Bohm-like (red) and self-consistent (green) diffusion. Galactic CRs are ignored.

thumbnail Fig. 4.

Proton number-spectra for the Alfvénic diffusion scenario at different times. Solid lines represent the total, dashed the downstream and dotted the upstream spectra.

The total proton spectrum coincides well with that deduced from fits of cosmic-ray data with Galactic propagation models, where the required injection spectrum tends to be slightly harder than E−2 below a few GeV, transitioning to a softer index s ≃ 2.4 above 10 GeV (Putze et al. 2009; Trotta et al. 2011). Our results satisfy that expectation for the entire energy range above 10 GeV where direct measurements are possible, and solar modulation is weak or absent, even though we conduct test-particle calculations for which the spectra of CRs at the shocks are always standard DSA s = 2.0 spectra. They also provide a natural explanation for the break in the cosmic-ray source spectrum that was introduced ad-hoc in the fits of the cosmic-ray data and appears here on account of the explicit treatment of turbulence driving upstream of the shock. All changes in the spectral index arise from the time evolution of SNRs, the turbulence spectra, and consequently the maximum energy of accelerated particles. Simple scaling relations such as the assumption of Bohm diffusion significantly modify the model expectations, as do steady-state descriptions of the turbulence level and the diffusion coefficient. The spectral indices obtained for Alfvénic diffusion depend on the particulars of wave growth, damping, and cascading. They reflect the balance between the reduction rate of the maximum energy of CRs at the shock and the escape rate of high-energy CRs from the SNR (Ohira et al. 2010).

Similar soft spectra to those caused by the decay of turbulence can be obtained, if the maximum particle energy significantly varies along the shock surface. Malkov & Aharonian (2019) discuss the expansion of an SNR in a homogeneous magnetic field, assuming that efficient acceleration only takes place at two polar caps, as the shock normal has to be sufficiently parallel to the magnetic field for DSA to operate. The spectra observed along some lines of sight through the acceleration region will then sample regions with different maximum energies and be softer than those predicted by standard DSA. Both methods rely on a strong variation of the maximum energy, either spatially or temporally, and the softer spectra arise from superposition of spectra at different locations or times.

Our total production index above the break energy differs from the results obtained by Celli et al. (2019). Their model predicts a spectral index of s = 2 for the total production spectrum as the spectrum we obtain through acceleration at the shock is always a s = 2 power-law. However, they assumed that a constant fraction of the SNRs kinetic energy is transferred into CRs. As a result, the fraction of thermal particles injected as CRs is not a constant in their calculation and depends on the shock ram pressure and the maximum particle energy5. The effect of an increasing shock surface – meaning more particles injected into the acceleration process at later times – is canceled by the impact on the proton spectral normalization of the decreasing shock speed, leading to s = 2. In our case, the fraction of thermal particles injected as CRs is a constant that leads to a different time dependence of the spectral normalization. Thus, our resulting spectra can be softer than s = 2, even if our acceleration spectrum has the standard DSA spectral index.

Figure 5 shows the distribution of the spectral index of the CR electrons compared to that of the protons. What should be noted from the figure is that beyond an age of about 300 yr, the electron spectra are softer than the proton spectra on account of synchrotron cooling, which is in agreement with the findings of Diesing (2019). However, the difference in spectral index is a strong function of the turbulence model, in particular its time dependence, and cannot be described as a simple shift in the spectral index. We emphasize again, that the decay of turbulence and the subsequent escape of particles is essential for the formation of soft particle spectra and spectral breaks, even if the acceleration mechanism at the shock produces a standard DSA spectrum. Steady-state models and simple scaling relations lack some of these features.

thumbnail Fig. 5.

Running power-law index of the electron (black) and proton (red) spectrum as function of energy displayed at different ages of the remnant. The Bohm-like diffusion scenario is shown in the top panel, and the Alfvénic-turbulence scenario in the bottom panel.

3.2. Reacceleration

It has been proposed by several authors (e.g., Tang & Chevalier 2015; Cardillo et al. 2016; Caprioli et al. 2018) that a part of the particles accelerated at an SNR might be reaccelerated, preexisting CRs. They may even dominate the gamma-ray output of W44 and IC443 (Ackermann et al. 2013).

To test this scenario, we included background CRs into our simulations. They should become relevant when the injection rate of particles into DSA at the shock is very low. To numerically control this threshold effect, we artificially varied the injection efficiency over time according to Eq. (5). We adopted a = −1.0 for the power-law dependence of the injection with time, ηi ∝ ta, and took 20 years as our reference point. After 2500 years, we injected with an efficiency that is 0.8% of that used initially. From this point in time, the normalization of the CRs injected at the shock is low enough that background CRs can dominate the acceleration and emission. The level of background CRs that we have in our simulation domain was negligible for the results presented in the previous section. It has been argued that the level of background CRs might be increased inside the molecular clouds that some remnants are interacting with. However, the possible level of background CRs inside the clouds is constrained by observations of the gamma-ray emission from these clouds and found to be comparable to that measured outside of the clouds (Yang et al. 2014).

Again, we investigate the two regimes of Bohm-like and self-generated Alfvénic diffusion. All other simulation parameters are kept the same as in the previous section. Figure 6 illustrates the CR spectra inside the remnant for the two diffusion models for the three phases of the SNR’s evolution.

thumbnail Fig. 6.

Volume-averaged downstream proton spectra for Bohm-like (red) and self-consistent (green) diffusion at different times. We include background CRs and artificially decrease, with time, the injection fraction of thermal particles at the shock. The CR-density at 300 yr for nondecreasing injection (direct) is shown as a dotted gray line for comparison.

Again, for Bohm-like diffusion, the spectrum of accelerated CRs remains a s = 2.0 power law with only a modest change in the maximum energy. After 2500 yr, the background CRs fill the entire simulation volume, and the normalization of the CR-spectrum does not decrease further than it did without background CRs (cf. Fig. 2). Additionally, the low-energy part of the spectrum is dominated by freshly injected CRs, as our imposed background CR-spectrum is very hard below 1 GeV (see Eq. (2)), and freshly injected CRs dominate the spectrum up to ≈100 MeV. However, the shape of the low-energy part of the spectrum strongly depends on the parametrization of the CR spectrum below 1 GeV, which is not directly measurable due to the effects of solar modulation.

The situation looks similar for self-generated Alfvénic turbulence. The main difference is, that the accelerated background CRs provide too little turbulence to be contained by the remnant. A higher level of background CRs would enhance the amplification of turbulence, and thus the confinement of CRs, but a much higher background flux is in contradiction with direct observations (Yang et al. 2014). The particle spectrum extends to about 100 GeV before it cuts off with a spectral index s >  −2.7. This break energy coincides with that deduced for W44 and IC443 (Ackermann et al. 2013).

In fact, both reacceleration of CRs from the sea of background CRs and the escape of CRs from radiative remnants provide similar signatures at higher energies. Both produce a spectral break at energies between 10 and 100 GeV with soft spectra at higher energies, with a slightly stronger turnover in the case of reacceleration. The emission signature will differ mostly at radio energies where either freshly injected CRs or background CRs dominate, which yield different spectra as shown in Fig. 7.

thumbnail Fig. 7.

Ratio of the radio flux produced by directly accelerated and reaccelerated electrons, normalized above 75 GHz. The assumed magnetic-field strength in the downstream region is 16 μG.

If freshly accelerated electrons dominate, the radio spectrum will be a featureless power law with index α = 0.5. The situation is different for reaccelerated electrons, as the diffuse Galactic radio emission indicates electron spectra harder than s = 2.0 below 4 GeV. This hard spectrum is retained when electrons get reaccelerated, whereas the parts of the spectrum that are softer than s = 2.0 change to s = 2.0 during the reacceleration process. If reaccelerated electrons are supposed to be more important than freshly accelerated electrons above a few GeV, then we would observe freshly accelerated electrons at very low energies with a p−2 spectrum, until reaccelerated electrons take over with a harder spectrum reflecting that of electrons in the ISM, until around a few GeV the electron spectrum turns over to p−2 again. The latter transition would appear in the radio synchrotron spectra right in the observable frequency band between 100 MHz and 50 GHz, depending on the strength of the magnetic field. We note that remnants that are discussed to reaccelerate CRs, like W44, show remarkably featureless power-law radio spectra (Castelletti et al. 2007).

Cristofari & Blasi (2019) recently argued that the gamma-ray emission from the SNRs SN1006, Vela Jr. and RX J1713.7−3946 might be explained in terms of IC-emission from reaccelerated galactic electrons without the need for a contribution from freshly accelerated electrons. Their ansatz and assumptions for the electron background spectrum are similar to those used in this paper, but they did not investigate the radio emission that would arise from the reacceleration of galactic electrons. As the magnetic fields in these three remnants are likely between 10 μG and 100 μG and thus comparable to the 16 μG we assumed in Fig. 7 (H.E.S.S. Collaboration 2018a; Sushch et al. 2018; Berezhko et al. 2003), the radio spectra of these remnants should feature turnovers described above. While the radio data for Vela Jr. and RX J1713.7−3946 are too sparse to draw any conclusions, it is clear that there is no evidence for a softening in the radio spectrum of SN1006; if anything, there might be a hardening with increasing energy (Allen et al. 2008).

Figure 8 shows the CR pressure in our simulations. The CR pressure stays well below 10% of the shock ram pressure (gray line), and so our test-particle assumption is valid most of the time. Only for Bohm-like diffusion do CRs become dynamically important at the end of the evolution, when the SNR enters the dissipative phase. It is likely that at this point the magnetic field is also dynamically relevant (Petruk et al. 2016).

thumbnail Fig. 8.

Evolution of CR pressure over time for the two diffusion scenarios and with/without the Galactic CR background. The gray line indicates 10% of the shock ram pressure.

Background CRs change the time evolution of the CR pressure, which starts to fall less rapidly as soon as the reaccelerated background CRs begin to dominate the particle distribution at the shock. The initial rapid decline of the CR pressure in simulations with background CRs arises from our choice of a decreasing injection efficiency, ηi ∝ t−1, in these simulations. In any case, CRs start to be dynamically important at very late time, as also suggested by Simpson et al. (2016).

3.3. Illuminated clouds

Escaping cosmic rays do not only replenish CRs in the Galaxy, but may also illuminate molecular clouds close to their mother SNRs, leading to intense gamma-ray emission (Gabici & Aharonian 2007). This illumination effect is a clear indicator of the acceleration of hadrons and has been observed in W28 (Cui et al. 2018). The spectra of escaping CRs are typically different from those of CRs inside the remnant. Figure 9 illustrates the spectra of escaped CRs 7.5 pc from the SNR shock in a 2 pc thick region at an age of 37 kyr for the four setups (Bohm-like/Alfvénic diffusion, acceleration/reacceleration) discussed before.

thumbnail Fig. 9.

CR spectra about 7.5 pc ahead of the forward shock at an age of 37 kyr. Red (black) lines refer to Bohm-like diffusion without (with) background CRs, and green (blue) lines stand for Alfvénic diffusion without (with) background CRs. Note that the y-axis is scaled with p2.75.

The spectra for Bohm-like diffusion resemble narrow bumps that can be described as log-parabolas. If the maximum energy of accelerated particles evolves quickly, as it does for diffusion in the time-dependent Alfvénic turbulence, the spectra are much wider and have a lower amplitude. Over a wide range, the far upstream spectrum of reaccelerated CRs is harder than s = 2.75 – which corresponds to a horizontal line in Fig. 9 – with a cutoff near 10 TeV. Freshly accelerated CRs have a higher flux and a roughly symmetric spectrum. Apart from the normalization and the slightly different spectral index, there is no spectral signature that would clearly distinguish reaccelerated CRs from particles injected from the thermal pool and accelerated at the shock.

3.4. Gamma-ray emission

At least in the TeV band, it is still unclear which emission process – inverse-Compton radiation or the decay of neutral pions – dominates the gamma-ray emission from SNRs. We evaluated the time evolution of the very-high energy (VHE) luminosity for both emission processes over the entire lifetime of an SNR and compared it to data of the H.E.S.S. Galactic SNR survey (H.E.S.S. Collaboration 2018b). Figure 10 shows our prediction for the time evolution of the gamma-ray luminosity for inverse-Compton (IC) and Pion-decay (PD) emission in the 1 − 10 TeV band for Bohm-like and Alfvénic diffusion.

thumbnail Fig. 10.

Time evolution of the gamma-ray luminosity in the TeV band. The curves reflect models with Bohm-like (red) and Alfvénic (green) diffusion, and they indicate IC (solid lines) and PD (dashed lines) emission. The curves are calculated for an ambient gas density of 0.4 cm−3 assuming the same injection-efficiency for electrons and protons.

Our models are calculated for a ambient density of 0.4 cm−3 and assume the same injection efficiency for electrons and protons6. The contribution of background CRs is ignored in this section.

There is a fundamental difference in the behaviour of IC and PD emission at later times that arises from the energy loss of electrons via synchrotron radiation, which suppresses the IC emission in the VHE-band after roughly 3000 yr. The maximum energy of the electrons is given approximately by a balance between the energy-loss rate and the acceleration rate,

E max 25 TeV ( v sh 1000 km s 1 ) 5 μ G B 0 , $$ \begin{aligned} E_\mathrm{max} \approx 25\,\mathrm{TeV}\, \left(\frac{v_\mathrm{sh} }{1000\,\mathrm{km\,s^{-1}}}\right)\sqrt{\frac{5\,\mu \mathrm{G}}{B_0}}, \end{aligned} $$(10)

where me, q, vsh, and B0 denote the electron mass, the elementary charge, the shock speed, and the upstream magnetic-field strength, respectively7. If the shock is slower than 2000 km s−1, the maximum electron energy falls below 50 TeV, and hence the IC-cutoff energy approaches 1 TeV. Our estimate suggests that the IC flux in the TeV band should start to decrease at an age of 3600, 1700, and 750 years for low (n = 0.04 cm−3), medium (n = 0.4 cm−3), and high (n = 4.0 cm−3) density of the ambient medium, respectively. This roughly fits to our simulations in which the actual peak luminosity is reached later in all three cases, as the expansion, and thus the increase in the number of radiating particles initially compensates for the decreasing maximum energy.

Protons do not efficiently lose energy and can produce VHE gamma rays throughout the lifetime of the SNR. Thus, the hadronic gamma-ray brightness keeps increasing with time for Bohm-like diffusion. In contrast, for Alfvénic diffusion, the faster escape of CRs leads to a roll-off in the gamma-ray luminosity that is faster for leptonic emission and lower for hadronic channels but evident in both. The luminosity peak is reached at earlier times than for Bohm-like diffusion, as any weaker driving of turbulence limits the acceleration efficiency.

The H.E.S.S. Collaboration published a study examining the VHE luminosity of eight detected and several undetected SNRs with known distances, ages, and ambient densities (H.E.S.S. Collaboration 2018b, and references therein)8 In Fig. 11, we compare our results for Bohm-like diffusions with these observations. To account for the different densities observed for these remnants, we ran two additional models with ambient densities of 0.04 cm−3 and 4.0 cm−3. All other parameters remained unchanged.

thumbnail Fig. 11.

Observed gamma-ray luminosity of galactic SNRs in the TeV band; blue markers indicate upper limits. The red curves reflect the models with Bohm-like diffusion and indicate inverse-Compton (solid lines) and pion-decay (dashed lines) emission. The thick curves are calculated for a ambient gas density of 0.4 cm−3. The shaded uncertainty bands represent the range of luminosity expected for a density in the range 0.04 − 4.0 cm−3. The references for the data points are given in the text.

To reproduce the luminosity of the detected SNRs, we had to lower the injection efficiency for electrons and protons, and we kept it the same for all ambient densities. The electron-to-proton ratio is 0.025 m e / m p $ 0.025\sqrt{{m_{\mathrm{e}}}/{m_{\mathrm{p}}}} $ at high energies.

High ambient densities result – as expected – in a much brighter PD emission, but suppress the IC emission at the same time on account of the rapid decline of the shock speed. Thus, the maximum electron energy already starts to decrease after one kyr as opposed to 2 − 4 kyr in the medium and low-density cases. Low ambient densities thus allow for a higher IC peak luminosity as synchrotron cooling starts to decrease the maximum electron energy only after a larger total number of electrons got accelerated, compared to the medium and high-density cases.

Initially, the density-dependence of the IC emission is much weaker than that of the hadronic emission. The VHE luminosity is proportional to the number of CR particles, which for a fixed injection efficiency is only weakly dependent on the ambient density. In the adiabatic phase, one finds, for the enclosed volume, V ρ 0 3/5 $ V\propto \rho_0^{-3/5} $. The total number of particles injected into the acceleration process is thus ρ 0 2/5 $ {\propto} \rho_0^{2/5} $. However, the PD-luminosity additionally depends on the ambient density, and the resulting scalings L IC ρ 0 2/5 $ L_{\rm IC} \propto \rho_0^{2/5} $ and L PD ρ 0 7/5 $ L_{\rm PD} \propto \rho_0^{7/5} $ are roughly observed in our simulations.

Our findings naturally explain why the brightest observed SNRs (Vela Jr., RX J1713) expand in very low density environments and remnants with similar ages – like Puppis A, which expands in a much denser medium – have been missed by the IACTs (H.E.S.S. Collaboration 2015).

We note that our predictions in Fig. 11 are only valid up to SNR ages of roughly 4 kyr, as measurements indicate that the diffusion coefficients in the precursor of these remnants start to become higher than Bohm-like (Sushch et al. 2018). In our Alfvénic scenario, we already observe this transition after about 2 kyr. Afterward, the escape of high-energy particles from the far-downstream region becomes significant. As a consequence, the PD luminosity should start to fall as the particle spectra get softer at high energies (see Fig. 10). We do observe a few bright SNR with ages around 104 yr. Their soft spectra are in line with the Alfvénic-turbulence model, but their fluxes fit better to the Bohm predictions. This indicates that interaction with dense clouds, additional turbulence amplification mechanisms, or other processes may play a significant role.

There are additional uncertainties like the filling factor or the type of supernova explosion that we have not taken into account in this simple model. For SN1006, it is for instance known that the nonthermal emission has a bipolar structure, indicating that just a fraction of the shock-surface is accelerating CRs. Thus, the brightness of an SNR can easily be a factor of a few below our model predictions. Further, most SNRs originate from core-collapse SNs, and these remnants will, at least initially, expand in a density with a r−2-dependence. Here, the total number of particles will be higher by about a factor three than predicted by measurements of the present-day post-shock density. For the density profiles:

ρ ( r ) = { ρ u for Type 1 a ρ u ( R sh / r ) 2 for CC , $$ \begin{aligned} \rho (r) ={\left\{ \begin{array}{ll} \rho _u\qquad&\mathrm{for}\,\,\mathrm{Type}\,\,\mathrm{1a} \\ \rho _u\left({R_{\rm sh}}/{r}\right)^2\qquad&\mathrm{for}\,\,\mathrm{CC}, \end{array}\right.} \end{aligned} $$(11)

the total number of particles passed through the shock is

N tot = 4 π ρ r 2 v sh d t = { 4 3 π R sh 3 ρ u for Type 1 a 4 π R sh 3 ρ u for CC . $$ \begin{aligned} N_\mathrm{tot}&= 4\pi \int \rho r^2 v_{\rm sh} \mathrm{d} t\nonumber \\&= {\left\{ \begin{array}{ll} \frac{4}{3}\pi R_{\rm sh}^3 \rho _u \qquad&\mathrm{for}\,\,\mathrm{Type}\,\,\mathrm{1a}\\ 4\pi R_{\rm sh}^3 \rho _u \qquad&\mathrm{for}\,\,\mathrm{CC}. \end{array}\right.} \end{aligned} $$(12)

We discuss in detail the differences between type-Ia and core-collapse SNR in Appendix A. Similar variations arise when the shock recently encountered a density jump, like the edge of a wind-blown bubble. In that case, the total number of particles would be lower than expected on the grounds of density measurements, and the remnant would consequently be dimmer. Kepler’s SNR might be in this situation, which indicates that dedicated modeling for each remnant is still needed.

4. Conclusions

We developed a model for the particle acceleration in SNRs up to a point where the forward shock of the remnant is radiative by solving the time-dependent transport equation for the CRs, together with the gas-dynamical equations for the plasma flow. We investigated two diffusion regimes – Bohm-like and Alfvénic diffusion. For the Alfvénic diffusion, we solved the time-dependent transport equation for the evolution of the magnetic turbulence together with the gas-dynamical and the CR equations. The ratio of CR to shock ram pressure increased in all scenarios, but remained well below 10% at all times, except for the very end of our Bohm-diffusion simulation. Hence our simulations were conducted in the test-particle limit.

We showed that inefficient confinement in the Alfvénic diffusion regime of high-energy particles leads to a rapid reduction in the maximum energy reachable at the shock, and to the eventual formation of a spectral break around 10 − 100 GeV. The spectra of CRs inside the SNR show no simple power-law structure above the break, but can be reasonably well-approximated by a power law with an index of s = 2.7. This spectral structure is similar to that observed from W44 and IC443.

The evaluation of the total production spectra, including CRs that reside outside of the SNR, showed also broken power law spectra with a spectral index of s ≈ 2.4 at high energies. This is in rough agreement with the injection spectra required by galactic propagation models, namely around s = 2 at low energies and a break to s ≃ 2.4 around 10 GeV. The spectral index obtained for the release spectrum is thus harder than the spectrum of the CRs that still reside inside the SNR.

We investigated the reacceleration of galactic CRs at the shock. For efficient injection from the thermal pool, reaccelerated CRs would be too few to be noted, even for core-collapse supernovae. To render reaccelerated CRs visible, we artificially reduced the thermal injection of CRs at the shock. For an injection efficiency scaling as η ∝ t−1, galactic CRs became the dominant source for particles accelerated at the shock after ≈2500 yr. We showed that the resulting spectra above a few GeV are similar to those obtained by the acceleration of particles injected from the thermal pool. As a result, the emission signature at high energies is indistinguishable from that in the direct-acceleration scenario. However, the harder spectra of galactic electrons below a few GeV leave an imprint in the radio emission, as standard DSA can not soften spectra with an spectral index below s = 2 but hardens spectra with a spectral index larger than s = 2. Thus, the resulting radio spectra show a hardening below a few GHz – a feature that could be observable by radio and infrared telescopes. The emission resulting from CRs leaving the SNR and possibly illuminating a molecular cloud ahead of the shock shows narrow log-parabola spectra if the diffusion is Bohm-like, but wide spectra if the diffusion is Alfvénic. Again, the difference between directly accelerated and reaccelerated CRs offers few, if any, options to distinguish between both scenarios.

We examined the luminosity evolution of type-Ia SNRs assuming that the VHE gamma rays originate either from pion decay or inverse-Compton radiation. For Bohm diffusion, we found that the IC-luminosity is initially only weakly density-dependent in contrast to the PD luminosity, which scales roughly linearly with the ambient density. The total number of CRs accelerated in an SNR is proportional to the remnant’s size and the ambient density. However, the higher number of particles available for acceleration in high-density environments is largely compensated by the smaller size of the remnants compared to low-density environments. As a result, the density-dependence of the gamma-ray luminosity is roughly the same as the density dependence of the emission process itself.

In the case of Bohm-like diffusion, the remnant’s IC luminosity peaks at an age around a few thousand years, after which cooling suppresses the emission at high energies, while the PD luminosity keeps increasing with time. A study of different ambient densities revealed that SNRs expanding in low-density environments show a higher maximum IC-luminosity in accordance with observations – all the brightest VHE gamma-ray remnants RCW86, Vela Jr, HESS J1731−347 and RX J1713.7−3946 – expand in low-density environments. For a low ambient density, the shock speed remains high for a long time, and acceleration remains faster than synchrotron losses for electrons.

If the diffusion coefficient is time-dependent as in the Alfvénic scenario, the PD-emission peaks as well. In this case both luminosities reach their maximum after roughly 1 kyr. Afterwards, cooling for electrons or the escape of high-energy protons reduce the VHE luminosity. This could explain why old SNRs are only detected in the VHE gamma ray when interacting with massive molecular clouds, as the high target density for PD-emission significantly enhances their luminosity. However, the position of the peak indicates that either the interaction with the dense material began late in the evolution of the SNR or is otherwise restricted, or a more efficient amplification of turbulence needs to be considered.

This study and thus its conclusions are based on the modeling of the type-Ia SNRs. However, we did also model core-collapse SNRs evolving in the wind zone with the spatially dependent magnetic field and found that results are comparable to type-Ia. The initially higher magnetic field might be important, as it causes severe synchrotron losses, but later on after around 500 years as the remnant expands into the shocked wind region, the evolution of the core-collapse SNR basically mimics the evolution of the type-Ia SNR due to the roughly constant density and magnetic field.


1

It is merely conjecture that the morphology of the nonthermal emission from SN1006 is often attributed to the local magnetic-field angle (Völk et al. 2003).

2

We consider a magnetic field with equally strong parallel and perpendicular components. The parallel direction is not compressed at the shock, hence the compression-factor is 11 $ \sqrt{11} $ instead of 4.

3

This is the time t* when radiative losses of a fluid element during t* are comparable to its initial thermal energy (Cox 1972), so the shock cannot be considered adiabatic anymore.

4

If all acceleration stopped at this given time so the particle spectrum would not be modified any more.

5

Their spectral normalization is ξ CR ρ 0 v sh 2 $ {\propto} \xi_{\rm CR} \rho_0 v_{\rm sh}^2 $

6

The same injection efficiency will result in the same total number of electrons and protons in the simulation domain and an electron-to-proton ratio of K ep m e / m p $ K_{\mathrm{ep}}\simeq\sqrt{{m_{\mathrm{e}}}/{m_{\mathrm{p}}}} $ at relativistic energies (Pohl 1993).

7

The formula was derived assuming Bohm diffusion, a magnetic-field compression at the shock by 11 , $ \sqrt{11,} $ and only considering the upstream diffusion coefficient as important for the acceleration time.

8

The flux for G69.7 is taken from Abeysekara et al. (2018), and a conservative distance of 7.1 kpc was adopted (Kilpatrick et al. 2016).

Acknowledgments

The authors gratefully thank the referee Mikhail Malkov for his constructive input.

References

  1. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 718, 348 [NASA ADS] [CrossRef] [Google Scholar]
  2. Abeysekara, A. U., Archer, A., Aune, T., et al. 2018, ApJ, 861, 134 [NASA ADS] [CrossRef] [Google Scholar]
  3. Acciari, V. A., Aliu, E., Arlen, T., et al. 2009, ApJ, 698, L133 [NASA ADS] [CrossRef] [Google Scholar]
  4. Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  5. Allen, G. E., Houck, J. C., & Sturner, S. J. 2008, ApJ, 683, 773 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bell, A. R., & Lucek, S. G. 2001, MNRAS, 321, 433 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bell, A. R., Matthews, J. H., & Blundell, K. M. 2019, MNRAS, 488, 2466 [NASA ADS] [CrossRef] [Google Scholar]
  9. Berezhko, E. G., Ksenofontov, L. T., & Völk, H. J. 2003, A&A, 412, L11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Blasi, P., Gabici, S., & Vannoni, G. 2005, MNRAS, 361, 907 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bohdan, A., Niemiec, J., Pohl, M., et al. 2019, ApJ, 878, 5 [NASA ADS] [CrossRef] [Google Scholar]
  12. Brose, R., Telezhinsky, I., & Pohl, M. 2016, A&A, 593, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Bykov, A. M., Ellison, D. C., Osipov, S. M., & Vladimirov, A. E. 2014, APJ, 789, 137 [NASA ADS] [CrossRef] [Google Scholar]
  14. Caprioli, D., Blasi, P., Amato, E., & Vietri, M. 2009, MNRAS, 395, 895 [NASA ADS] [CrossRef] [Google Scholar]
  15. Caprioli, D., Zhang, H., & Spitkovsky, A. 2018, J. Plasma Phys., 84, 715840301 [CrossRef] [Google Scholar]
  16. Cardillo, M., Amato, E., & Blasi, P. 2016, A&A, 595, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Castelletti, G., Dubner, G., Brogan, C., & Kassim, N. E. 2007, A&A, 471, 537 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Celli, S., Morlino, G., Gabici, S., & Aharonian, F. 2019, MNRAS, 490, 4317 [NASA ADS] [CrossRef] [Google Scholar]
  19. Cox, D. P. 1972, ApJ, 178, 159 [NASA ADS] [CrossRef] [Google Scholar]
  20. Cristofari, P., & Blasi, P. 2019, MNRAS, 489, 108 [NASA ADS] [CrossRef] [Google Scholar]
  21. Cui, Y., Yeung, P. K. H., Tam, P. H. T., & Pühlhofer, G. 2018, ApJ, 860, 69 [Google Scholar]
  22. Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18 [Google Scholar]
  23. Diesing, R. 2019, Int. Cosmic Ray Conf., 36, 59 [Google Scholar]
  24. Diesing, R., & Caprioli, D. 2019, Int. Cosmic Ray Conf., 36, 59 [Google Scholar]
  25. Dwarkadas, V. V., & Chevalier, R. A. 1998, ApJ, 497, 807 [NASA ADS] [CrossRef] [Google Scholar]
  26. Gabici, S., & Aharonian, F. A. 2007, ApJ, 665, L131 [NASA ADS] [CrossRef] [Google Scholar]
  27. Giacalone, J., & Jokipii, J. R. 2007, ApJ, 663, L41 [NASA ADS] [CrossRef] [Google Scholar]
  28. Giuliani, A., Cardillo, M., Tavani, M., et al. 2011, ApJ, 742, L30 [NASA ADS] [CrossRef] [Google Scholar]
  29. H.E.S.S. Collaboration (Abramowski, A., et al.) 2015, A&A, 575, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. H.E.S.S. Collaboration (Abdalla, H., et al.) 2018a, A&A, 612, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. H.E.S.S. Collaboration (Abdalla, H., et al.) 2018b, A&A, 612, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Hanusch, A., Liseykina, T. V., & Malkov, M. 2019, ApJ, 872, 108 [NASA ADS] [CrossRef] [Google Scholar]
  33. Jaffe, T. R., Banday, A. J., Leahy, J. P., Leach, S., & Strong, A. W. 2011, MNRAS, 416, 1152 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kilpatrick, C. D., Bieging, J. H., & Rieke, G. H. 2016, ApJ, 816, 1 [NASA ADS] [CrossRef] [Google Scholar]
  35. Li, R., Zhou, C. T., Huang, T. W., et al. 2018, Phys. Plasmas, 25, 082103 [NASA ADS] [CrossRef] [Google Scholar]
  36. López-Coto, R., & Giacinti, G. 2018, MNRAS, 479, 4526 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lucek, S. G., & Bell, A. R. 2000, MNRAS, 314, 65 [NASA ADS] [CrossRef] [Google Scholar]
  38. Malkov, M. A. 1998, Phys. Rev. E, 58, 4911 [NASA ADS] [CrossRef] [Google Scholar]
  39. Malkov, M. A., & Aharonian, F. A. 2019, ApJ, 881, 2 [NASA ADS] [CrossRef] [Google Scholar]
  40. Malkov, M. A., Diamond, P. H., & Sagdeev, R. Z. 2011, Nat. Commun., 2, 194 [NASA ADS] [CrossRef] [Google Scholar]
  41. Malkov, M. A., Diamond, P. H., & Sagdeev, R. Z. 2012, Phys. Plasmas, 19, 082901 [NASA ADS] [CrossRef] [Google Scholar]
  42. Matsumoto, Y., Amano, T., Kato, T. N., & Hoshino, M. 2017, Phys (Lett: Rev) [Google Scholar]
  43. Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
  44. Moskalenko, I. V., & Strong, A. W. 1998, ApJ, 493, 694 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  45. Moskalenko, I. V., Strong, A. W., Ormes, J. F., & Potgieter, M. S. 2002, ApJ, 565, 280 [NASA ADS] [CrossRef] [Google Scholar]
  46. Ohira, Y., Murase, K., & Yamazaki, R. 2010, A&A, 513, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Orlando, E. 2018, MNRAS, 475, 2724 [NASA ADS] [CrossRef] [Google Scholar]
  48. Petruk, O., & Kopytko, B. 2016, MNRAS, 462, 3104 [NASA ADS] [CrossRef] [Google Scholar]
  49. Petruk, O., Kuzyo, T., & Beshley, V. 2016, MNRAS, 456, 2343 [NASA ADS] [CrossRef] [Google Scholar]
  50. Pohl, M. 1993, A&A, 270, 91 [NASA ADS] [CrossRef] [Google Scholar]
  51. Ptuskin, V. S., & Zirakashvili, V. N. 2003, A&A, 403, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Ptuskin, V. S., & Zirakashvili, V. N. 2005, A&A, 429, 755 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
  53. Putze, A., Derome, L., Maurin, D., Perotto, L., & Taillet, R. 2009, A&A, 497, 991 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Reville, B., & Bell, A. R. 2013, MNRAS, 430, 2873 [NASA ADS] [CrossRef] [Google Scholar]
  55. Simpson, C. M., Pakmor, R., Marinacci, F., et al. 2016, ApJ, 827, L29 [NASA ADS] [CrossRef] [Google Scholar]
  56. Skilling, J. 1975, MNRAS, 172, 557 [NASA ADS] [CrossRef] [Google Scholar]
  57. Sushch, I., Brose, R., & Pohl, M. 2018, A&A, 618, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [NASA ADS] [CrossRef] [Google Scholar]
  59. Tang, X., & Chevalier, R. A. 2015, ApJ, 800, 103 [NASA ADS] [CrossRef] [Google Scholar]
  60. Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2012a, Astropart. Phys., 35, 300 [NASA ADS] [CrossRef] [Google Scholar]
  61. Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2012b, A&A, 541, A153 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Telezhinsky, I., Dwarkadas, V. V., & Pohl, M. 2013, A&A, 552, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Trotta, R., Jóhannesson, G., Moskalenko, I. V., et al. 2011, ApJ, 729, 106 [NASA ADS] [CrossRef] [Google Scholar]
  64. van Marle, A. J., Casse, F., & Marcowith, A. 2018, MNRAS, 473, 3394 [NASA ADS] [CrossRef] [Google Scholar]
  65. Völk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2003, A&A, 409, 563 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Yang, R.-Z., de Oña Wilhelmi, E., & Aharonian, F. 2014, A&A, 566, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Zeng, H., Xin, Y., & Liu, S. 2019, ApJ, 874, 50 [NASA ADS] [CrossRef] [Google Scholar]
  68. Zirakashvili, V. N., & Ptuskin, V. S. 2012, Astropart. Phys., 39, 12 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Core-collapse SNRs

We also modeled the acceleration of CRs at a core-collapse SN expanding in a wind zone, using the initial conditions described in Sushch et al. (2018) and a magnetic field that is ∝1/r. The parameters are given in Table A.1. They are roughly consistent with those expected for the environment of red-supergiants (RSG). After ≈1800 yr the density and magnetic-field strength are comparable to those used for the medium-density type-Ia simulations.

Table A.1.

Parameters of the SNR models.

The remnant size for an RSG is typically smaller than that for a type-Ia explosion, as the material around the star is initially much denser. However, the shock speed is comparable after 1800 years, and we calculated the gamma-ray emission for that time. The spectra are presented in Fig. A.1.

thumbnail Fig. A.1.

Top panel: gamma-ray spectra of IC (solid) and PD (dashed) emission for a RSG (red) and a type-Ia (black) SNR at 1800 years. Bottom panel: gamma-ray flux ratio of an RSG vs. a type-Ia environment for PD (dashed) and IC (solid) radiation.

It can be seen, that the PD gamma-ray flux is about three times higher in the CC case compared to the type-Ia scenario. Although the volume of a type-Ia SNR is roughly three times larger than that of its core-collapse sibling, Eq. (12) indicates that both the number of cosmic rays and the total amount of target material are increased by a factor three, hence the overall gain by a factor three. For similar reasons, the IC gamma-ray intensities are comparable between CC and type-Ia SNRs, as the differences caused by the size and the number of cosmic rays cancel out. The initially higher magnetic field in the RSG-scenario lowers the maximum energy of electrons, because synchrotron cooling is more efficient (see also Eq. (10)). As a result, the leptonic gamma-ray luminosity in the 1 − 10 TeV band is about six times lower in the RSG case than that in the type-Ia scenario. This indicates that the spatial structure around the SNR is one of the factors that determine the IC gamma-ray luminosity. All the above applies only as long as the remnant expands in the unshocked wind. As soon as the shock enters the shocked wind, it is again propagating in an area of approximately constant density. This transition usually happens earlier – after around 500 years for a RSG – than the age of 1800 years we are considering here. After entering the shocked wind, the impact of the initial wind-zone will fade, and our solutions for type-Ia will become valid again, because core-collapse and type-Ia explosions have a comparable explosion energy, which is – besides the ambient density – the crucial parameter that defines the evolution during the Sedov stage. A detailed modeling of SNRs expanding in the various environments around massive stars is important to understanding the gamma-ray emission of early SNRs, but that is clearly beyond the scope of this paper.

The role of background CRs is still negligible in the RSG scenario, even after the remnant entered the low-density environment of the shocked wind. In this case, the contribution of background CRs can be of the order of 1% when effects of stellar modulation are neglected. In reality, low-energy CRs will have been pushed out of the circumstellar medium by the wind of the progenitor star, and the background level of CRs will be even lower.

Appendix B: Alternative CR background spectrum

The background CR spectrum given in Eq. (2) is rather hard below 1 GeV, and it has been argued that the true spectrum might be softer (Orlando 2018). There has to be some low-energy cutoff, otherwise the energy deposited in low-energetic cosmic rays would not be finite, but the position of this cutoff is unknown. We investigated to what extent our spectra are affected by our choice of background CR spectrum. We repeated our reacceleration runs with the alternative background spectrum

N PR ( E ) = { N 0 E 2.75 for E 1 GeV N 0 E 1.5 for E < 1 GeV . $$ \begin{aligned} N_\mathrm{PR} (E) = {\left\{ \begin{array}{ll} N_0 E^{-2.75}\,\, \mathrm{for}\, E\ge \mathrm{1\,GeV} \\ N_0 E^{-1.5}\,\,\, \text{ for} E < \mathrm{1\,GeV} \end{array}\right.}\!\!\!\!\!\!. \end{aligned} $$(B.1)

As a result, more low-energetic CRs are present to drive the amplification of magnetic turbulence and the confinement of CRs around the SNR. The difference in the spectra is seen most significantly in the escape spectra. Figure B.1 compares the volume-averaged spectra in a shell 7.5 pc ahead of the forward shock after 37 kyr.

thumbnail Fig. B.1.

CR spectra about 7.5 pc ahead of the forward shock at an age of 37 kyr. The black line refers to the acceleration of particles from the thermal pool, where the red (blue) line refers to the reacceleration of the hard (soft) background spectrum below 1 GeV. It should be noted that the y-axis is scaled with p2.75.

The higher number of low-energy particles enables a stronger driving of turbulence and a better confinement of these particles around the SNR. As a result, the escape spectrum is softer for the softer CR background spectrum and approaches the spectral index of the CRs accelerated from the thermal pool of particles. In this sense, our original background spectrum represents a lower limit to number of low-energy particles that are available to drive the magnetic turbulence and confine the CRs in and around the SNR. The softer the CR spectrum below ≈1 GeV, the closer the resemblance of the reacceleration spectra to those arising from acceleration from the thermal pool.

All Tables

Table A.1.

Parameters of the SNR models.

All Figures

thumbnail Fig. 1.

Curve shows the radial position of the forward shock as a function of time. The transitions between free expansion and the Sedov phase, and between Sedov and the post-adiabatic phase, appear roughly at 1.3 kyr and 35 kyr, respectively.

In the text
thumbnail Fig. 2.

Volume-averaged downstream proton spectra for Bohm-like (red) and self-consistent (green) diffusion. The background of Galactic CRs is neglected.

In the text
thumbnail Fig. 3.

Total volume-averaged proton spectrum in the case of Bohm-like (red) and self-consistent (green) diffusion. Galactic CRs are ignored.

In the text
thumbnail Fig. 4.

Proton number-spectra for the Alfvénic diffusion scenario at different times. Solid lines represent the total, dashed the downstream and dotted the upstream spectra.

In the text
thumbnail Fig. 5.

Running power-law index of the electron (black) and proton (red) spectrum as function of energy displayed at different ages of the remnant. The Bohm-like diffusion scenario is shown in the top panel, and the Alfvénic-turbulence scenario in the bottom panel.

In the text
thumbnail Fig. 6.

Volume-averaged downstream proton spectra for Bohm-like (red) and self-consistent (green) diffusion at different times. We include background CRs and artificially decrease, with time, the injection fraction of thermal particles at the shock. The CR-density at 300 yr for nondecreasing injection (direct) is shown as a dotted gray line for comparison.

In the text
thumbnail Fig. 7.

Ratio of the radio flux produced by directly accelerated and reaccelerated electrons, normalized above 75 GHz. The assumed magnetic-field strength in the downstream region is 16 μG.

In the text
thumbnail Fig. 8.

Evolution of CR pressure over time for the two diffusion scenarios and with/without the Galactic CR background. The gray line indicates 10% of the shock ram pressure.

In the text
thumbnail Fig. 9.

CR spectra about 7.5 pc ahead of the forward shock at an age of 37 kyr. Red (black) lines refer to Bohm-like diffusion without (with) background CRs, and green (blue) lines stand for Alfvénic diffusion without (with) background CRs. Note that the y-axis is scaled with p2.75.

In the text
thumbnail Fig. 10.

Time evolution of the gamma-ray luminosity in the TeV band. The curves reflect models with Bohm-like (red) and Alfvénic (green) diffusion, and they indicate IC (solid lines) and PD (dashed lines) emission. The curves are calculated for an ambient gas density of 0.4 cm−3 assuming the same injection-efficiency for electrons and protons.

In the text
thumbnail Fig. 11.

Observed gamma-ray luminosity of galactic SNRs in the TeV band; blue markers indicate upper limits. The red curves reflect the models with Bohm-like diffusion and indicate inverse-Compton (solid lines) and pion-decay (dashed lines) emission. The thick curves are calculated for a ambient gas density of 0.4 cm−3. The shaded uncertainty bands represent the range of luminosity expected for a density in the range 0.04 − 4.0 cm−3. The references for the data points are given in the text.

In the text
thumbnail Fig. A.1.

Top panel: gamma-ray spectra of IC (solid) and PD (dashed) emission for a RSG (red) and a type-Ia (black) SNR at 1800 years. Bottom panel: gamma-ray flux ratio of an RSG vs. a type-Ia environment for PD (dashed) and IC (solid) radiation.

In the text
thumbnail Fig. B.1.

CR spectra about 7.5 pc ahead of the forward shock at an age of 37 kyr. The black line refers to the acceleration of particles from the thermal pool, where the red (blue) line refers to the reacceleration of the hard (soft) background spectrum below 1 GeV. It should be noted that the y-axis is scaled with p2.75.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.