On the Spectrum of Electrons Accelerated in Supernova Remnants

TL;DR

This paper models the spectrum of cosmic ray electrons accelerated by supernova remnants, highlighting how synchrotron losses cause electrons to have steeper spectra than protons, impacting theories of Galactic cosmic ray transport.

Contribution

It introduces a semi-analytic model accounting for synchrotron losses, revealing the spectral difference between electrons and protons in supernova remnants.

Findings

Electrons have a steeper spectrum than protons by 0.1 to 0.4 in spectral index.

The spectral difference affects models of Galactic cosmic ray transport.

Implications for the observed positron excess in cosmic rays.

Abstract

Using a semi-analytic model of non-linear diffusive shock acceleration, we model the spectrum of cosmic ray (CR) electrons accelerated by supernova remnants (SNRs). Because electrons experience synchrotron losses in the amplified magnetic fields characteristic of SNRs, they exhibit substantially steeper spectra than protons. In particular, we find that the difference between the electron and proton spectral index (power law slope) ranges from 0.1 to 0.4. Our findings must be reckoned with theories of Galactic CR transport, which often assume that electrons and protons are injected with the same slope, and may especially have implications for the observed "positron excess."

Equations

$$\eta =\frac{4}{3\sqrt{\pi }}(R_{\mathrm{sub}}-1){\xi }_{\mathrm{inj}}^3{e}^{-{\xi }_{\mathrm{inj}}^2},$$

$$L(t_0,t)\equiv [\rho (t)/{\rho }_0{]}^{1/3}=[v_{\mathrm{sh}}(t)/v_{\mathrm{sh}}(t_0){]}^{\frac{2}{3\gamma }},$$

$$P_{\mathrm{w}}(x)=\frac{B(x{)}^2}{8\pi {\rho }_{\mathrm{ISM}}{v}_{\mathrm{sh}}^2}=\frac{1+u(x)}{4M_{\mathrm{A}}(x)u(x)}{P}_{\mathrm{cr}}(x),$$

$$B_1=\frac{\sqrt{4\pi {\rho }_1}{u}_1{v}_{\mathrm{sh}}}{M_{\mathrm{A},1}}=\sqrt{\pi {\rho }_{\mathrm{ISM}}}{v}_{\mathrm{sh}}{\frac{P_{\mathrm{cr},1}(2-P_{\mathrm{cr},1})}{1-P_{\mathrm{cr},1}}}^{3/2},$$

$$f_{\mathrm{e}}(p)=K_{\mathrm{ep}}{f}_{\mathrm{p}}(p){\left[1+0.523(p/p_{\mathrm{e},\max }{)}^{9/4}\right]}^2{e}^{-p^2/p_{\mathrm{e},\max }^2},$$

$$\frac{dE}{dt}=-\frac{4}{3}\sigma_{\rm T}c\bigg{(}\frac{E}{m_{\rm e}c^{2}}\bigg{)}^{2}\frac{B_{2}^{2}}{8\pi}-\frac{E}{L}\frac{dL}{dt},$$

$$q\equiv -⟨\frac{d\log \varphi (E)}{d\log E}⟩,$$

Full text

††preprint: APS/123-QED

Introduction.— Developing a complete paradigm for the origin of Galactic cosmic rays (CRs) with energies up to $\sim 10^{8}$ GeV requires a detailed understanding of their acceleration and propagation. The best source candidates for such acceleration are supernova remnants (SNRs), which provide sufficient energetics and an efficient acceleration mechanism (Hillas, 2005; Caprioli et al., 2010a). Namely, particles are scattered by magnetic field perturbations, resulting in diffusion across the SNR forward shock and an energy gain with each crossing (Fermi, 1954; Krymskii, 1977; Axford et al., 1977; Bell, 1978; Blandford and Ostriker, 1978). This mechanism, known as diffusive shock acceleration (DSA), predicts power law energy distributions of CRs, $\propto E^{-q}$, where $q$ depends only on the shock dynamics and $q\to 2$ for strong shocks.

Once accelerated, CR protons and electrons diffuse through the Galaxy such that their spectrum is modified by escape from the Galaxy and, in the case of electrons, energy losses due to synchrotron and inverse-Compton scattering. Moreover, CR protons interact with protons in the interstellar medium (ISM) to produce secondary particles, most notably positrons and antiprotons. Thus, in this standard picture, one would expect the positron and antiproton spectra to follow that of their parent protons, modulo effects due to their subsequent escape and energy loss.

To put this analysis into more quantitative terms, consider CR protons and electrons injected by SNRs with spectra $\propto E^{-q_{\rm p}}$ and $\propto E^{-q_{\rm e}}$, respectively. In the standard paradigm for CR transport, the Galactic residence time of these particles scales as $E^{-\delta}$, with $\delta\sim 0.2-0.4$ (e.g., Aguilar et al. [AMS Collaboration], 2016; Lipari, 2019). Thus, we would expect the observed proton spectrum to go as $E^{-(q_{\rm p}+\delta)}$. Leptons also experience energy losses due interactions with the Galactic magnetic field (synchrotron) and radiation fields (inverse Compton). We therefore expect $N_{\rm e}\propto E^{-(q_{\rm e}+\beta)}$, where $\beta\gtrsim\delta$ reflects the effective spectral steepening due to a combination of escape and energy loss. Assuming positrons and antiprotons are secondaries produced in interactions between CR and ISM protons, the positron spectrum should scale as $E^{-(q_{\rm p}+\delta+\beta)}$, and the antiproton spectrum as $E^{-(q_{\rm p}+2\delta)}$.

A notable observation that appears to be in conflict with this paradigm is the “positron excess” observed by PAMELA (Adriani, [PAMELA Collaboration], 2013) and AMS-02 (Accardo et al. [AMS Collaboration], 2014). Both collaborations report a positron fraction, $\chi\equiv\Phi_{e^{+}}/(\Phi_{e^{-}}+\Phi_{e^{+}})$ where $\Phi_{e^{+}}$ and $\Phi_{e^{-}}$ are the positron and electron fluxes, that rises with energy. In the picture described above, $\chi\simeq\Phi_{e^{+}}/\Phi_{e^{-}}\propto E^{-(q_{\rm p}-q_{\rm e}+\delta)}$. Thus, under the standard assumption that $q_{\rm p}=q_{\rm e}$, $\chi\propto E^{-\delta}$ and should decrease with energy (see, e.g., Amato and Blasi, 2018, for a thorough review).

This discrepancy may be at least partially resolved if electrons are injected into the Galaxy with a steeper spectrum than protons (i.e., $q_{\rm e}>q_{\rm p}$). Such steepening is physically motivated, as electrons experience synchrotron losses during the acceleration process. Although the lifetime of a typical SNR is much shorter than the CR galactic residence time, CR acceleration leads to magnetic field amplification (Skilling, 975a; Bell, 1978, 2004; Amato and Blasi, 2009), producing magnetic fields hundreds of times stronger than that of the Galaxy (e.g., Völk et al., 2005; Caprioli et al., 2009). The result is that the synchrotron loss time in SNRs is generally shorter than the DSA timescale and the effects of synchrotron emission are non-negligible.

In this Letter, we use a semi-analytic model based on the solution of the Parker equation for the CR transport to calculate the CR proton and electron spectrum accelerated by typical SNRs, accounting for the effects of magnetic field amplification. We then use these spectra to estimate $q_{\rm e}$ and $q_{\rm p}$. This work represents the first calculation of the CR electron acceleration spectrum that self-consistently accounts for magnetic field amplification and thus the resulting synchrotron losses that take place within SNRs (see, e.g., Ohira et al., 2012; Berezhko and Ksenofontov, 2013, for examples of previous estimates). Our findings may have significant bearing on CR propagation models and the interpretation of observations such as the “positron excess.”

Let us now introduce the formalism that we use to model SNR evolution and CR acceleration.

Remnant Evolution— SNRs are evolved using the formalism described in Diesing and Caprioli (2018). More specifically, SNR evolution can be understood in terms of four stages: the ejecta-dominated stage, in which the mass of the swept-up ambient medium is less than that of the SN ejecta, the Sedov stage, in which the swept-up mass dominates the total mass and the SNR expands adiabatically, the pressure-driven snowplow, in which the remnant cools due to forbidden atomic transitions but continues to expand because its internal pressure exceeds the ambient pressure, and, finally, the momentum-driven snowplow, in which the internal pressure falls below the ambient pressure and expansion continues due to momentum conservation.

While we model SNRs through the end of the pressure-driven snowplow, the majority of CRs are accelerated during the transition between the ejecta-dominated and Sedov stages. The DSA timescale for CRs of energy $E=E_{\rm max}$ is given by $\tau_{\rm DSA}\approx D/v_{\rm sh}^{2}$ where $D$ is the diffusion coefficient and $v_{\rm sh}$ is the shock speed. Assuming Bohm diffusion (Caprioli and Spitkovsky, 2014a), $D(E)\propto r_{\rm L}\sim E/B_{2}$ where $r_{\rm L}$ is the Larmor radius and $B_{2}$ is the post-shock magnetic field. Thus, $E_{\rm max}\sim B_{2}v_{\rm sh}^{2}t$ and, during the ejecta dominated stage characterized by roughly constant velocity, $E_{\rm max}$ increases. During the Sedov stage, the shock slows down such that $v_{\rm sh}\propto t^{-3/5}$, meaning that $E_{\rm max}$ decreases with time, i.e., $E_{\rm max}\propto B_{2}(t)t^{-1/5}$ (Cardillo et al., 2015; Bell et al., 2013). Our results are therefore most sensitive to the adiabatic SNR stages.

To model SNR evolution, we use the analytical approximation for the ejecta dominated stage presented in Truelove and Mc Kee (1999). Once the swept up mass exceeds the ejecta mass and the Sedov stage begins, we transition to the thin-shell approximation, in which we assume that most of the mass resides in a thin layer that expands due to pressure in the hot cavity behind it (Bisnovatyi-Kogan and Silich, 1995; Ostriker and McKee, 1988; Bandiera and Petruk, 2004).

All SNRs are assumed to eject $M_{\rm ej}=1M_{\odot}$ (1 solar mass) with $E_{\rm SN}=10^{51}\rm erg$ into a uniform ambient medium of density $n_{\rm ISM}\in[10^{-2},10^{1}]\ \rm cm^{-3}$.

*Proton Acceleration—*Instantaneous proton spectra are calculated using the Cosmic Ray Analytical Fast Tool (CRAFT) a semi-analytical formalism described in Caprioli et al. (2010b); Caprioli (2012) and references therein (in particular, Amato and Blasi, 2005, 2006). CRAFT self-consistently solves the diffusion-convection equation (e.g., Skilling, 975a) for the transport of non-thermal particles in a quasi-parallel, non-relativistic shock, including the dynamical backreaction of accelerated particles and of CR-generated magnetic turbulence. CRAFT is quick and versatile, but achieves the same degree of accuracy as Monte Carlo and numerical methods (Caprioli et al., 2010).

Particles are injected into the acceleration mechanism following the prescription in Blasi et al. (2005), namely that ions with momentum greater than $\xi_{\rm inj}$ ($\sim$ a few) times the post-shock thermal momentum are promoted to CRs (“thermal leakage,” see Malkov, 1998; Kang et al., 2002). While kinetic simulations show that protons are injected via specular reflection and shock drift pre-acceleration rather than via thermal leakage (Caprioli et al., 2015), such a prescription is calibrated with self-consistent kinetic simulations to ensure continuity between the thermal and non-thermal distributions (Caprioli and Spitkovsky, 2014b). $\xi_{\rm inj}$ can be mapped onto $\eta$, the fraction of particles crossing the shock injected into DSA, via

[TABLE]

where $R_{\rm sub}$ is the subshock compression ratio (i.e., the ratio of the density immediately behind the shock to that immediately in front of it) Blasi et al. (2005). In this analysis, $\xi_{\rm inj}$ (and thus $\eta$) is left as a free parameter, which allows us to span a range of shocks where CRs are either test-particles or dynamically important.

Once the proton spectrum has been calculated at each timestep of SNR evolution, particle momenta are weighted by $1/L(t_{0},t)$ and spectra are summed, with $L(t_{0},t)\geq 1$ accounting for adiabatic losses (see Caprioli et al., 2010a; Morlino and Caprioli, 2012, for more details). Thus, we obtain a cumulative spectrum over the lifetime of the SNR. More specifically, $L(t_{0},t)$ can be written in terms of the time-dependent decompression of a fluid element with initial density $\rho(t_{0})\equiv\rho_{0}$. Since $\rho^{\gamma}(t)\propto v_{\rm sh}^{2}(t)$,

[TABLE]

where $4/3\lesssim\gamma\lesssim 5/3$ is the adiabatic index of the plasma and CRs (e.g., Caprioli et al., 2010a; Diesing and Caprioli, 2018).

*Magnetic Field Amplification—*The propagation of energetic particles ahead of the shock is expected to excite different flavors of streaming instability (Bell, 1978, 2004; Amato and Blasi, 2009), driving magnetic field amplification and enhancing CR diffusion (Caprioli and Spitkovsky, 2014c, a). The result is magnetic field perturbations with magnitudes that can exceed that of the ordered background magnetic field. This magnetic field amplification has been inferred via the X-ray emission of many young SNRs, which exhibit narrow X-ray rims due to synchrotron losses by relativistic electrons (e.g., Parizot et al., 2006; Bamba et al., 2005; Morlino et al., 2010; Ressler et al., 2014).

We model magnetic field amplification as in Caprioli et al. (2009); Caprioli (2012). Here, we assume that the pressure in Alfvén waves saturates at $P_{\rm w}\simeq P_{\rm cr}/(2M_{\rm A})$, where $P_{\rm w}$ and $P_{\rm cr}$ are the pressures in Alfvén waves and CRs normalized to the ram pressure $\rho_{\rm ISM}v_{\rm sh}^{2}$ and $M_{\rm A}\equiv v_{\rm sh}/v_{\rm A}=v_{\rm sh}\sqrt{4\pi\rho_{\rm ISM}}/B$ is the Alfvénic Mach number calculated in the amplified magnetic field. In the limit in which the fluid and Alfvénic Mach numbers $\gg 1$, we obtain

[TABLE]

where $u(x)$ is the fluid velocity normalized to $v_{\rm sh}$. Following the prescription described in Morlino and Caprioli (2012), we find an expression for the magnetic field in front of the shock,

[TABLE]

where the subscript 1 denotes quantities immediately in front of the shock. Behind the shock (denoted with subscript 2), the magnetic field strength is assumed to be $B_{2}\simeq\sqrt{1+2R_{\rm sub^{2}}/3}B_{1}$, since magnetic field components perpendicular to the shock normal are compressed. For $P_{\rm cr,1}\approx 10\%$, the shock parameters described above give $B_{2}$ near a few hundred $\mu$G, in good agreement with X-ray observations of young SNRs (Völk et al., 2005; Caprioli et al., 2008).

*Electron Spectrum—*Once the instantaneous proton spectrum, $f_{\rm p}(p)$, has been calculated, the instantaneous electron spectrum, $f_{\rm e}(p)$ is calculated as in Morlino et al. (2009) using the analytical approximation provided by Zirakashvili and Aharonian (2007):

[TABLE]

where $p_{\rm e,max}$ is the maximum electron momentum determined by equating the acceleration and synchrotron loss timescales. $K_{\rm ep}$ is the normalization of the electron spectrum relative to that of protons; its value ranges between $10^{-2}$ and $10^{-4}$ (Völk et al., 2005; Park et al., 2015; Sarbadhicary et al., 2017) but has no bearing on the spectrum slope.

To determine the cumulative spectrum over the lifetime of an SNR, the electron energy $E$ is evolved by integrating

[TABLE]

where the first and second terms account for synchrotron and adiabatic losses respectively (inverse Compton losses are subdominant). The weighted instantaneous spectra are then summed to determine a cumulative spectrum, an example of which shown in Figure 1. It is worth noting that electron escape upstream will have a negligible impact on this result, as escape is only important at the highest energies, where the acceleration time becomes comparable with the age of the system (Caprioli et al., 2010a). When losses are important, the diffusion length of electrons will not allow them to escape.

Results.— Having calculated the cumulative proton and electron spectra $\phi(E)$ at the end of the SNR lifetime, i.e., when the shock becomes subsonic and the remnant merges with the ISM, we estimate the power law slope as

[TABLE]

where $q$ is averaged between $10-10^{4}$ GeV for protons and between $10-100$ GeV for electrons. The energy ranges are chosen to ensure that particles are fully relativistic and to exclude high-energy cut-offs. The uncertainty in $q$ is estimated as the standard deviation over the range of calculation.

Figure 2 shows the resulting average slopes of electron and proton spectra as a function of $\eta$ for multiple values of $n_{\rm ISM}$. We find that the electron spectrum is consistently steeper than that of protons, regardless of acceleration efficiency or ISM density. For typical parameters ($n=1\ \rm{cm^{-3}}$ and $\eta\simeq 3\times 10^{-3}$, which returns the canonical $P_{\rm cr}\simeq 10\%$ (Caprioli and Spitkovsky, 2014b)), $\Delta q\equiv q_{\rm{e}}-q_{\rm{p}}\simeq 0.4$.

Figure 2 also implies that $q_{\rm p}$, $q_{\rm e}$, and $\Delta q$ depend on $\eta$; when the number of particles injected into DSA increases, proton and electron spectra steepen. More specifically, in the low $\eta$ limit, we recover the “test-particle scenario,” in which the CR pressure is small and magnetic field amplification inefficient. The result is a proton slope consistent with the standard DSA prediction ($q_{\rm p}\simeq 2$). As $\eta$ increases, so too does the efficiency of magnetic field amplification and thus the velocity of magnetic perturbations responsible for scattering CRs. Since Alfvén waves generated by CRs tend to travel against the fluid, this increase in magnetic field corresponds to a decrease in the effective compression ratio felt by CRs, resulting in a steepening of their spectrum (see Zirakashvili and Ptuskin, 2008; Caprioli, 2011, 2012).

An increase in $\eta$ also increases $\Delta q$, since larger $P_{\rm cr}$ and hence larger $P_{\rm B}$ lead to more severe synchrotron losses. Increasing $n_{\rm ISM}$ has a similar effect; the fraction of the bulk momentum flux converted to magnetic pressure is roughly constant at the few percent level such that $B_{2}\propto\rho_{\rm ISM}^{1/2}v_{\rm sh}$ (see Figure 4 for a clear illustration of this effect).

Figure 3 provides a more detailed picture of our calculated spectra and further illustrates the impact of synchrotron losses. The color scale indicates the magnitude of the instantaneous proton flux (top panel) and electron flux (bottom panel), weighted to account for energy losses, as a function of energy (x-axis) and time (y-axis). Note that the electron fluxes are multiplied by a normalization constant $1/K_{\rm ep}$ for display purposes. As suggested in the preceding section, the largest contribution to the proton spectrum occurs near the onset of the Sedov stage ($t\sim 200$ yr).

To understand how Figure 3 characterizes the synchrotron losses of CR electrons, recall that, in the case of protons, $E_{\rm max}$ is determined by equating the acceleration timescale, $\tau_{\rm DSA}$, with the lifetime of the remnant, $\tau_{\rm SNR}$, giving $E_{\rm max}\propto B_{2}(t)t^{-1/5}$ during the Sedov stage. This effect is apparent in Figure 3 (top panel), which shows a clear cutoff at high energies that decreases with time. Electrons, on the other hand, experience synchrotron losses on a timescale $\tau_{\rm syn}<\tau_{\rm SNR}$, meaning that their $E_{\rm max}$ is set by $\tau_{\rm DSA}\simeq\tau_{\rm syn}$. Moreover, steepening (or rollover) of the electron spectrum will occur at an even lower energy, $E_{\rm roll}$, above which $\tau_{\rm syn}\lesssim\tau_{\rm SNR}$. Since $\tau_{\rm syn}\propto B_{2}^{-2}E^{-1}$, we find that $E_{\rm roll}\propto B_{2}^{-2}t^{-1}$. Assuming that $B_{2}\propto v_{\rm sh}\propto t^{-3/5}$ during the Sedov stage, the result is that $E_{\rm roll}\propto t^{1/5}$. Again, this effect is apparent in Figure 3 (bottom panel), which exhibits a sharp steepening in energy, the position of which increases with time.

Since the difference between the CR electron and positron spectra arises from synchrotron losses, it is best understood in terms of the post-shock magnetic field, $B_{2}$ (see Figure 4). Namely, as $B_{2}$ increases, so does $\Delta q$, since an increase in magnetic field strength leads to more severe synchrotron losses, and thus a steepening of the electron slope. Most notably, the post-shock magnetic field strengths inferred, e.g., for the Tycho SNR correspond to $0.3\lesssim\Delta q\lesssim 0.4$ (gray band in Figure 4), in perfect agreement with its multi-wavelength emission (Morlino and Caprioli, 2012). In general, the typical magnetic fields that we estimate with CRAFT are consistent with those inferred from observations, implying that our conclusion that $\Delta q\gtrsim 0.3$ does not depend on modeling details of magnetic field amplification.

Discussion.— In summary, we used a semi-analytic model of non-linear diffusive shock acceleration to model the spectra of CR protons and electrons accelerated by SNRs. We find that electrons are injected into the Galaxy with spectra that are consistently steeper than those of protons, with a difference in slope of $\Delta q\simeq 0.1-0.4$. This steepening is the result of synchrotron losses in the large magnetic fields inferred in SNRs; therefore, it does not depend on the microphysics embedded in our model.

Our result may have significant implications for models of CR propagation in the Galaxy, which typically assume that protons and electrons are injected with the same spectrum. In particular, it must be reckoned with the “positron excess” reported by PAMELA (Adriani, [PAMELA Collaboration], 2013) and AMS-02 (Accardo et al. [AMS Collaboration], 2014). Comparisons between recent AMS-02 positron and electron data (Aguilar et al. [AMS Collaboration], 2019a, b), and DAMPE and CALET electron+positron data (Ambrosi et al. [DAMPE Collaboration], 2017; Adriani et al. [CALET Collaboration], 2017) suggest that the positron fraction increases with energy as $\chi\propto E^{0.3}$ between 10 and $\sim 300$ GeV. In the standard propagation paradigm, $\chi(E)\propto E^{-(q_{\rm p}-q_{\rm e}+\delta)}\propto E^{\Delta q-\delta}$, hence $\Delta q\gtrsim\delta+0.3$ could reproduce the rising in the positron fraction without introducing any source of primary positrons, be it astrophysical (e.g., pulsars) or exotic (dark matter). Although measurements of CR lithium, beryllium, and boron suggest that $\delta\sim 0.2-0.4$ (e.g., Aguilar et al. [AMS Collaboration], 2018), the antiproton to proton ratio is instead consistent with $\bar{\delta}\lesssim 0.1$ (Aguilar et al. [AMS Collaboration], 2016). This discrepancy is likely due to the increase in the proton-proton cross section with energy, which partially compensates for diffusive steepening with a hardening of the antiproton spectrum, parameterized by $\epsilon$ (Donato and Serpico, 2011; Korsmeier et al., 2018). Since positrons are produced in proton-proton interactions, too, taking $\bar{\delta}\sim\delta-\epsilon\lesssim 0.1$ also for the positrons implies that $\Delta q\simeq 0.4\gtrsim\bar{\delta}+0.3$ can entirely account for the “positron excess.” Note that, since the positron fraction is the ratio of lepton fluxes, radiative losses do not affect this conclusion. Moreover, the expected secondary production due to propagation saturates the normalization of the positron spectrum (e.g., Blum et al., 2013), while the positron to antiproton ratio is consistent with proton-proton branching ratios (e.g. Lipari, 2017, 2019; Blum et al., 2018). Intriguingly, these findings suggest that positrons may be of secondary origin after all, and that the “positron excess” may in fact be an electron deficit. This picture will be investigated more quantitatively in a forthcoming paper.