© ESO 2021

1 Introduction

The detection of extended gamma-ray emission around two nearby pulsars – Geminga and PSR B0656+14 (Abeysekara et al. 2017a,b) – spawned an extensive discussion of TeV-halos around pulsars (Giacinti et al. 2020). Escaping high-energy electrons and positrons can be confined close to the accelerating pulsar-wind nebulae (PWN) by self-amplified magnetic turbulence and enhance the TeV gamma-ray emission in the vicinity of the PWN (Evoli et al. 2018). The resulting halos have an energy-dependent morphology, with a smaller extension toward the highest energies (Principe et al. 2020).

However, the interpretation of the observational data is not trivial. Simple uniform-diffusion models can explain the observed morphology of the gamma-ray emission on account of confining all electrons and positrons close to the pulsar. A contribution to the positron excess observed in the AMS data (Chang et al. 2008; Adriani et al. 2009) would thus be ruled out. In reality, a spatially nonuniform diffusion coefficient has to be expected for self-amplified turbulence and could both explain the observed morphology and allow for sufficient escape flux to support the local positron-flux (Profumo et al. 2018).

A similar self-regulation of the diffusion coefficient by escaping particles is known to exist around supernova remnants (SNRs) as well. There, the scattering turbulence is created mainly by escaping hadrons, and the electrons accelerated at the SNR blastwave are trapped as a side effect (Ohira et al. 2010; Nava et al. 2016).

So far, there has not been a direct measurement of cosmic rays (CRs) escaping from a SNR. A possible scenario would be CRs illuminating molecular clouds close to the SNRs where the enhanced target density boosts hadronic gamma-ray emission (Yan et al. 2012). There is evidence of larger extension of the gamma-ray emission around RXJ1713.7-3946 compared to the X-ray emission that could indicate CR escape (H.E.S.S. Collaboration 2018a). However, there is indirect evidence that CR escape has to happen around SNRs. Recent studies of particle acceleration in SNRs showed that typically soft, broken power-law spectra of aged SNRs (Zeng et al. 2019) can be produced by the escape of particles at the highest energy from the interior of the SNR. Here, the most energetic particles escape once the SNR is not capable of accelerating them any further, creating soft spectra inside the SNR where most of the emission is produced (Brose et al. 2020; Celli et al. 2019).

These concepts potentially resolve the tension between the s = 2 spectra predicted by shock-acceleration theory and the typically soft emission spectra of s ≥ 2.7 observed inside evolved SNRs, especially in the gamma-ray domain. However, the question of whether SNRs are the sources of Galactic CRs has yet to be answered. The earlier studies show that there is no contradiction between the observed soft emission spectra and the somewhat harder total-production spectra with s ≈ 2.2−2.4 that are predicted for the sources of Galactic CRs by Galactic propagation models (Trotta et al. 2011). The accumulated proton spectrum is approximately the same as the time integral over the SNRs lifetime and represents the CR yield that finally gets released into the sea of Galactic CRs once the SNR shock faded.The aim of this paper is to explore the observational signatures that can be expected from escaping CRs around SNRs.

2 Basic equations and assumptions

The methodology is similar to that of Brose et al. (2020), and here we shall only give a short summary. We combine 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 simulation of the hydrodynamical flow profiles.

2.1 Cosmic rays

We solve the kinetic equation for the differential number density of CRs, N, $$\begin{array}{ll}\frac{\partial N}{\partial t}=\hfill & \nabla (D_r\nabla N-uN)\hfill \\ \hfill & -\frac{\partial }{\partial p}\left((N\dot{p})-\frac{\nabla \cdot u}{3}Np\right)+Q\hfill \end{array}$$(1)

in the testparticle limit, where D_r denotes the spatial diffusion coefficient, u the advective velocity, ṗ energy losses (see Sect. 2.1.2), and Q the source of thermal particles (Skilling 1975).

We rewrite Eq. (1) by transforming the radial coordinate, r, to a new spatial coordinate that is co-moving with the SNR shock and provides a very high numerical resolution at the shock front (Brose et al. 2020, and references therin). The sea of Galactic CRs is neglected in this simulations. As long as the fraction of injected particles is constant over time, the emission from freshly injected particles is always dominant over the contribution from background CRs.

2.1.1 Injection

We use a thermal leakage model (Blasi et al. 2005; Malkov 1998) for the injection of particles. Here, the efficiency of injection η_i is given by $${\eta }_i=\frac{4}{3\sqrt{\pi }}(\sigma -1){\psi }^3{\text{e}}^{-{\psi }^2},$$(2)

where σ is the shock compression ratio, and ψ is the multiple of the thermal momentum, at which we inject particles.

Several authors noted that the bipolar morphology in the nonthermal emission of SN 1006 can be attributed to effects of the shock-obliquity on the injection or acceleration efficiency (Völk et al. 2003; Petruk et al. 2009; Beshley & Petruk 2012; Pais & Pfrommer 2020). This notion seem to be supported by the hybrid-simulations of Caprioli & Spitkovsky (2014), where the acceleration at quasi-perpendicular shocks is strongly suppressed. Reville & Bell (2013) on the other hand found a quasi-universal behavior of shocks irrespective of the magnetic field orientation very far upstream of the shocks by using a spherical-harmonics expansion of the CR Fokker–Planck equation. This suggests that injection may only be weakly dependent on the shock orientation. We have to ignore possible obliquity effects in our spherically symmetric model and assume a quasi-parallel configuration across the shock-surface.

This injection scenario is a simplification, in particular 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.

The multiple of the thermal momentum, ψ, determines the fraction of thermal particles that get turned into CRs. However, ψ and hence η are only weakly constrained. Low values for ψ are required in scenarios featuring nonlinear modifications of the shock structure by CR pressure. The originally proposed lower boundary ψ ≈ 3.5 ensures that thermal particles have a mean-free path larger than the shock thickness and can participate in the DSA-process. For this work we choose ψ = 4.2, which guarantees a CR pressure of less than 2.5% of the shock ram pressure during the entire simulation. Likewise, this value is close to the amount of injected particles seen in SN 1006, where the likely leptonic origin of the emission allows a reasonably sound estimate of η_i (see Appendix A for details). We have chosen an equal amount of injected electrons and protons.

2.1.2 Inverse-Compton losses

In addition to the synchrotron losses that were considered already in earlier versions of RATPaC, electrons will also suffer losses by collisions with photons from background photon-fields. Usually these inverse-Compton (IC) losses can be neglected as synchrotron losses in amplified magnetic field will dominate. However, high-energy electrons spend a large fraction of time upstream of the shock, where IC and synchrotron losses are of similar strength.

We used the approximations derived by Reimer et al. (2006) to account for IC losses by collisions with the cosmic microwave background.We note that additional photon fields might have to be taken into account when modeling core-collapse supernova remnants where strong local infrared and optical photon fields have to be expected or when specific Type-Ia SNRs are modeled for which estimates of the background photon-fields are available based on their location in the Galactic disk.

2.2 Magnetic turbulence

In parallel to the transport equation for CRs, we solved a transport equation for the magnetic-turbulence spectrum, assuming Alfvén waves only. The temporal and spatial evolution of the spectral energy-density per unit logarithmic bandwidth, E_w, is described by $$\begin{array}{l}\frac{\partial E_w}{\partial t}+\cdot \nabla (u{E}_w)+k\frac{\partial }{\partial k}\left(k^2{D}_k\frac{\partial }{\partial k}\frac{E_w}{k^3}\right)\\ =2({\Gamma }_g-{\Gamma }_d)E_w.\end{array}$$(3)

Here, u denotes the advection velocity, k the wavenumber, D_k the diffusion coefficient in wavenumber space, and Γ_g and Γ_d the growth and damping terms, respectively (Brose et al. 2016).

We calculated the diffusion coefficient from E_w using $$\begin{array}{ll}{D}_r\hfill & =\frac{4v}{3\pi }{r}_g\frac{U_m}{E_w},\hfill \end{array}$$(4)

where U_m denotes the energy density of the large-scale magnetic field, v is the particle velocity, and r_g the gyro-radius of the particle.

As initial condition, we used a diffusion coefficient, and hence a turbulence spectrum, as suggested by Galactic propagation modeling (Trotta et al. 2011), but reduced by a factor ten on account of numerical constraints, $$\begin{array}{ll}{D}_0\hfill & ={10}^{28}{\left(\frac{pc}{10\text{GeV}}\right)}^{1/3}{\left(\frac{B_0}{3\text{μG}}\right)}^{-1/3}\text{. }\hfill \end{array}$$(5)

This choice is a factor of 10 higher than in our previous works.

We use growth-rate based on the resonant streaming instability (Skilling 1975; Bell 1978), $$\begin{array}{ll}{\Gamma }_g\hfill & =A\cdot \frac{v_{\text{A}}{p}^{2}v}{3E_{\text{w}}}\left|\frac{\partial N}{\partial r}\right|,\hfill \end{array}$$(6)

where v_A is the Alfven-velocity. We introduced a linear scaling factor, A, to artificially enhance the amplification. We used A = 10 throughout this paper to mimic the more efficient amplification due to the nonresonant streaming instability (Lucek & Bell 2000; Bell 2004). The particulars of the nonresonant amplification are beyond the present capabilities of our code since the back-reaction of CR streaming that terminates the wave-growth, a modification of the bulk flow (Riquelme & Spitkovsky 2009; Niemiec et al. 2010; Kobzar et al. 2017), cannot be accounted for. Likewise difficult to handle, and in fact quite unclear, is the scattering efficiency of the nonresonant modes. The value we choose for A guarantees a cut-off energy in the gamma-ray spectrum between 1 and 10 TeV as observed in young SNRs.

We calculate the total magnetic field strength as $$\begin{array}{ll}{B}_{\text{tot}}\hfill & =\sqrt{B_0^2+B_{\text{Turb}}^2},\hfill \end{array}$$(7)

where B₀ is the large-scale magnetic field. We solve the induction equation to model the transport of the frozen-in large-scale magnetic field (Telezhinsky et al. 2013). The far-upstream field is assumed to be uniform with strength 5 μG and the field is assumed to be fully turbulent, resulting in a magnetic field-compression of $B_{\text{d}}=\sqrt{11}{B}_{\text{u}}$.

Since we exceed the growth rate of the resonant streaming instability (Bell 1978) by a factor of ten, the turbulent field is amplified to δB ≫ B₀ during the initial phases of SNR evolution. The peak amplitude of the field is reached right at the shock. Downstream the field strength quickly falls due to efficient cascading of turbulence. The resulting magnetic-field profiles resemble the profiles suggested by Pohl et al. (2005). In the early phases of SNR evolution peak amplitudes of ≈ 100 μG are reached, whereas after 1000 yr one finds ≈40 μG – a value compatible with the 30 μG estimated for the 1000 yr-old SNR SN1006 (Acero et al. 2010).

The growth of the magnetic turbulence and hence the magnetic field is balanced by cascading. This process is described as a diffusion process in wavenumber space, and the diffusion-coefficient is given by (Zhou & Matthaeus 1990; Schlickeiser 2002) $$\begin{array}{ll}{D}_{\text{k}}\hfill & =k^3{v}_{\text{A}}\sqrt{\frac{E_{\text{w}}}{2B_0^2}}.\hfill \end{array}$$(8)

This phenomenological treatment will result in a Kolmogorov-like spectrum, if cascading is dominant. Since v_A ∝ B_tot, the cascading rate will depend on the level of magnetic turbulence in two different regimes $$D_k\propto \{\begin{array}{ll}\sqrt{E_{\text{w}}}\hfill & \text{for}{\text{E}}_{\text{w}}<{\text{B}}_{\text{0}}^{\text{2}}\text{/8}\pi \hfill \\ E_{\text{w}}^{3/2}\hfill & \text{for}{\text{E}}_{\text{w}}>{\text{B}}_{\text{0}}^{\text{2}}\text{/8}\pi .\hfill \end{array}$$(9)

Once the turbulent field dominates over the background field, the cascading rate depends more sensitively on the energy density of magnetic turbulence. Without requiring other damping mechanisms, the enhanced cascading efficiently limits the maximum level of turbulence to a level commensurate with that derived from SNR observations and PIC simulations (Vink 2006; Riquelme & Spitkovsky 2009; Niemiec et al. 2010).

2.3 Thermal plasma

In the test-particle limit, the evolution of an SNR can be described with the standard gas-dynamical equations: $$\begin{array}{ll}\frac{\partial }{\partial t}\left(\begin{array}{c}\rho \\ m\\ E\end{array}\right)+\nabla {\left(\begin{array}{c}\rho v\\ mv+PI\\ (E+p)v\end{array}\right)}^T\hfill & =\left(\begin{array}{c}0\\ 0\\ L\end{array}\right)\hfill \\ \frac{\rho v^2}{2}+\frac{P}{\gamma -1}\hfill & =E\text{,}\hfill \end{array}$$

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. We solve this system of equations 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 the results for Type-Ia supernova explosions. We initiate the simulations with exponential-ejecta profiles: $$\begin{array}{l}{\rho }_{\text{SN}}=A\exp (-v/v_e)t_i^{-3,}\text{ and }v=r/t_i,\\ v_e={\left(\frac{E_{\text{ex}}}{6M_{\text{ej}}}\right)}^{1/2,}\text{ and }A=\frac{6^{3/2}}{8\pi }\frac{M_{\text{ej}}^{5/2}}{E_{\text{ex}}^{3/2}}\end{array}$$

as initial conditions (Dwarkadas & Chevalier 1998). Here, t_i = 2.5 yrs is the start time of our simulation, M_ej = 1.4M_sol the ejecta mass, E_ex = 10⁵¹ erg the explosion energy, and r the spatial coordinate. The density of the ambient medium was chosen to be 0.4 cm⁻³.

3 Results

We followed the evolution of the remnant for 25 000 yr. The remnant enters the Sedov–Taylor phase after 1300 yr¹ and would enter the post-adiabatic phase after 35 000 yr.The total magnetic field reaches 90 μG after 300 yr and drops to 40 μG after 1000 yr. After 10 000 yr, the turbulent field is weaker than the compressed large-scale field. In the following, we first describe the morphology of the TeV-halo over the lifetime of the remnant and then examine its detectability with today’s and future gamma-ray experiments.

Fig. 1 Maps of normalized surface brightness of a Type-Ia SNR after 300, 1000, 2000 and 10 000 yr (from top left to bottom right). Left hemispheres are for IC emission and right hemispheres for PD emission.

3.1 Halo evolution

We calculated intensity maps for inverse-Compton and Pion-decay (PD) emission at three energies over the lifetime of the SNR. The results are presented in Fig. 1.

Initially, both radiation mechanisms produce a shell-like morphology in all energy bands. From roughly 1000 yr on, the shell thickness of IC emission exceeds that of PD radiation, because the latter is boosted by the high gas density immediately downstreamof the shock. Contrary, the IC emission reflects only the distribution of electrons, and already after 1000 yr we notice IC emission outside the SNRs shell. The age of 1000 yr roughly corresponds to the time when the SNR luminosity peaks regardless of the gamma-ray emission mechanism, on account of the transition to the Sedov–Taylor stage. For an ambient density of 0.4 cm⁻³ and an average free-expansion velocity of 5000 km s⁻¹, five solar masses of material will have been swept-up by the shock within the first 1000 yr of expansion². The transition to the Sedov-phase usually also marks the time of the highest maximum energy of particles (Ptuskin & Zirakashvili 2003). At later times, particles of the highest energy start to escape from the remnant (Brose et al. 2020).

Even after 10 000 yr, the PD morphology remains shell-like whereas that of IC emission becomes center-filled. Even low-energy CRs propagated into the center of the remnant, and the projection enhances the brightness towards the center of the remnant in the IC channel. The low density of the thermal plasma in the center keeps suppressing PD emission from the central region. As more and more high-energy electrons escape with increasing age, an extensive halo of > 100-GeV electrons is formed around the remnant. There is also a faint halo in the PD channel after 10 000 yr, but it is much weaker than the IC halo on account of the low gas density.

Synchrotron losses most strongly modify the distribution of high-energy electrons. They cause the thinner IC shell at 1 TeV compared to lower energies that is visible in the intensity maps. The magnetic field peaks at the shock, and so synchrotron losses are strongest there. The same effect is responsible for the relative smaller extend of the bright IC-region at 2000 yr. In the unshocked ejecta the field strength is 0.2 μG, making IC scattering the dominant energy-loss mechanism. Generally, the IC-loss timescale exceeds that for synchrotron losses where the magnetic field surpasses 5 μG.

Figure 2 shows the projected profiles of gamma-ray intensity at four stages of the remnant evolution. It is clearly visible that the extension of the IC halo increases with energy, in contrast to the case of PWNs. The size of the halo is not defined by a balance between acceleration of the highest energetic electrons and their synchrotron cooling, as all high-energetic particles have been accelerated at early times. More important is that the high-energy particles experience the largest diffusion coefficient and thus can fill a larger halo.

The morphology of the remnant also strongly depends on irregularities in the ambient medium. Deviations from spherical symmetry may be caused if the coherence length of the turbulent ambient magnetic field is larger than the size of the remnant (Pais & Pfrommer 2020). Similar effects have been described earlier for the hadronic and leptonic emission morphology of remnants expanding in a uniform ambient field (Petruk et al. 2009; Beshley & Petruk 2012). Our results qualitatively agree with these conclusions in regions where the field is parallel to the shock or where the coherence-scale of the ambient field is smaller than the size of the remnant. However, accelerated diffusion into the center of the remnant with the onset of the Sedov-stage produces a more center-filled morphology for the IC-emission than can be obtained by models relying solely on advection.

A comparison of our predictions with measurements is difficult as there is no firm detection of an emission halo around a SNR (see also Sect. 3.3 for details). However, extended gamma-ray emission around the two Pulsars Geminga and PSR B0656+14 has been detected (Abeysekara et al. 2017a). The measurements show a roughly exponential decrease of the surface-brightness with increasing radial distance from the central pulsar. This behavior is also reproduced in simple models using one or two zones with spatially constant diffusion coefficients around the pulsars (Di Mauro et al. 2020). We obtain a similar behavior in our halo-profiles (see Fig. 2) for those evolutionary stages during which CRs are still accelerated to the highest energies. However, near the cut-off energy the profiles transition to a more linear trend reflecting that those particles simply escape. This behavior is more pronounced for electrons, on account of the synchrotron cooling of the high-energy electrons in the strong magnetic field inside the remnant.

Fig. 2 Projected intensity profiles normalized to the surface brightness at the shock at four stages of SNR evolution.

Fig. 3 Normalized diffusion coefficient in the upstream region as function of distance to the shock and the SNR age. The left panel is for a particle energy of 100 GeV, the right panel for 3 TeV. The dashed lines correspond to 1%, 10% and 50% of DISM.

3.2 Reduction of the diffusion coefficient

The escape of particles from the SNR and consequently the amplification of turbulence change the diffusion coefficient in the vicinity of the SNR (Fujita et al. 2010, 2011). The escape of high-energy particles from their acceleration sites gained new attention with the detection of an extended halo around the Pulsars Geminga and PSR B0656+14 (Abeysekara et al. 2017a).

Figure 3 shows the spatial variation of the diffusion coefficient relative to the assumed diffusion coefficient in the ISM, which is set to 10% of the conventional Galactic diffusion coefficient (Trotta et al. 2011). This choice was made to keep the timestep in our simulations reasonably large, otherwise the simulation of very old SNRs would be impossible.

For both electrons and ions the reduction of the diffusion coefficient is strongest after 5000 yr (10%-level) and 2000 yr (1%-level), around the time of, or shortly after, the transition from free expansion to the Sedov–Taylor phase, when the maximum energy of the accelerated particles starts to decrease. Then there are too few freshly accelerated particles at high energy to sustain the level of turbulence, and consequently the diffusion coefficient starts to increase. The main damping mechanism of turbulence is cascading which, in contrast to other studies in the context of Pulsars (Evoli et al. 2018), we treat not as a simple loss term but account for the energy transfer to smaller scales. Consequently, the distributions in Fig. 3 show an energy-dependent time evolution. The spatial extent is similar though, because the suppression of the diffusion coefficient for 100 GeV-particles is governed mainly by down-cascading from larger scales, that is by turbulence driven by particles at higher energy. Turbulence driving by the low-energy particles is only important very close to the shock.

To better understand the energy-dependent time evolution of upstream diffusion after the end of the free-expansion phase, it is instructive to closely inspect the cascading time (Schlickeiser 2002), $$\begin{array}{ll}\tau (k)\hfill & =\frac{1}{v_{A}k}\sqrt{\frac{2U_B}{E_w}},\hfill \end{array}$$(14)

which depends on the wavenumber, k, the spectral-energy density of the turbulent field, E_w, the Alfvén speed, v_A, and the energy density in the large-scale magnetic field, U_B. At all wavenumbers, except at k_min resonant with the most energetic particles, turbulence energy is cascaded from large scales (small k) to small scales (large k). So, in a quasi-steady state, the level of turbulence at small scales is sustained by a continuous influx from larger scales.

At an age of 4000 yr, the cascading times for turbulence resonant with 3 TeV (100 GeV) particles 5 pc ahead of the shock is about 500 yr (30 yr). It takes 4300 yr for the shock toarrive in this region, suggesting that cascading is efficient. However, the full spectral transport of turbulence is in many cases slower than the cascading time suggests, implying that the variation in the level of turbulence arises from a changing CR density gradient. As the freshly accelerated particles become fewer and less energetic with time, particles escaping from deep downstream become more important for turbulence driving.

During the initial growth of turbulence in the precursor, for the two particle energies shown in Fig. 3 the diffusion coefficient shows the same trend, until the supply of freshly accelerated 3 TeV-particles is exhausted. However, the precursor scale of 100 GeV-particles is significantly smaller than for 3 TeV-particles, because the turbulence scattering the lower-energy particles results from cascading of modes resonant with the more energetic particles. Later, the escape of particles, that were trapped in the interior of the remnant, becomes the main driver of turbulence, and in fact the maximum extent of the turbulence precursor, and hence the reduced diffusion coefficient, is reached well after the beginning of the Sedov-phase.

It is important to emphasize again that the diffusion coefficient for low-energy particles is determined by the turbulence that the high-energy particles provide. This has consequences beyond SNRs, as for example low-energy particles in the halos of PWNs will reside closer to their acceleration sites than is suggested, if one treats cascading as a simple damping mechanism (Evoli et al. 2018).

3.3 Detectability

Figure 4 shows the emission spectra expected from the SNR itself and the halo at different times. We assumed a generic distance of 1 kpc to the remnant. The contribution from the halo is shown twice, once as volume-integrated emission from r > r_sh, and once only the component that in projection appears to come from beyond the projected radius of the SNR. We modeled 3D objects, but we observed them in their 2D appearance in the sky. Thus, part of the emission from the halo will in projection appear to come from the SNR-interior, as the point of origin is in front of or behind the remnant.

The evolution of the spectral energy distributions (SEDs) for PD emission clearly shows the mechanism for the spectral softening: as highly energetic protons escape from the inside of the remnant, the low-energy particles remain inside, leading to a soft spectrum (Brose et al. 2020). The high-energy particles outside the SNR lack target material for the production of significant emission, and hence hadronic halo emission is unlikely detectable with current-generation gamma-ray observatories, even if the emission from the remnant and the halo could be disentangled.

In case of the IC-emission, cooling becomes relevant after about 1000 yr, and the spectra from the interior of the remnant are modified by it. The halo emission from IC scattering is always much brighter than the PD component. Overall, the flux from the halo is only 20% to 30% of that of the SNR itself, making a detection possible only for the brightest known Galactic SNRs. As low-energy particles reside closer to the shock, they contribute little to the projected halo emission whose spectrum is consequently harder than the overall emission from the remnant.

Figure 4 demonstrates that leptonic emission detected from an SNR by current-generation gamma-ray observatories probably includes emission from the halo. Most of the halo emission is produced close to the remnant, and in projection it is very difficult to observationally separate it from emission from the interior. It will probably require the tenfold higher sensitivity of the future Cerenkov Telescope Array (CTA) to directly detect halo emission from the brightest TeV-remnants. What may have been seen already in RXJ1713.7-3946 is that the leptonic halo causes larger extent of the very-high gamma-ray emission compared to that of the X-ray emission.

A special case is the remnant SN 1006 that has a bipolar morphology of nonthermal emission (Pye et al. 1981; Acero et al. 2010), possibly originating from the variation of shock obliquity in a dominant large-scale magnetic field (Völk et al. 2003). Regions with a quasi-perpendicular magnetic field are associated with inefficient ion acceleration (Caprioli & Spitkovsky 2014). The apparent alignment of the bipolar morphology with the plane of the sky means that the emission from SN 1006 is much less distorted by projection effects than are more spherical remnants. If so, the contribution from the halo to the high-energy gamma-ray emission could be a factor of a few higher than that from other remnants of similar age. A direct measurement of the halo emission might still be difficult due to the limited spatial resolution of current-generation Imaging Air Cerenkov telescopes (IACTs) and the low flux from SN 1006. In any case, a large fraction of the most energetic electrons resides in the upstream region (or halo). Due to the lower magnetic-field strength in that region, they contribute relatively little to the X-ray emission of SN 1006 (see also Sect. 3.5) and a lot to the gamma-ray emission. Figure 4indicates that this could lead to a spectral discrepancy between X-rays and high-energy gamma rays that is difficult to reproduce in a one-zone model. In the past, such a finding for SN 1006 was interpreted as effect of nonlinear shock modification (Acero et al. 2010). We suggest that this discrepancy is simply a natural result of halo emission.

Fig. 4 Comparison of emission spectra for IC (red) and PD (black) emission. Emission from the SNR is filled-black and filled-red forPD and IC-emission respectively. The halo-emission is filled-orange and filled-gray for IC and PD-emission respectively at 300, 1000, 2000 and 10 000 yr (from top left to bottom right). The spectra including (and excluding) the projection effect constitute the upper (lower) boundaries of the SNR emission and the lower (upper) limit of the area for emission from the halo.

Fig. 5 Maps of spectral-index variation in three energy bands at 300, 1000, 2000 and 10 000 yr (from top left to bottom right), relative to the average spectral index of the remnant (including the halo). Left hemispheres are inverse-Compton emission and right hemispheres Pion-decay emission.

3.4 Spectral index distribution

The angular resolution of gamma-ray observatories is key to a study of spatial variations of the gamma-ray spectra (H.E.S.S. Collaboration 2018a,b; Humensky & VERITAS Collaboration 2015). So far, no significant variation across the remnant has been detected.

The spectrum of emission from outside the remnant is typically harder than that from the interior. Figure 5 illustrates spectral index variations across the remnant at different ages. For emission from the interior of the remnant, the deviations are moderate. Depending on the energy-band, |Δs| can reach values of up to 0.8, and in general the emission from the interior of the remnant is softer than the average spectral index from the remnant. However, |Δs|≲ 0.2 is usually fulfilled for all but the highest energies of IC-emission. There, the mean spectral index is close to that seen upstream of the shock due to the large contribution of upstream-electrons to the gamma-ray flux. Still, the spectral index is quite uniform across the interior of the remnant, whereas the deviations of the spectral index of the halo emission are typically larger.

The detection of hadronic halo emission is unlikely given the current observational sensitivities. In the H.E.S.S. data for RXJ1713.7-3946 the 1-σ statistical uncertainty in the spatially-resolved spectral index maps is Δs ≈ ± 0.1 (H.E.S.S. Collaboration 2018a), likewise Δs ≈ ± 0.15 for Vela Junior (H.E.S.S. Collaboration 2018b). The systematic and the statistical uncertainties are comparable. In the VERITAS results for IC443, the spectral uncertainty ranges between Δs = ±0.11 and Δs = ±0.42 (Humensky & VERITAS Collaboration 2015). Considering the limited spatial resolution and the poor statistics in some regions, a significant detection of spectral-index deviations is unlikely with current-generation IACTs.

The higher sensitivity of CTA should reduce the uncertainty, Δs, by roughly a factor of three, which should allow the detection of spectral-index variation with 3σ significance. We note that deviations from a spherical symmetry or asymmetric environments, such as density gradients in the ISM or fast motion of the progenitor (Meyer et al. 2020, 2021), could lead to additional significant variations in the measured spectral index across the remnant.

Fig. 6 Synchrotron flux from the remnant at various times. The upper boundaries of the filled areas mark the total emission from the remnant while the lower boundaries indicate emission from downstream only. The green band indicates the 1–10 GHz range.

3.5 Synchrotron emission

We evaluated the nonthermal synchrotron emission from the SNR over its lifetime. Figure 6 shows the time evolution of the total synchrotron flux, as well as the flux from the downstream region only. Contrary to the IC-emission, there is little difference between the total and the downstream-only flux, except for very old remnants, because a strong magnetic field boosts the emission from the downstream. At later times, most of the high-energy electrons were able to escape the remnant and the magnetic field is not amplified any more, the halo emission starts to dominate.

The flux in the X-ray band decreases with time, because the magnetic-field strength decreases and cooling takes its toll in the electron spectra. Consequently, the remnant is brightest and best detectable during the first 2000 yr after the supernova – matching the census of Galactic SNRs (Vink 2012).

In contrast, the radio flux increases throughout the simulations. Energy losses are not important, and the steady accumulation of GeV-scale particles compensates for the weakening magnetic field. Interestingly, the escape of electrons affects the spectral index in the radio band. Between 1 GHz and 30 GHz (Roughly 10⁻⁸ keV to 10⁻⁷ keV), the radio spectra are fairly soft after 10 000 yrs, with a spectral index α ≈ 0.75 (S_ν ∝ ν^−α).

Figure 7 shows the emission morphology at radio and X-ray energies for four stages of SNR evolution. Throughout the SNR’s lifetime the magnetic-field strength peaks right downstream of the forward shock. Initially, when the turbulent magnetic field exceeds the compressed large-scale field, cascading causes a decline toward the interior. Later, the trend persists due to the flow structure in the downstream of the shock. Consequently, the synchrotron emissivity peaks immediately downstream to the shock-location.

The X-ray morphology consists of a thin shell throughout the entire lifetime even though a part of the soft X-ray emission originates from the halo after 2000 yr and might be detectable for bright sources. The cutoff of the synchrotron spectrum is in or slightly below the X-ray band, and so the X-ray output is very sensitive to the structure of the magnetic field (Pohl et al. 2005). A slight shift in the characteristic synchrotron frequency due to a a spatial variation of the magnetic field can impose a strong variation in the synchrotron flux.

The radio emission shows the same shell-like structure. However, the spatial variation of the magnetic field has a moderate impact on the synchrotron emissivity, typically ∝ B^1.5. The radio shells appear thicker as a consequence. Moreover, diffusion of electron towards the center of the remnant partially compensates the weak magnetic field in the interior. Thus, the remnant appears somewhat center-filled in the radio band at later stages. Additionally, 4% of the radio emission at 10 000 yr originates from the halo. The magnetic field is dominated by the large-scale field at that time, which boosts the radio emission from a downstream electron by a factor of 6 compared to that of an upstream electron. This means that 20% of the low-energy electrons reside upstream of the shock and are escaping the remnant already.

To explain the center-filled radio morphology of the remnant, Sushch & Hnatyk (2014) argued that additional, inward-moving shocks are needed to almost uniformly fill the interior of the Vela SNR with radio-emitting electrons. Our findings suggest that the same can be achieved by diffusive transport of electrons, in particular when in later stages of evolution the diffusion coefficient inside and outside of the remnant increases. However, Vela is likely the result ofa core-collapse SN and thus additional factors need to be taken into account. Core-collapse SNRs initially evolve in the freely expanding wind of the progenitor star, in which the magnetic field will follow a B ∝ 1∕r dependence.Beyond the wind-termination shock, the magnetic-field strength is almost constant out to the contact discontinuity (Sushch et al. 2018, and in prep.), which strongly affects the X-ray and radio maps and likely enhances the trend towards a center-filled radio morphology.

Fig. 7 Maps of normalized synchrotron surface brightness of a Type-Ia SNR after 300, 1000, 2000 and 10 000 yr (from top left to bottom right).

4 Conclusions

We performed numerical simulations of particle acceleration in SNRs, solving time-dependent transport equations of CRs and magnetic turbulence in the test-particle limit alongside the standard gas-dynamical equations for Type-Ia SNRs. We derived the CR diffusion coefficient from the spectrum of magnetic turbulence, that evolves through driving by the resonant streaming instability as well as cascading and wave damping.

The gamma-ray morphology of SNRs depends on the radiation mechanism. Leptonic gamma-ray emission tends to produce a center-filled appearance with increasing remnant age, whereas hadronic gamma-ray emitters retain a shell-like morphology, on account of the distribution of target material for p–p interactions.

We showed that a gamma-ray halo has to be expected for leptonic emitters, despite the projection effects that place part of the halo emission inside the projected SNR radius. Whereas a firm detection of halo emission is in sight for the leptonic case, the halo component of hadronic emission likely remains undetectable, unless additional target material is present in the vicinity of the SNR.

Our simulations show that the most energetic particles tend to reside outside the SNR. The relatively weak magnetic field in the upstream region leads to a low intensity of X-ray emission from these high-energetic electrons, making the X-ray spectrum appear softer than expected for a uniform magnetic field. It is likely that the observed spectral between X-ray and gamma-ray emission from SN 1006 is based on this effect instead of nonlinear shock-modification.

Upstream of the forward shock, the diffusion coefficient is reduced, confining the accelerated particles close to the SNR. There is a strong spatial and temporal evolution of the diffusion coefficient in the entire precursor. The diffusion coefficient of low-energy particles is determined by cascading of turbulence downward from larger scales. In contrast to PWNs, for which neither the abundance of relativistic ions nor the time profile of particle acceleration are known, for SNRs we can self-consistently determine the injection spectra of CRs that will create the region of reduced diffusion around SNRs.

An investigation of the distribution of gamma-ray spectral indices across the interior of SNRs showed that no deviation beyond a |δs| ≈ 0.2 has to be expected for the known and resolvable Galactic SNRs. Detecting this small deviations is beyond the capabilities of current-generation IACTs but might be possible with CTA. Strong deviations from spherical symmetry can further enhance the detectability.

Initially the remnant emits nonthermal X-rays above 0.1 keV but cooling and a decreasing magnetic field reduce the X-ray luminosity and make a detection unlikely after a few thousand years. The radio flux from the remnant increases throughout its lifetime, and the radio morphology evolves from shell-like to a more center-filled appearance.

Acknowledgements

Robert Brose acknowledges funding from an Irish Research Council Starting Laureate Award (IRCLA/2017/83).

Appendix A Electron injection derived from SN1006

The injection fraction, η_i, is a crucialparameter for our simula- tions. The normalization of the proton spectrum determines the gradient of the CR distribution in the upstream region and hence the growth rate of magnetic turbulence. At late times, η_i can be understood as the number ratio of particles injected into shock acceleration and the total number of thermal particles in the system, $$\begin{array}{ll}{\eta }_i\hfill & =\frac{{\displaystyle {\int }_{E_i}^{\infty }N}(E)\text{d}E}{{\displaystyle \int {\rho }_0}/m\text{d}V}\text{ , }\hfill \end{array}$$(A.1)

where E_i is the injection energy, N(E) the differential number density of CRs, and ρ₀∕m the number density of the ambient medium. If η_i is known, equation (2) can be used to determine ψ for the setup of the simulation.

We note that η_i can also be determined from SNR observations. We pointed out that the luminosity-age distribution of Galactic SNRs can be related to the total number of acceler- ated electrons in these remnants, if one considers the density profile for the Sedov solution. Then, the gamma-ray emission above 1 TeV is dominated by the IC channel (Brose et al. 2020).

If the ambient density is very low, as it is for SN 1006, η_i can also be inferred directly, because it is likely that all of the high-energy photons are produced via IC emission. The analysis of the H.E.S.S.-collaboration revealed a total energy in electrons of W_e = 3.3 ⋅ 10⁴⁷ erg, an electron spectral-index of s_e = 2.1, and a electron cutoff-energy of E_c = 10 TeV (Acero et al. 2010). These values can now be used to determine the total number of electrons present in SN 1006. We assume for simplicity that energetic electrons are present with constant spectral index be- tween E₀ = 10 MeV and E_c = 10 TeV, so that $$\begin{array}{ll}{W}_e\hfill & ={\displaystyle {\int }_{E_0}^{E_c}E}\cdot N(E)\text{d}E\text{ , }\hfill \\ N(E)\hfill & =N_0\cdot E^{-s}\text{for}{\text{E}}_{\text{0}}<\text{E}<{\text{E}}_{\text{c}}\text{ , }\hfill \\ \frac{W_e}{{\displaystyle {\int }_{E_0}^{E_c}N}(E)\text{d}E}\hfill & =\frac{E_c^{2-s}-E_0^{2-s}}{E_c^{1-s}-E_0^{1-s}}\frac{1-s}{2-s}=\alpha \text{ , }\hfill \\ {\displaystyle {\int }_{E_0}^{E_c}N}(E)\text{d}E\hfill & =\frac{W_e}{\alpha }=6.4\cdot {10}^{51}\text{. }\hfill \end{array}$$

The total number of thermal particles can be calculated by integrating the density of the ambient medium over the volume filled by the SNR. SN 1006 has an reported radius of R = 9.6 pc in a medium with ρ₀ ∕m = 0.085 cm⁻³. $$\begin{array}{ll}{\displaystyle \int \frac{{\rho }_0}{m}}\text{d}V\hfill & =\frac{4}{3}\pi R^3\frac{{\rho }_0}{m}=9.26\cdot {10}^{57}\text{. }\hfill \end{array}$$(A.6)

This gives η_i = 6.9 ⋅ 10⁻⁷ for electrons. A correction is in order, because not the entire shock surface of SN 1006 is efficiently accelerating particles. Estimating that only  20 % of SN 1006’s surface are contributing to the high-energy emission, we find η_i = 3.5 × 10⁻⁶, which corresponds to ψ = 4.205. This is very close to the value ψ = 4.2 that we used in our simulations.

Experimental data suggests that for particles arriving at earth the electron-to-proton ratio is K_ep = 0.01. We inject 1.2 times more electrons than protons to account for a composition with 10 % helium. Since electrons and protons are injected at the same ψ but at different energies, this results in K_ep ≈ 0.01.

References

All Figures

	Fig. 1 Maps of normalized surface brightness of a Type-Ia SNR after 300, 1000, 2000 and 10 000 yr (from top left to bottom right). Left hemispheres are for IC emission and right hemispheres for PD emission.
In the text

	Fig. 2 Projected intensity profiles normalized to the surface brightness at the shock at four stages of SNR evolution.
In the text

	Fig. 3 Normalized diffusion coefficient in the upstream region as function of distance to the shock and the SNR age. The left panel is for a particle energy of 100 GeV, the right panel for 3 TeV. The dashed lines correspond to 1%, 10% and 50% of DISM.
In the text

Fig. 4 Comparison of emission spectra for IC (red) and PD (black) emission. Emission from the SNR is filled-black and filled-red forPD and IC-emission respectively. The halo-emission is filled-orange and filled-gray for IC and PD-emission respectively at 300, 1000, 2000 and 10 000 yr (from top left to bottom right). The spectra including (and excluding) the projection effect constitute the upper (lower) boundaries of the SNR emission and the lower (upper) limit of the area for emission from the halo.

In the text

	Fig. 5 Maps of spectral-index variation in three energy bands at 300, 1000, 2000 and 10 000 yr (from top left to bottom right), relative to the average spectral index of the remnant (including the halo). Left hemispheres are inverse-Compton emission and right hemispheres Pion-decay emission.
In the text

	Fig. 6 Synchrotron flux from the remnant at various times. The upper boundaries of the filled areas mark the total emission from the remnant while the lower boundaries indicate emission from downstream only. The green band indicates the 1–10 GHz range.
In the text

	Fig. 7 Maps of normalized synchrotron surface brightness of a Type-Ia SNR after 300, 1000, 2000 and 10 000 yr (from top left to bottom right).
In the text