© P. Cristofari et al. 2021

1. Introduction

Although there is strong evidence that particle acceleration takes place in supernova remnants (SNRs), it is still debated whether these objects can be the sources of all Galactic cosmic rays (CRs). The problem has various levels: First, the processes of particle acceleration and magnetic field amplification in an individual SNR depend on the type of explosion and the type of environment where it occurs. Second, the spectrum of particles accelerated at the shock and the one that the SNR releases in the surrounding interstellar medium (ISM) are, in general, quite different. The latter is typically made of two components, the one that escapes the remnant at any given time from the upstream region and the one that is advected downstream and is eventually liberated when the SNR dissipates. The particles trapped in the downstream region are affected by losses, which in general act differently on protons and electrons. Even this simple line of thought leads to two conclusions: (1) Although the instantaneous spectrum of particles accelerated at the SNR shock is a power law in momentum, the released spectrum does not need to be; (2) the spectra of protons and electrons from an SNR can, in general, be different. It is worth keeping in mind that nonlinear effects might lead to a slight deviation from the perfect power laws predicted in linear theory (Reynolds & Ellison 1992; Malkov & Drury 2001).

As previously discussed by Cristofari et al. (2020), the spectrum of CR protons from different types of SNRs is hardly a pure power law, and it extends to a maximum momentum that depends rather critically on the type of SNR and on the environment in which the explosion takes place. The structure in the spectrum is due to the contribution of different periods during the evolution of the remnant as well as the overlap of the advected and escaped fluxes. For type Ia supernovae (SNe) the proton spectrum extends to ≲100 TeV, but it is characterized by a marked steepening at around teraelectronvolt energies. Such a steepening is associated with the transition between the spectrum of advected particles and that of escaped CRs. The maximum energy at any given time was calculated using the growth rate of the nonresonant hybrid instability (Bell 2004). For core collapse SNe, higher maximum energies can in principle be reached, but the steepening at the transition discussed above is so pronounced that in fact the flux at the Earth gets strongly suppressed at energies of a few TeV. For very energetic SN explosions (≳5 × 10⁵¹ erg) taking place in a dense pre-SN wind (here called type II* SNe) and having an estimated rate of occurrence on the order of a few per 10⁴ years, the normalization of the flux approaches the measured CR flux (at all energies) and the maximum energy is around the knee. The spectral shape of CRs contributed by each of these SN types is characterized by structures that appear qualitatively similar to the bumps recently observed by DAMPE (An et al. 2019).

In the present work we focus on the description of magnetic field amplification and its implications for the maximum energy reached by protons and electrons as well as radiative losses of electrons trapped in the downstream plasma of an SNR. In particular, we provide a detailed description of the adiabatic losses and escape of protons from SNRs of different types, as well as the transport of electrons subject to synchrotron losses. This last calculation was recently presented by Diesing & Caprioli (2019) for a typical remnant from a type Ia progenitor expanding in a uniform ISM; the authors concluded that the spectrum of electrons is systematically steeper than that of protons, in line with the requirement arising from calculations of the transport of nuclei and electrons in the Galaxy (Evoli et al. 2020). However, the finding of Diesing & Caprioli (2019) was based on a recipe for magnetic field amplification derived by Amato & Blasi (2006) for resonant streaming instability (but without the natural saturation to δB ∼ B that is appropriate for these modes) and modified with a phenomenological recipe, such that the Alfvénic Mach number is replaced by the same quantity calculated in the amplified field; such a prescription was used, for example, to reproduce the multiwavelength emission from the Tycho SNR (Morlino & Caprioli 2012; Slane et al. 2014). This prescription leads to relatively large magnetic field amplification at late times and a maximum energy of electrons that remains loss-dominated even for old SNRs. We show that the difference in spectral shape between protons and electrons is highly sensitive to the strength of the magnetic field in such late phases of SNR evolution, and the effect virtually disappears if magnetic field amplification is described solely by Bell instability. We carry out this calculation for three types of SNRs, for both protons and electrons, in order to assess the role of the environment around SNRs for the shaping of CR spectra.

The article is organized as follows: In Sect. 2 we summarize our understanding of CR-induced magnetic field amplification at SNR shocks and compare the predicted magnetic field with a compilation of observational results for a number of SNRs. In Sect. 3 we describe how we follow protons and electrons in the downstream region of the shock, and we summarize the description of the escape of these particles from the upstream region. In Sect. 4 we briefly summarize our treatment of the temporal evolution of different types of SNRs in the surrounding environment. In Sect. 5 we illustrate the main result of our calculations in terms of the injection spectra of protons and electrons from different types of SNRs, integrated in time over the whole temporal evolution of the SNR shock through the surrounding medium. Finally, in Sect. 6 we discuss the implications of our results.

2. Magnetic field amplification in SNRs

In this section we summarize our current understanding of CR-induced magnetic field amplification at SNR shocks. Magnetic field perturbations can be produced in the shock proximity due to a variety of processes, but only some of them lead to important effects in terms of particle scattering. For instance, the propagation of a shock front in a medium with density inhomogeneities can excite a Richtmeier-Meshkov instability (Giacalone & Jokipii 2007) that leads to the growth of perturbations downstream of the shock on a timescale of order ℓ/v_A, where ℓ is the spatial scale of the density inhomogeneities upstream and v_A the Alfvén speed. Although such a magnetic field may be important in terms of determining the morphology of synchrotron emission from a remnant and the strength of the downstream magnetic field, it does not appreciably affect the diffusion time of accelerated particles upstream, and hence it does not lead to a substantial increase in the maximum energy that can be reached through diffusive shock acceleration (DSA). A fundamental step forward in the investigation of the interaction between CRs and the surrounding medium has been made with the discovery of the nonresonant hybrid instability (Bell 2004), which is expected to be excited upstream of a shock due to the accelerated particles themselves. This is a current-driven instability, excited both by CR particles that are leaving the acceleration region and by CRs diffusively confined in the vicinity of the shock front. If the density of CRs with momentum > p at the shock is n_CR(> p), the electric current that these particles carry is J_CR(> p) = eD(∂n_CR/∂z)_shock = ev_shn_CR(> p) and it extends over a precursor distance ∼D(p)/v_sh. On the other hand, the density of CRs escaping toward the upstream infinity is limited to the highest energy particles and can be estimated as ∼n_CR(> p)(v_sh/c). Hence, the corresponding current is J_CR(> p) ≈ en_CR(v_sh/c)c = en_CRv_sh, numerically equivalent to that estimated above for the same momentum despite the fact that escaping particles are assumed to be streaming away ballistically (namely, moving at roughly the speed of light, c). If the differential spectrum of accelerated particles at the shock is , the normalization can be easily found by requiring that the CR pressure be a fraction, ξ_CR, of the ram pressure at the shock location, :

(1)

which implies that , with . We notice that the normalization constant defined in this way depends very weakly on the minimum and maximum momenta, p_min and p_max, provided 4 ≤ α < 5, as expected for particle acceleration by DSA. In particular, for α = 4 one has I(α) ≈ ln(p_max/m_pc). The spectrum of accelerated particles may be outside this range only if nonlinear effects due to the CR pressure lead to the formation of a precursor upstream, which in turn may lead to spectra harder than p⁻⁴ (e.g., Jones et al. 2001; Malkov & Drury 2001). However, in practice, strong spectral modification should not be expected, because of numerous other effects (Berezhko & Ellison 1999; Caprioli et al. 2009a). Moreover, recent self-consistent kinetic simulations of strong shocks suggest that the formation of a shock post-cursor naturally leads to spectra with α between 4 and 5 (Caprioli et al. 2020), consistent with γ-ray observations of SNRs (Caprioli 2011). In general, the slope depends on the shock Mach number (if the shock is not strong), but even at the end of the Sedov-Taylor (ST) phase the Mach number remains much larger than unity (see Sect. 4).

The current carried by accelerated particles with momentum > p is

(2)

Here we especially focus on the scenario where nonresonant hybrid instability is excited by escaping particles since this channel leads to the formation of magnetic perturbations far upstream. As discussed by Bell (2004), the fastest growing mode is associated with a wavenumber, k_max, and can be written as γ_max = k_maxv_A, where is the Alfvén speed in the unperturbed field, B₀. The wavenumber k_max is determined based on the condition:

(3)

and the excitation of the instability occurs only if k_max > r_L(p) = pc/eB₀, which translates to the following constraint:

(4)

In other words, the instability is excited if the energy density in the form of escaping particles is larger than that of the preexisting magnetic field. As discussed in many previous works, there are different approaches to the saturation of the instability. The most intuitive one, based on comparing the plasma displacement due to the J_CR × δB/c force with the Larmor radius of particles in the amplified field, δB, leads to the conclusion that the field stops growing when the energy density in the form of escaping particles equals that in the amplified field:

(5)

At saturation, the spatial size of the perturbations becomes comparable with the Larmor radius in the amplified field, such that the current is now disrupted because of efficient CR scattering off the self-generated perturbations, thereby causing the drive for magnetic field amplification to stop.

Some comments on Eq. (5) are in order: (1) From the point of view of the scattering of particles with momentum p, the diffusion coefficient is D(p) ∝ p/δB. If the spectrum of accelerated particles is ∼p⁻⁴, then the diffusion coefficient is Bohm-like (linear in momentum). In a general case, α = 4 + ϵ, the diffusion coefficient turns out to be D(p) ∼ p^1 + ϵ/2. Since typically ϵ ∼ 0 ÷ 0.3, the expected deviations from Bohm-like behavior are small. (2) With the exception of the case ϵ = 0, Eq. (5) shows a weak dependence upon the momentum p where the current is calculated. (3) In terms of the magnetic field immediately upstream of the shock, all particles should be included in estimating the magnetic field in Eq. (5). For α = 4 (ϵ = 0), the resulting field does not depend on this choice, but for ϵ > 0 the resulting magnetic field shows a weak dependence on the minimum momentum. Assuming a minimum momentum χm_pv_sh, with χ > 1, one has:

(6)

If the turbulent field upstream of the shock is roughly isotropic and the perpendicular components are compressed at the shock, with compression factor r, the mean value of the compressed downstream magnetic field is . For r = 4, the compression factor is . It follows that:

(7)

For the parameter values appropriate for the Tycho SNR, this expression returns a downstream magnetic field that ranges between 150 μG (α = 4) and 400 μG (α = 4.3), in good agreement with the value inferred from multiwavelength observations (e.g., Morlino & Caprioli 2012).

In Fig. 1 we compare the quantity in Eq. (7) with the corresponding value measured in several SNRs as reported by Vink (2012). The solid and dotted lines show the result of Eq. (7) for the cases α = 4.3 and α = 4, respectively, assuming χ = 3. When the amplified field upstream becomes smaller than B₀, the magnetic field relevant for synchrotron losses becomes of order . For typical values of parameters in the ISM, this is reflected as erg g⁻¹ in Fig. 1. This condition also identifies the end of the stage where nonresonant modes get excited (see also Amato & Blasi 2009) and typically occurs when v_sh ≲ 1000 km s⁻¹.

Fig. 1. Magnetization of the downstream region of SNR shocks as a function of shock velocity. The data points are from Vink (2012). The lines refer to the growth of nonresonant modes with an injection spectrum with slope α = 4 (solid red) and α = 4.3 (dotted red) and to the recipe adopted by Diesing & Caprioli (2019; dashed green), which is independent of the spectrum of accelerated particles.

Equation (5) clearly shows how, for this type of instability, when saturation is reached the quantity δB²/8πρ scales with the shock velocity as ( for α = 4.3) and is independent of the strength B₀ of the preexisting magnetic field. Other recipes for saturation (see for instance the one proposed by Riquelme & Spitkovsky 2009) suggest some weak dependence of the magnetic field at saturation on B₀. These recipes lead to a somewhat less amplified magnetic field, a point to keep in mind in the discussion below.

The dashed green curve in Fig. 1 represents the same quantity calculated following the recipe originally put forward by Caprioli (2011) and used in Morlino & Caprioli (2012). A few comments are in order concerning the genesis of this approach. Resonant streaming instability was originally included in nonlinear theories of DSA by Amato & Blasi (2006), where the equation for the growth of these modes was solved analytically. The saturation level derived from this simple approach can be written as:

(8)

where M_A = v_sh/v_A, 0 is the Alfvén Mach number calculated with respect to the Alfvén speed in the field B₀, . It follows that the compressed field downstream reads:

(9)

It is worth noting that, since M_A = v_sh/v_A, 0, effectively . In Fig. 1, this is shown as a dashed purple line. We should stress that this amplified field should be used only when the Bell modes are not allowed to grow (low shock velocity). In the opposite condition, the growth rate of the resonant instability is proportional to and leads to a lower magnetic field than the nonresonant instability, and usually δB/B₀ ≲ 1 (Zweibel 1979; Amato & Blasi 2009).

Morlino & Caprioli (2012) proposed a formal modification of Eq. (9) so as to mimic the onset of nonlinear effects: The proposed modification consists in attempting to interpret the Mach number in Eq. (8) as the one calculated with respect to the Alfvén speed in the amplified field, δB. This assumption results in a modified expression for the amplified field that, after compression at the shock, reads:

(10)

This quantity is plotted in Fig. 1 for r = 4 (dashed green line). The most distinctive characteristic of this trend is the scaling with , which is quite different from the typical of nonresonant hybrid modes (see Eq. (7) for α = 4). Equation (7) results in a much larger magnetic energy density downstream at late times, when the shock velocity drops.

Though being a viable phenomenological approach, this is based on a few assumptions that are worth keeping in mind: (1) Although based on a perturbative approach, Eq. (8) leads to B₁/B₀ ≫ 1 for high shock speeds, which is unphysical if modes are strictly resonant. (2) As pointed out above, the growth rate used to describe resonant modes, and which led to Eq. (8), is only valid in the case in which the energy density of particles is much smaller than , which is certainly not the case for fast shocks (in fact, this is the very regime where nonresonant hybrid modes grow). This was discussed at length by Amato & Blasi (2009, see also Blasi 2013).

To add to this rather confused situation, hybrid kinetic simulations of particle acceleration at shocks lead to a prescription for the amplified field that reads (Caprioli & Spitkovsky 2014):

(11)

which is formally similar to Eq. (8) but with a different numerical factor. It is not easy to encapsulate these simulation results in the above theory: The simulations were run for a range of parameters in which Bell modes would grow and are seen to be growing. On the other hand, these simulations are nonrelativistic, which implies that the anisotropy of the accelerated particles (which for an astrophysical shock is ∼v_sh/c ∼ 10⁻² ÷ 10⁻³) is ∼v_sh/v ∼ 0.1. This might lead to a larger fraction of the CR energy being channeled into CRs compared with an SNR shock. Moreover, if accelerated particles are forced to remain nonrelativistic, the growth of the Bell instability can be shown to saturate to a value that is larger than the one in Eq. (6). In any case, no dependence upon the background field B₀ is expected, based on semi-analytical arguments, unless more complex phenomena come into play.

Finally, the magnetic field energy density in Eq. (11) shows a scaling , but the nonrelativistic simulations discussed by Caprioli & Spitkovsky (2014) were run for a given shock velocity and changing B₀ so as to achieve different Alfvén Mach numbers. This procedure is optimal for unraveling the scaling of B² with v_sh/v_A but would not reveal the additional scaling with the actual CR speed, c/v_A. These aspects definitely deserve further investigation.

In the following we focus on the investigation of the two prescriptions for magnetic field amplification described in Eqs. (7) and (10). We also comment on the case in which the field is amplified through nonresonant streaming instability when the appropriate condition is satisfied and by resonant streaming instability (Eq. (8)) at later times.

A quick inspection of Fig. 1 shows that data points are too sparse to allow us to infer a clear dependence of the amplified field on shock speed, v_sh, although there is a mild preference for a trend. As one can see from Eq. (7), the nonresonant hybrid modes are expected to lead to , which compares well with the observed trend for α = 4 ÷ 4.3. As discussed below, most of the modification of the electron spectrum due to energy losses occurs at late times (low shock speeds). Hence, for the purpose of calculating the difference in spectrum between electrons and protons, it is of critical importance to understand what the downstream magnetic field for older SNRs is.

The importance of the nonresonant hybrid instability for particle scattering is well known and will be briefly summarized here only for the sake of completeness. Following Bell et al. (2013) and Cristofari et al. (2020), the maximum energy of protons can be estimated by requiring the condition , which leads to:

(12)

Introducing the expression for the total magnetic field upstream, B₁, one can rewrite this expression as

(13)

where is the cyclotron frequency and is the Alfvén speed, both calculated in the amplified magnetic field, B₁. It is worth pointing out that for α = 4.3 the total magnetic field at the shock is slightly larger than for α = 4, but the corresponding p_max is lower, because there is less power available at the scales resonant with particles with momentum p_max. The maximum momentum of electrons is computed by equating the acceleration time to the minimum between the synchrotron time and the age of the SNR (see, e.g., Blasi 2010).

3. Cumulative spectra of CRs at SNRs

The contribution of CRs from an individual SNR can be written as the sum of two components: (1) particles accelerated at the SNR shock and trapped downstream of the expanding shell until the time when the shock dissipates away and the particles are released into the ISM after having suffered the effect of losses, N_loss, and (2) particles accelerated at the shock up to the highest energy that can be achieved at that time and leaving the acceleration region moving upstream, N_esc. Although the latter contribution is strongly peaked around the maximum energy reached at that time, the integration over the whole expansion history of the remnant leads to a continuous spectrum of CRs released into the ISM (Caprioli et al. 2009b).

The spectrum of particles accelerated at a strong shock in the test-particle limit (Caprioli et al. 2009b), assuming a free escape boundary condition at some location upstream, reads:

(14)

For p ≪ p_max(t), f^p(p) ∝ (p/p_inj)^{−3r/(r − 1)}, and for p ≫ p_max(t), .

Working under the usual assumption that a fraction of the ram pressure of the SNR shock is converted into CRs, it is easy to write the spectrum of CRs accelerated at the SNR shock at any given time. In a given time interval, dt, the number of particles of momentum p dn_acc(p, t) accelerated at the shock reads:

(15)

where v_sh/r is the velocity downstream of the SNR shock and r is the compression factor at the shock.

If accelerated particles could escape the acceleration region immediately after penetrating downstream, the total number of particles integrated over the entire active lifetime of an SNR, from t₀, typically the beginning of the free expansion phase, to T_SN, the end of the ST phase, would read:

(16)

This quantity could be interpreted as the flux of CR particles contributed by an SNR in the absence of adiabatic energy losses. The departures from such a spectrum will be used later as an index of the effect of losses on the spectrum of CRs from each SNR. In principle, if the SNR shell is broken in some locations where escape is allowed, the time-integrated CR flux would be somewhat similar to that estimated using Eq. (16).

The spectrum of electrons at the SNR shock is calculated as in Morlino et al. (2009), using the approximated expression proposed by Zirakashvili & Aharonian (2007):

(17)

where K_ep is the electron-to-proton ratio, typically in the range 10⁻⁴ − 10⁻². The cumulative spectrum of electrons, , can then be calculated as in Eq. (16), using f^e instead of f^p.

Below, we calculate the spectrum of protons and electrons that leaves an SNR after accounting for the effect of adiabatic and radiative energy losses.

3.1. Adiabatic and radiative losses of CRs downstream of the shock

As the SNR shock expands, the particles produced at the shock, and trapped inside the SNR downstream of the shock, suffer adiabatic and radiative losses. The latter are dominated by the emission of synchrotron photons, while inverse-Compton scattering, although included in our calculation, is typically negligible. In fact, the rate of inverse-Compton losses on cosmic microwave background photons is the same as that of synchrotron losses in a ∼3 μG magnetic field, and post-shock fields are much larger due to amplification and compression. On the other hand, inverse-Compton losses might become important for SNRs located in star-forming regions if the energy density of optical-to-infrared photons exceeds ∼100 eV cm⁻³.

The number of particles (protons or electrons) that are liberated by an individual SNR at the end of the evolution can be easily written in terms of the conservation of the total number of particles. In fact, the number of particles with momentum p at the end of the SN evolution (t = T_SN) is the result of all the particles produced at earlier times (t < T_SN) with momentum p′> p, such that in a time T_SN the momentum has degraded down to p. We can then write:

(18)

The change in momentum of a particle injected at time t′ with momentum p′ due to losses is:

(19)

where σ_T is the Thompson cross section and ℒ accounts for adiabatic energy losses in terms of the change in volume between the two times, t′ and t:

(20)

where ρ_down is the density downstream of the shock. If the expansion is adiabatic, , where ρ is the gas density upstream of the shock. For protons, synchrotron losses are negligible, while for electrons both adiabatic and radiative losses are important.

3.2. Escaping particles

The flux of particles escaping the accelerator from the upstream region can be written following Caprioli et al. (2009b, see also Cristofari et al. 2020):

(21)

where we introduced the function

(22)

which describes the spectral shape of the particles that can escape upstream at a given time. For protons, the function G(p, t) is strongly peaked around the maximum momentum of protons, p_max(t), at the given time, t. For electrons, when the maximum energy is determined by diffusion, the meaning of this function is the same as for protons. When the magnetic field is large enough that synchrotron losses dominate the maximum energy of accelerated electrons (namely for young SNRs), the function G is vanishingly small and no escape upstream is possible.

4. Evolution of the shock in the circumstellar environment

The evolution in time of SNR shocks directly impacts the spectra of particles injected by these sources into the ISM, when integrated on the SNR life span. In this section we briefly summarize the calculations adopted here for the description of such evolution, which is discussed in more detail in Cristofari et al. (2013).

For the sake of calculating the contribution to the CR flux, SNe are broadly classified into two groups, depending on the mechanism triggering the explosion: thermonuclear SNe (typeIa) and core-collapse SNe (type II). In addition, we consider a peculiar type of very energetic core-collapse SNe, which we call type II*. This type is introduced here just to illustrate the wide range of physical parameters that apply to type II SNe and to demonstrate that these energetic events can accelerate particles up to the petaelectronvolt range (when E_SN and Ṁ are sufficiently high and M_ej is sufficiently low; see the discussion on the parameter space presented by Cristofari et al. 2020). Schematically, SNR shocks from type Ia SNe expand in a uniform ISM. The evolution of the shock radius, R_sh, and shock velocity, v_sh, in the ejecta-dominated phase and the ST phase are well described by self-similar solutions (see, e.g., Chevalier 1982; Truelove & McKee 1999; Ptuskin & Zirakashvili 2005). Here we rely on the approach presented in Truelove & McKee (1999, 2000), adopting the same formalism as in Cardillo et al. (2015):

(23)

where a = −5, , λ_FE = 4/7, λ_ST = 2/5, and . This description holds until the end of the ST phase (i.e., the beginning of the radiative phase) of SNR evolution. This transition typically occurs when the age of the SNR becomes on the order of the cooling time: kyr (Blondin et al. 1998). In this work, we assume E_SN = 10⁵¹ erg, M_ej = 1 M_⊙, and n₀ = 1 cm⁻³ as reference values for type Ia SNe. It is worth recalling that in the presence of efficient CR acceleration at the SNR shock, as discussed by Diesing & Caprioli (2018), the evolution in time of the shock in the final stages of the evolution may be affected rather remarkably by the CR pressure, and the beginning of the radiative phase may be delayed. We do not include these effects here.

In general, type II SNR shocks expand in a complex medium structured by the evolution of the massive progenitor star before the explosion of the SN. Throughout its main sequence, the stellar wind produces a low-density hot-temperature bubble that is in pressure equilibrium with the ISM. Later, when the massive star reaches the final stages of stellar evolution, typically entering the red super giant (RSG) stage, a slow dense wind forms. Therefore, after the SN explosion, the shock expands through a dense wind and then through a low density bubble before finally reaching the ISM. The density of the dense wind created by the RSG is typically n_w = Ṁ/(4πmu_wr²), were Ṁ ∼ 10⁻⁵ M_⊙ yr⁻¹ is the mass-loss rate, u_w ∼ 10⁶ cm s⁻¹ the velocity of the wind, and m = m_p(1 + 4f_He)/(1 + f_He) ∼ 1.27m_p the mean mass of the ISM per hydrogen nucleus (Weaver et al. 1977). The density of the low density bubble is cm⁻³, with t_Myr the duration of the main sequence, which is on the order of a megayear (Longair 1994). The transition between the dense RSG wind and the low density cavity, r₁, is set by equating the RSG wind pressure to the thermal pressure of the hot cavity: , where k is the Boltzmann constant. The radius of the hot bubble is pc, where L₃₆ is the main sequence star wind power in units of 10³⁶ erg s⁻¹. In the case of type II SNe, the evolution in the structured medium described above does not allow for direct self-similar solutions. It is, however, possible to work under the “thin-shell” approximation, considering that the swept-up gas is located in a think layer behind the shock wave (see, e.g., Ostriker & McKee 1988; Bisnovatyi-Kogan & Silich 1995). In the case of spherically symmetric distribution, Ptuskin & Zirakashvili (2005) obtained:

(24)

and

(25)

For a typical type II SN, we assume E_SN = 10⁵¹ erg, M_ej = 5 M_⊙, and Ṁ = 10⁻⁵ M_⊙ yr⁻¹. For a type II*, we assume E_SN = 5 × 10⁵¹ erg, M_ej = 1 M_⊙, and Ṁ = 10⁻⁴ M_⊙ yr⁻¹.

5. Results

In this section we discuss our results in terms of the spectrum of protons and electrons injected into the ISM by different types of SN explosions. As discussed above, the CR spectrum released by an individual source is generally the sum of (1) the contribution due to particles that escaped upstream at any given time and (2) the contribution of CRs trapped downstream, where energy losses are in action, and liberated in the final stages of the SN evolution, which we assume coincide with the beginning of the radiative phase. While the proton escape flux is always present, the electron escape flux vanishes when radiative losses limit the maximum energy. It is this phase of radiative and adiabatic losses that can potentially make the overall spectrum of electrons different from that of protons if radiative losses are sufficiently severe.

The spectrum of protons liberated by type Ia, II, and II* SNRs is shown in Fig. 2, with the corresponding slopes shown in the bottom panels. As discussed by Cristofari et al. (2020), the highest energies are reached at very early times, but they do not reflect the equally large energies in the final spectrum of protons because the amount of mass that the SN shock processes at such early times is very small. In fact, the effective maximum energy is the one reached at the beginning of the ST phase. The easiest case is that of type Ia SNRs, where the explosion occurs in the normal ISM, with a spatially constant gas density and background magnetic field. For type Ia SNRs the effective maximum energy is a few tens of TeV (left panel of Fig. 2). There is an additional spectral steepening at somewhat lower energies due to the temporal evolution of the maximum energy. More specifically, the steepening occurs at the maximum energy reached at the end of the ST phase, typically a few TeV. The flux of escaping CR protons starts at about the same energy, as is clearly visible in Fig. 2.

Fig. 2. Spectra of protons produced at SNRs from type Ia (top), type II (center), and type II* (bottom) SNRs for α = 4 (thick lines) and α = 4.3 (thin lines) if they were instantaneously liberated into the ISM (broken shell assumption). The dashed curves illustrate the effect of adiabatic losses in the downstream region, while the dotted lines refer to the escape flux from the upstream region. In the bottom part of each panel we also show the local slope of the spectrum q(p) at a given momentum.

For a strong shock, such as the one expected for a young SNR expanding in the normal ISM, the spectrum of accelerated particles at the shock location has a slope very close to 4 (thick lines in Fig. 2). Nevertheless, as recently discussed by Caprioli et al. (2020), the spectrum can be steeper if the finite velocity of scattering centers in the downstream plasma is taken into account. For this reason, in Fig. 2 we also show the case α = 4.3 (thin lines). In all cases of interest, the spectra of CR protons that are injected into the ISM (as the sum of the two contributions) are quite close to the spectrum at the shock in terms of slope, with the exception of the highest energies, as discussed above.

For type II SNRs, the spectrum of CR protons is shown in the middle panel of Fig. 2. For the sake of making a fair comparison between the three types of SN explosions, here we used an acceleration efficiency of ξ_CR = 0.1 for all of them. As discussed by Cristofari et al. (2020), because of the different rates of occurrence of these events in the Galaxy, for type II SNRs the efficiency is required to be somewhat lower than for type Ia, which is also reflected in a lower value of the maximum energy of particles accelerated at the shock (see Eq. (12)). Despite this bias, the maximum achievable energy for type II SNRs remains on the order of ∼10⁵ GeV and falls short of the knee by a large amount, as already pointed out by Cristofari et al. (2020).

Only when parameters are pushed to the extreme (what we have called here type II* SNRs) can the maximum energy reach the knee, as shown in the right plot of Fig. 2. As already pointed out by Caprioli et al. (2009b), the superposition of the escape flux from the different stages of shock evolution in the complex environment around these SNRs may lead to the appearance of bumps in the overall CR spectrum that might be related to the feature recently measured by DAMPE in the 10−100 TeV region of the proton spectrum (An et al. 2019).

The corresponding spectra of electrons injected by SNRs of different types into the ISM are shown in Fig. 3. The thick and thin curves refer to α = 4 and α = 4.3, respectively. The dash-dotted line identifies the spectrum of particles accelerated at the shock, as if they were immediately liberated into the ISM, without energy losses. The solid lines are the spectra of electrons liberated into the ISM after adiabatic and synchrotron losses downstream of the shock, while the upstream escape flux, limited to the times when the maximum energy of electrons is not determined by energy losses, is shown in the form of dotted lines. If the SNR shell were broken or if confinement in the downstream region were energy-dependent (e.g., due to turbulence damping), the actual contribution would lie between the dash-dotted and solid lines.

Fig. 3. Spectra of electrons produced at SNRs from type Ia (top), type II (center), and type II* (bottom) progenitors for α = 4 (thick lines) and α = 4.3 (thin lines). The dash-dotted line indicates the spectrum of particles accelerated at the shock, the solid lines the spectra of electrons liberated into the ISM after losses downstream of the shock, and the dotted lines the upstream escape flux. In the bottom part of each panel we also show the slope of the spectrum q(p) at a given momentum.

The rate of synchrotron losses is larger when the condition for the growth of the magnetic field through the excitation of the nonresonant hybrid instability is fulfilled. As discussed in Sect. 2, for this instability, and hence the mechanism becomes less effective or even ineffective in the late stages of SNR evolution; these stages are, however, crucial for the production of low energy electrons. As a consequence, the effect of radiative energy losses is only important at energies at or above teraelectronvolt levels, while it is minor at lower energies, as clearly shown in Fig. 3, independent of the type of SNR considered. It follows that if the magnetic field amplification is mainly due to the excitation of Bell modes, then losses cannot be the main reason for a difference in the spectra of protons and electrons in SNRs.

This result is illustrated even more clearly in Fig. 4, where we show the difference in slope between protons and electrons as a function of energy for type Ia (left), type II (middle), and type II* (right) SNRs. The thick and thin curves refer to the cases α = 4 and α = 4.3, respectively. The dashed lines were obtained using the growth of nonresonant hybrid modes to calculate the magnetic field, while the dotted lines are based on Eq. (10), the same as in Diesing & Caprioli (2019). The reason for the very different results is that the two recipes lead to quite different predictions for the downstream magnetic field at late times, which are the most important times for determining the electrons’ spectrum: The nonresonant instability leads to an intense field in the early phases when the shock is faster, but the instability virtually shuts off at later times when the shock velocity drops below ∼1000 km s⁻¹. On the other hand, the prescription used by Diesing & Caprioli (2019) leads to larger magnetic fields being retained even in late stages, such that the acceleration of electrons remains loss-dominated at such late times. This is reflected in steeper electron spectra. The main lesson to be learned from this exercise is that the spectrum of electrons liberated by SNRs into the ISM is strongly dependent upon knowledge of the magnetic field in the phases of SNR evolution about which we know the least.

Fig. 4. Difference between the proton and electron spectral index at SNRs from type Ia (top), type II (center), and type II* (bottom) progenitors for α = 4 (thick lines) and α = 4.3 (thin lines). The dash-dotted line indicates the spectrum of particles accelerated at the shock, the solid lines the spectra of electrons liberated into the ISM after losses downstream of the shock, and the dotted lines the upstream escape flux. In the bottom part of each panel we also show the local slope of the spectrum q(p) at a given momentum.

The issue boils down to whether we have observational bounds or theoretical arguments that can be used to assess the credibility of different scenarios. From the observational point of view, to our knowledge, there is no evidence in favor of or against having relatively large magnetic field amplification at late times. From the theoretical point of view, the only prejudice that we can mention is that the recipe based on Eq. (10) is somewhat weak, for the reasons explained above. To be somewhat more conservative and limit our attention to scenarios that are based on some sort of theoretical ground, one could say that the magnetic field downstream is provided by Eq. (7) when the energy density of accelerated particles satisfies the condition for the growth of the nonresonant instability, and by Eq. (9) when the nonresonant instability stops growing and only resonant Alfvén modes get excited. In both cases, when the amplified magnetic field becomes much smaller than B₀, the field to be used for the purpose of calculating the electrons’ energy losses reduces to B₀.

To further address the importance of the late times in SNR evolution for shaping the electron spectrum, we considered a rather extreme situation in which particle acceleration is efficient until 80 kyr. Indeed, it has been shown that the duration of the ST phase can be substantially extended, for instance because of the CR pressure on the SNR shock (Diesing & Caprioli 2018). In Fig. 5, we illustrate the results for a type Ia SN (for a slope 4.3 of the spectrum of accelerated particles): Although the effect of losses is more pronounced (in the Bell and Bell-plus-nonresonant cases), it still remains limited to high energies, while no systematic difference is visible in the range 10−1000 GeV.

Fig. 5. Spectra of electrons produced at SNRs from type Ia progenitors for α = 4.3. We illustrate the effect of the different prescriptions for the magnetic field amplification and an increased duration of the ST phase, TSN = 80 kyr (thin lines). The dash-dotted line indicates the spectrum of particles accelerated at the shock, the solid lines the spectra of electrons liberated into the ISM after losses downstream of the shock, and the dotted lines the upstream escape flux. In the bottom part of each panel we also show the slope of the spectrum q(p) at a given momentum.

6. Discussion

Supernova remnants contribute a spectrum of CRs comprising two terms: one is the flux of particles that escape from upstream of the shock and the other is the flux of particles that are liberated into the ISM when the shock is dissipated and the particles initially trapped downstream are free to leave the remnant. For protons, the escape flux, at any given time, is strongly peaked around the maximum energy at that time. The integration over the whole history results in a broad spectrum that roughly extends from the maximum momentum reached at the end of the ST phase and the maximum momentum reached at the beginning of the same phase. On the other hand, protons that are trapped in the downstream region lose energy adiabatically, and these losses are reflected in the normalization and shape of the CR spectrum liberated into the ISM.

For electrons, the situation is somewhat more complex: Electrons are sensitive to energy losses, both through the maximum energy that they can achieve and because of the radiative losses that they suffer while being advected downstream. Their escape from the upstream region is similar to that of protons when their maximum energy is not limited by losses, while it drops to zero when radiative losses are the main limitation to electron acceleration.

We have described both protons and electrons, and both their escape flux from the upstream region and at the end of the evolution of the SNR, for different types of SNRs. We calculated the magnetic field in the acceleration region using the excitation rate of the nonresonant hybrid instability induced by CR protons and its saturation level. The magnetic field produced by the streaming of CR protons also determines the maximum energy and the radiative losses of electrons.

Our main conclusions can be summarized as follows:

(1) In the context of nonresonant streaming instability, which is the mechanism expected to drive magnetic field amplification in the early stages of the evolution of SNRs, only very powerful and rare events, resulting from the explosion of energetic SNe (what we have called type II* SNRs), can accelerate protons up to the PeV range. In these cases, a pronounced feature appears in the total spectrum around ∼100 TeV, where the spectrum liberated at the end of the SNR evolution merges with the escape flux integrated over time. A weaker feature is present at ∼10 TeV because of the complex environment in which these SN explosions take place. For type Ia SNe, the maximum energy of protons as derived in our calculations is ∼100 TeV, but there is a pronounced dip in the spectrum at ∼10 TeV. For type II SNe, the spectrum is strongly suppressed at energies around a few tens of TeV. In all cases, the spectra are power laws at low energies.

(2) The maximum energies of protons as discussed above are mainly the result of the self-confinement of CRs in the shock proximity due to the excitation of the current-driven nonresonant streaming instability (see also Cristofari et al. 2020). We calculated the strength of the magnetic field downstream of the shock that would be expected based on this process as a function of the shock velocity, and we compared the results with the observed trend, as reported by Vink (2012). The scaling , typical of this instability, is in good agreement with the sparse data available and provides the correct normalization. For old SNRs, with typical velocities of v_s ≲ 1000 km s⁻¹, the nonresonant instability is quenched and the strength of the magnetic field is comparable with that in the surrounding medium, δB ≲ B₀. A recipe for magnetic field evolution often used in the literature was put forward by Caprioli (2011) and Morlino & Caprioli (2012). For typical speeds of SNR shocks, such a recipe leads to larger values of the amplified field, especially for older SNRs. As discussed in Sect. 2, this recipe has numerous caveats that are important to keep in mind.

(3) The amplified magnetic field described above plays a crucial role in shaping the spectrum of the accelerated electrons that are advected downstream because of radiative losses. In the cases we investigated, a substantial steepening of the spectrum was only obtained when the phenomenological prescription of Caprioli (2011) was adopted, which leads to results that are in qualitative agreement with those of Diesing & Caprioli (2019). When the magnetic field is assumed to be amplified according to the growth of the nonresonant (Bell 2004) and resonant streaming instability (Amato & Blasi 2006), the effect of losses on the electron spectrum typically reduces to a cutoff in the overall spectrum rather than a broad steepening.

The appearance of a steepening in the electron spectrum (with respect to the proton spectrum) depends on the time dependence of the maximum electron energy, . In the initial stages of SNR evolution (ejecta-dominated to the early ST phase), the maximum energy is determined by energy losses: In the few hundred μG typical of an SNR at the beginning of the ST phase, synchrotron losses degrade the particle energy to ∼10 GeV by the end of the remnant life. Later times (in the STphase) may contribute electrons with increasingly larger energies if the amplified magnetic field decreases sufficiently fast with the shock velocity, as is typically the case for all the above prescriptions. When becomes age-limited, as for protons, contributions from different times pile up, eventually leading to an overall cutoff in the teraelectronvolt range. In summary, sustained magnetic field amplification may provide a steepening in the electron spectrum, but specific recipes for the amplification determine whether the steepening extends over several orders of magnitude in energy or is instead limited to a narrow energy range, thereby appearing as a cutoff in the electron spectrum. Current data suggest that a global steepening must be present between ∼10 GeV and a few TeV (Evoli et al. 2020), although it is not clear that synchrotron losses are the cause of this phenomenon (Ohira et al. 2012).

This result shows that if SNRs are responsible for the observed spectrum of electrons, the late phases of their evolution, usually considered to be of little interest for CR acceleration, are in fact crucial in the sense that synchrotron energy losses may be at work for long periods. On the other hand, the streaming instabilities that are thought to be most effective in amplifying the magnetic field may not be able to provide the magnetic fields that, through synchrotron losses, would cause the required steepening. Since such stages are crucial to shaping the spectrum of electrons, it is worth asking whether there are physical phenomena that we did not take into account that might lead to larger (or smaller) fields at late times: One point to keep in mind is that the total compression factor at CR-modified shocks might be somewhat larger than 4 if CR acceleration remains efficient; for instance, in the simulations of Haggerty & Caprioli (2020), one has r ∼ 6−7. This can indeed lead to somewhat larger magnetic fields, as one can infer from Eq. (7). However, in order to have a sizeable effect on the spectrum of electrons, this should happen for low shock velocities, while it is typically expected to be more of a concern for fast shocks. One could also speculate that, in addition to CR-induced magnetic field amplification (which takes place upstream), the field could be further amplified downstream of the shock, perhaps through Richtmeier-Meshkov instability. As discussed above, the duration of the STphase, which may also be controlled by the nonthermal CR pressure (Diesing & Caprioli 2018), may affect the spectrum of electrons as well. Finally, we note that if shock acceleration were efficient even when the shock reaches a Mach number of a few, the spectrum of accelerated particles, which is affected by the compression factor, would become steeper. For typical values of the parameters, when the shock reaches the radiative phase, its Mach number is still ∼40, so it is possible that the CR spectrum may still evolve during the early radiative stage.

Clearly, there are also physical processes that lead to smaller magnetic fields other than those predicted above, thereby making the spectral modification for electrons even less effective. For instance, when the shock moves in the ordinary ISM, as is typically the case in the late stages of the evolution of an SNR, neutral hydrogen induces a strong level of ion-neutral damping, which substantially limits the growth of perturbations. These effects are not included in our calculations, nor in numerical simulations of DSA.

From the discussion above, it is clear that if the difference in the spectrum of electrons and protons is to be attributed to radiative losses in the downstream region of an SNR shock, the required conditions appear to be rather at odds with those that we would typically assume to exist. In particular, some rather efficient mechanism for magnetic field amplification should come into effect for late, slow moving shocks. A dedicated effort to investigate these stages is needed.

On the other hand, if the results of our investigation are taken at face value, then the spectral shape of electrons and protons liberated into the ISM by an individual SNR should be very similar. Hence, it would follow that the observed difference should be attributed to phenomena occurring after the particles have been released into the ISM. If the diffusion coefficient describing transport in the Galaxy is the same for the two species, as one would expect, the only possibility left open is that electrons and protons may develop different spectral shapes while propagating in the neighborhood of the source due to the large perturbations induced by the escaping particles (Schroer et al. 2020). This possibility is currently being investigated.

Acknowledgments

We thank Rebecca Diesing and Rino Bandiera for interesting discussions about electron spectra and radio emission from SNRs.The research of PB was partially funded through Grant ASI/INAF No. 2017-14- H.O. DC was partially supported by NASA (grants NNX17AG30G, 80NSSC18K1218, and 80NSSC18K1726) and by NSF (grants AST-1714658, AST-1909778, and PHY-2010240).

References

All Figures

	Fig. 1. Magnetization of the downstream region of SNR shocks as a function of shock velocity. The data points are from Vink (2012). The lines refer to the growth of nonresonant modes with an injection spectrum with slope α = 4 (solid red) and α = 4.3 (dotted red) and to the recipe adopted by Diesing & Caprioli (2019; dashed green), which is independent of the spectrum of accelerated particles.
In the text

Fig. 2. Spectra of protons produced at SNRs from type Ia (top), type II (center), and type II* (bottom) SNRs for α = 4 (thick lines) and α = 4.3 (thin lines) if they were instantaneously liberated into the ISM (broken shell assumption). The dashed curves illustrate the effect of adiabatic losses in the downstream region, while the dotted lines refer to the escape flux from the upstream region. In the bottom part of each panel we also show the local slope of the spectrum q(p) at a given momentum.

In the text

Fig. 3. Spectra of electrons produced at SNRs from type Ia (top), type II (center), and type II* (bottom) progenitors for α = 4 (thick lines) and α = 4.3 (thin lines). The dash-dotted line indicates the spectrum of particles accelerated at the shock, the solid lines the spectra of electrons liberated into the ISM after losses downstream of the shock, and the dotted lines the upstream escape flux. In the bottom part of each panel we also show the slope of the spectrum q(p) at a given momentum.

In the text

Fig. 4. Difference between the proton and electron spectral index at SNRs from type Ia (top), type II (center), and type II* (bottom) progenitors for α = 4 (thick lines) and α = 4.3 (thin lines). The dash-dotted line indicates the spectrum of particles accelerated at the shock, the solid lines the spectra of electrons liberated into the ISM after losses downstream of the shock, and the dotted lines the upstream escape flux. In the bottom part of each panel we also show the local slope of the spectrum q(p) at a given momentum.

In the text

Fig. 5. Spectra of electrons produced at SNRs from type Ia progenitors for α = 4.3. We illustrate the effect of the different prescriptions for the magnetic field amplification and an increased duration of the ST phase, TSN = 80 kyr (thin lines). The dash-dotted line indicates the spectrum of particles accelerated at the shock, the solid lines the spectra of electrons liberated into the ISM after losses downstream of the shock, and the dotted lines the upstream escape flux. In the bottom part of each panel we also show the slope of the spectrum q(p) at a given momentum.

In the text