Nonthermal Radiation from Supernova Remnants

Nonthermal Radiation from Type Ia Supernova Remnants

Paul P. Edmon,, Hyesung Kang,, T. W. Jones,, Renyi Ma,,
Department of Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
Department of Earth Sciences, Pusan National University, Pusan 609-735, Korea
Department of Astronomy and Space Sciences, Chungnam National University, 305-764, Daejeon, Korea
E-mail: edmonpp@msi.umn.eduE-mail: twj@astro.umn.eduE-mail:
July 14, 2019

We present calculations of expected continuum emissions from Sedov-Taylor phase Type Ia supernova remnants (SNRs), using the energy spectra of cosmic ray (CR) electrons and protons from nonlinear diffusive shock acceleration (DSA) simulations. A new, general-purpose radiative process code, Cosmicp, was employed to calculate the radiation expected from CR electrons and protons and their secondary products. These radio, X-ray and gamma-ray emissions are generally consistent with current observations of Type Ia SNRs. The emissions from electrons in these models dominate the radio through X-ray bands. Decays of s from collisions mostly dominate the gamma-ray range, although for a hot, low density ISM case (), the pion decay contribution is reduced sufficiently to reveal the inverse Compton contribution to TeV gamma-rays. In addition, we present simple scalings for the contributing emission processes to allow a crude exploration of model parameter space, enabling these results to be used more broadly. We also discuss the radial surface brightness profiles expected for these model SNRs in the X-ray and gamma-ray bands.

radiation mechanisms: non-thermal, acceleration of particles, (ISM:) supernova remnants, (ISM:) cosmic rays
pagerange: Nonthermal Radiation from Type Ia Supernova RemnantsAcknowledgmentspubyear: 2010

1 Introduction

Diffusive Shock Acceleration (DSA) at astrophysical shocks has become the standard model for production of cosmic rays (CRs) (e.g.,  Malkov & Drury, 2001, and references therein). The CRs observed at the knee (eV) and below are most commonly explained through DSA in the blast waves of Galactic supernova remnants (SNRs) (Blandford & Ostriker, 1978; Lagage & Cesarsky, 1983; Blandford & Eichler, 1987; Drury et al., 2001). Several indirect arguments, including energetics and composition, support this explanation (e.g.,  Gaisser, 2005; Hillas, 2006; Ellison et al., 2007).

Direct evidence that at least some SNRs are capable of accelerating electrons and possibly hadrons to at least tens of TeV comes from observations of nonthermal X-ray filaments in some historical SNRs (e.g.,  Bamba et al., 2006; Parizot et al., 2006) and from TeV -rays observed in several shell SNRs (e.g.,  Funk et al., 2007; Lemoine-Goumard, 2007). The nonthermal X-rays are likely to be synchrotron emission from TeV electrons (e.g.,  Uchiyama et al., 2007). The TeV -ray origins, although they obviously require TeV or higher energy charged particles, have less certain origins and are sensitive to several model parameters. Inverse Compton scattering, secondary decays and bremsstrahlung may all contribute, but at levels that depend on details of the accelerated CR spatial and energy distributions as well as ambient photon, magnetic field and plasma properties. However, recent observations by Fermi of several Type II SNRs indicate that secondary decay from proton-proton interactions is the most likely explanation for the -ray emission there, which would observationally confirm the SNR origins of CR protons (Abdo et al., 2010).

In order to test the SNR origins of Galactic CRs, generally, and the accuracy and physical assumptions of DSA model calculations, specifically, it is important to determine from simulations not only the CR distributions, but also the related emissions, so they can be compared with current and anticipated observations (e.g.,  Blandford & Eichler, 1987; Hillas, 2005). Indeed several model SNR comparisons of this type have been published (e.g.,  Ellison & Cassam-Chenaï, 2005; Berezhko, Ksenofontov & Völk, 2009; Morlino, Amato & Blasi, 2009; Zirakashvili & Aharonian, 2010). Those demonstrate that observed emissions from some individual SNRs can be explained by the simulations with reasonable model choices, although they do not yet provide unambiguous model tests nor clear confirmation of the specific origins of TeV emissions.

At this early stage of our understanding it is especially useful to explore generally the nonthermal emissions produced by CRs in common classes of SNR models and how they depend on model parameters. In this spirit we present here an examination of the nonthermal radio to -ray emissions for Sedov-Taylor phase Type Ia SNRs, using a set of time dependent simulations that followed the evolution of spherical blasts, including nonlinear DSA at quasi-parallel shocks with thermal leakage injection of CR protons and electrons out of the bulk plasma. The pressure of CR protons is included in the gasdynamic equations, while CR electrons are treated as a test-particle component. The principal differences expected in the electron and proton CR populations come from differences in injection rates from the thermal population and from significant synchrotron/Compton radiative losses experienced by the highest energy electrons, but not the protons. We consider four different SNR models, including a range of uniform external environments and blast energies.

The structure of our paper is as follows. Section 2 provides a brief outline of the simulations that include both CR protons and electrons. A description of the code used to calculate nonthermal emissions as well as results from our emission calculations are presented in §3, while §4 summarizes our conclusions.

2 Nonlinear DSA simulations of SN Ia Sedov-Taylor remnants

2.1 The Spherical CRASH code

In the simulations discussed here the evolution of the CR modified SNR shock is followed using a one-dimensional, spherically symmetric version of the Cosmic RAy SHock (CRASH) gasdynamic-CR code. The basic gasdynamic equations modified to implement nonlinear DSA and full code details can be found in Kang & Jones (2006). CRASH applies adaptive mesh refinement (AMR) techniques and subgrid shock tracking to obtain high spatial resolution close to the blast shock, where it is crucial for converged solution of the coupled gas dynamics and CR transport equations. The spherical CRASH code also incorporates comoving coordinates expanding with the SNR blast. The momentum-dependent CR distribution is evolved through a time-dependent solution of the diffusion convection equation (DCE). CRASH incorporates dynamical backreaction from CRs onto the gas dynamics, including influence of the CR pressure, , and energy exchanges coming from injection of low energy CRs at the gas subshock and also dissipation of Alfvén waves stimulated by streaming CRs. We note that Berezhko and collaborators have extensively studied the similar problem, implementing a different, unconventional numerical method for the gasdynamics that normalizes the spatial variable by diffusion length at each momentum value and solves the CR kinetic equation iteratively to match the downstream and upstream solutions at the subshock (e.g.,  Berezhko et al., 1994; Berezhko & Völk, 1997; Berezhko, Ksenofontov & Völk, 2002). We find that the results of our simulations using the more conventional CRASH code are quite consistent with these previous studies.

Since DSA operates on relativistic electrons of a given in exactly the same way as it does on protons (i.e., it is exclusively rigidity dependent), the pitch-angle-averaged phase space distribution functions for CR proton and electron components, and both obey the DCE,


Here represents either or , and is the spatial diffusion coefficient in the radial direction (Skilling, 1975). Henceforth, particle momenta, , of both protons and electrons are expressed in units of , with the proton mass.

The cooling term is for protons and

for electrons combining synchrotron and IC cooling, where is the Thomson cross section, and is the effective magnetic field strength including the energy density of the ambient radiation field. We discuss the magnetic and radiation energy densities in §3.2. The code does not include direct backreaction from large scale magnetic fields, nor the pondermotive force of the wave turbulence.

The velocity represents the effective radial motion of scattering centers with respect to the bulk flow velocity, . Assuming that waves upstream of the gas subshock are predominantly resonantly generated through CR streaming and that the shock has a parallel magnetic field geometry, upstream of the gas subshock is set to the Alfvén speed, . The current version of CRASH does not follow the self-consistent evolution of the magnetic field strength through wave-particle interactions. So here we simply set , where is the upstream mean magnetic field strength. Downstream, the simulations assume , since the Alfvénic turbulence in that region is probably relatively balanced. The transition in scattering center motion, , reduces the effective velocity difference experienced by CRs across the shock compared to the bulk flow, , thus reducing the DSA efficiency. Gas heating due to Alfvén wave dissipation is represented by the term


where is a commonly used dimensionless parameter that controls the degree of dissipation. This dissipation term derives from a simple model in which Alfvén waves are resonantly amplified by streaming CRs in balance with local wave dissipation processes (e.g.,  Jones, 1993). As previously shown in SNR simulations (e.g.,  Berezhko & Völk, 1997; Kang & Jones, 2006; Caprioli et al., 2011; Ptuskin et al., 2010), and in more general contexts in earlier work (e.g.,  McKenzie & Völk, 1982; Markiewicz, Drury, & Völk, 1992; Jones, 1993), Alfvénic drift and precursor heating by wave dissipation both reduce DSA efficiency and associated shock modification. The significance of Alfv́en wave dissipation and drift can generally be characterized in terms of the Alfv́enic Mach number, i.e., , where is the upstream plasma flow speed relative to the shock. In these simulations , so the influence is modest.

Bohm-like spatial CR diffusion in the radial direction was used for both protons and electrons in the simulations; namely,


where , is the upstream magnetic field strength in micro-Gauss, and is the upstream density. The density dependence in this diffusion model represents enhancement of resonant Alfvén wave amplitudes through compression.

Low energy CR protons and electrons were injected at the shock in the same manner through thermal leakage (Gieseler, Jones, & Kang, 2000; Kang, Jones, & Gieseler, 2002). In this model, thermal ions in a Maxwellian distribution immediately downstream of the gas shock have a finite probability to escape upstream into the low energy CR population, provided they have sufficient upstream-directed momentum to overcome transverse MHD waves generated through the cyclotron instability in the shock transition (Malkov & Vök, 1998). This behavior is modeled numerically through a “transparency function”. The function contains one adjustable parameter, , which compares the upstream magnetic field, , with the amplitude of postshock, amplified MHD waves that interact with low energy particles, . The value was used for the models considered here. The resulting fraction of thermal protons injected into the CR pool, . It was shown previously that the DSA efficiency saturates at such a high injection rate (if for strong shocks), so except for details near start up, the SNR simulation results are insensitive to . The temperature used for the thermal population was computed from the gas pressure, , and density, , as with , so assuming equal proton and electron temperatures.

Because postshock thermal electrons have smaller gyroradii, compared to thermal protons, the injection rate of electrons should be much smaller. The ratio of CR electron to proton number injected at the shock, , is commonly assumed, since about 1% of the observed Galactic CR flux near a GeV is due to electrons (Schlickeiser, 2002; Reynolds, 2008). On the other hand, Morlino, Amato & Blasi (2009) point out that the CR electron flux observed here at Earth is actually a convolution of CR electrons from SNR and other electron accelerators of all ages. It is also entirely possible that CR electron injection, which is not well understood, and acceleration efficiency change dramatically for SNRs of different ages. On the other hand, it is commonly found in models of young SNRs that the best fit to the observational data requires (e.g.,  Berezhko, Ksenofontov & Völk, 2009; Morlino, Amato & Blasi, 2009; Zirakashvili & Aharonian, 2010), so we set in our calculations. However, since the electron population is passive and the SNRs are optically thin in emission bands of relevance here, our resulting emissions can be generally linearly scaled to another preferred value of .

In flows where CR backreaction is important, the formation of a CR precursor compresses and heats the inflowing plasma. These developments can lead to substantial changes in the strength of the dissipative, gas subshock and to the postshock conditions relative to those in a pure gas dynamic shock of the same Mach number. Given these flow modifications in front of the subshock, it is helpful in the following discussion to identify specifically the unmodified, upstream conditions by the subscript ’0’, the conditions immediately upstream of the gas subshock by the subscript ’1’, and the conditions immediately downstream of the full shock structure by the subscript ’2’.

2.2 Remnant model parameters

() (K) () () (pc) (years) ( km/s)

Note: The model ISM Alfvén speed, for S1,S2; for S3; for S4.

Table 1: SNR model parameters

Table 1 lists the dynamical parameters for the SNR models: the uniform ambient density, , the ISM temperature, , the upstream magnetic field strength, , and the SN explosion energy . The parameters for S1- S3 (S4) represent warm (hot) phase ISM environments. Models S1, S3, and S4 assume an explosion energy, erg, while model S2 begins with a blast containing four times this energy, so erg. All models assume there are no pre-existing CRs in the ambient ISM.

The SN ejecta mass is set to . The simulations are intended to follow evolution during the adiabatic, Sedov-Taylor (ST) evolution stage. Thus, convenient normalization variables are for density, along with ST similarity scales for length, and for time. The velocity and pressure scales are defined as and , respectively. In order to avoid confusion later, we note here a distinction between ST normalization subscripts, ’’ and initial or upstream condition subscripts ’’.

It is worth noting that the mass swept up by the blast at time and radius is given in the ST solution as


where is the ST similarity constant for a blast in a uniform medium with the gas adiabatic index . The ST shock speed can be expressed as


The second to last relation reveals that , emphasizing that the dynamical state of a ST blast is determined by the blast energy and the mass contained within the blast. Thus, except for numerical factors of order unity, many of the results we present below can be approximately represented independent of the ambient radial density distribution in terms of in combination with and rather than and .

Recent X-ray observations of young SNRs reveal in many cases magnetic fields stronger by at least an order of magnitude than the average ISM field (e.g.,  Bamba et al., 2006; Parizot et al., 2006). The existence of such strong fields is now commonly interpreted as the result of amplification within the shock precursor, either through resonant (Lucek & Bell, 2000), or nonresonant (Bell, 2004) streaming instability, or possibly hydrodynamic instability driven by the CR pressure gradient in the precursor (Beresnyak, Jones, & Lazarian, 2009). The simulations discussed here do not follow self-consistent amplification in the precursor of the magnetic field strength through these processes. Instead, to provide field values in the model SNRs consistent with observations, the upstream magnetic field strengths in SNR models S1 and S2 were set to the relatively large value, G. For comparison, the S3 and S4 models adopted, , which is similar to the mean ISM magnetic field.

In the simulation of the electron synchrotron energy losses and the resulting synchrotron emission, we assume that the field strength is modified during passage through the shock and into the SNR interior in a way that maintains consistency with the diffusion coefficient model (equation [4]). This model assumes in the relativistic limit that the scattering length, , along with a fixed ratio of wave field to total field strength. Consequently, it assumes , so that gas compression through the modified shock, (see Fig. 2), leads to postshock field values in these models, G in S1 and S2 and G in S3 and S4.

In all the simulations CR backreaction quickly increases compression from the initial to . Subsequently, decreases slowly in each case roughly as , with (see Fig. 2). Kang et al. (2009) found for strong, CR modified, plane shocks that approximately , where is the shock sonic Mach number. For these SNR simulations, , so their result would predict , which is reasonably consistent. For ST SNRs one has more generally that , so for scaling relations discussed in §3 we will assume that .

The SNR simulations start with ST similarity blast waves at a time , because the detailed evolution before that time does not strongly affect the later development. In fact, hydrodynamic simulations of SN blast waves (without the CR terms) show that the evolution of the outer shock speed can be approximated by for , although the true ST solution is established only after the inner reverse shock is reflected at the center at , (e.g.,  Kang, 2006). The simulations are carried out to , which is deemed sufficient to establish basic ST phase properties of the CR population and their dynamical impact until either the blast strength weakens substantially or the blast becomes radiative. Earlier studies showed that the highest momentum produced by DSA during the expansion of an SNR, , is achieved near the end of the free expansion stage, so when , and the transfer of explosion energy to the CR component occurs mostly during the early ST stage, so also when (e.g.,  Berezhko & Völk, 1997).

2.3 Evolution of the CR modified SNRs

Figure 1: The time evolution of gas density and temperature for model S1 along with CR pressure and the factor , which provides a scaling for thermal bremsstrahlung. Times shown are: (red dotted), 3 (green short dashed), 6 (blue long dashed), 10 (black solid), and 15 (magenta dot-dashed). (A colour version of this figure is available in the online journal.)

We now summarize the distinctive time evolutionary features of the SNR models used in this study, focusing on nonlinear DSA influences that impact radiative emissions. Fig. 1 shows radial profiles from model S1 of several useful quantities near the shock at times . The S1 model is representative in this regard of all the simulations. Gas density and temperature are shown in the top two panels. They are combined in the lower right to show as a proxy for bolometric thermal bremsstrahlung emissivity. The CR pressure profile is included in the lower left panel. To simplify the plots, initial ST conditions from are not included. At the first time shown in Fig. 1, , the shock is modified substantially through compression and heating in the CR precursor. In fact, the CR acceleration and nonlinear modification peak early in the evolution and then decrease in time as the shock slows down (see also Fig. 2).

CR modification of the shock enhanced the postshock density but reduced the gas pressure, leading to a postshock gas temperature lower than that expected for the ST shock near . The effect of CR pressure on the temperature distribution behind the shock can be seen by noting that the slope of temperature rise towards the interior steepens as begins to drop sharply inside the blast. This interior drop in is an artifact of the location of the contact discontinuity separating matter interior and exterior to the initial conditions shock. Its shape depends on the rate of spatial CR diffusion into the blast cavity.

In addition to raising the total density jump within the shock, compression through the precursor also preheats gas before it enters the subshock structure. This adiabatic effect by itself would raise the upstream temperature by a factor , or about 1.6 in the S1 model. In fact, Fig. 1 shows that in the S1 model . The additional heating comes from MHD wave dissipation (equation [3]). Compression through the precursor also decelerates upstream flow before it enters the subshock. All these effects greatly weaken the subshock. So the sonic Mach number of the subshock decreases quickly to in the S1 and S4 models and to in S3 model by and remains approximately the same afterward, even though the total shock sonic Mach number decreases roughly as . Since stays constant throughout the evolution, the CR injection fraction () via thermal leakage also stays constant after (see Fig. 2). With the adopted injection parameter (), the CR injection fraction is high, , and so the CR acceleration should be close to maximum efficiency.

These shock modifications also significantly influence the thermal bremsstrahlung emission in the SNR, whose distribution is represented by the profile in the lower right panel of Fig. 1. The thermal bremsstrahlung shell will appear thinner in the strongly modified SNRs than in the analogous gas dynamic SNR. The enhanced postshock density also enhances the thermal bremsstrahlung emissivity. The thinner shell and higher emissivity in the shell combine for a relatively small difference in the bolometric thermal bremsstrahlung luminosity compared to an unmodified shock. The thermal bremsstrahlung luminosity in a fixed spectral band below the spectral cutoff should generally be enhanced, however, since that has an approximate dependence with respect to the bolometric emission. Of course, the lowered temperatures behind the modified shock also move the upper cutoff of the thermal bremsstrahlung SED from the hard X-ray band into the soft X-ray band.

Figure 2: The time evolution of the compression factors, and , postshock pressures, and , the injection fraction, , and the volume integrated energy, and is shown for six models. The wave heating parameter is for S1r, S3r and S4r models, while for S1, S3, and S4 models. The r models are the dashed lines. (A colour version of this figure is available in the online journal.)

More insights to the evolutionary behaviors of the S1, S3 and S4 models can be obtained from Fig. 2. As already noted, both the precursor compression, , and the total shock compression, , as well as the CR pressure relative to increase quickly from their initial values 111 The very early evolution is affected by numerical start-up of the simulations., reach their maximum values around , and decrease slowly in time as the shock slows down. In case of the S2 model with higher (not shown in Fig. 2), the CR acceleration and the nonlinear modification are slightly more significant compared to the S1 model, because of a faster shock speed. If the Alfvénic drift and heating were not included, the temperature of the ambient medium (and so the sonic Mach number) would be the primary parameter that determines the CR injection and acceleration efficiencies. Compare, for instance, the S3 model in warm ISM to the S4 model in the hot ISM. Because of the ISM temperature differences, the sonic Mach number of the total shock in the S3 model is times higher than S4 at a given . So the S3 SNR would be the more efficient CR accelerator. The low density and high temperature ISM of the S4 model produce low sonic and Alfvénic Mach numbers (, ), making the CR acceleration much less efficient in that model.

The magnetic field strength also has a significant influence on the SNR evolution even though no direct MHD dynamical effects are modeled. Most important to these simulations, the field strength sets the CR spatial diffusion rate and thus, in concert with , the CR acceleration rate (see §2.4, below). Consequently, particles are accelerated to several times higher momenta () in the S1 model in comparison to the S3 model because of the enhanced magnetic field and smaller diffusion coefficient (see Fig. 3 below). On the other hand, the stronger field leads to a faster Alfvén speed, which tends to reduce DSA efficiency by adding entropy to the precursor. This reduces the effective velocity jump across the subshock and also reduces the subshock Mach number. Those effects are responsible for the moderate reductions in CR acceleration efficiency in the S1 model compared to the S3 model. The effect of Alfvén wave heating is further illustrated in Fig. 2, where we show three additional simulations, S1r, S3r, and S4r, in which the wave dissipation parameter is reduced from to . The smaller wave dissipation rate reduces non-adiabatic heating in the precursor, allowing greater compression through the precursor, so an increase in CR acceleration efficiency.

Fig. 2 also shows the volume integrated thermal energy () and CR energy in units of the SN explosion energy. The energy transfer to CRs seems to saturate for , and the fraction of the blast energy transferred to CRs () is about 0.7, 0.5, 0.45, 0.35 in the S3, S2, S1, and S4 models, respectively.

In summary, adopting simple models for Alfvén wave transport, these simulations demonstrate that the CR acceleration would be less efficient at SNRs in hot rarefied ISM where both the sound speed and Alfvén speed are faster, compared to those in a warm phase ISM.

We note that the evolution of CR modified SNRs shown in Fig. 2 (especially, , , and ) differ somewhat from the similar SNR models presented in Kang (2006). Several parameters were different there, e.g., , , K for S1-3 models, and G for the S4 model. The main differences result, however, from a coding error in those earlier simulations, which caused the Alfvénic drift and heating effects to not be fully implemented. The results in the present study are, however, consistent with independent simulations of similar SNRs (Berezhko, Ksenofontov & Völk, 2002), in which a Type Ia SNR with erg and in warm phase ISM with , K, and G was calculated.

2.4 The CR proton and electron populations

Figure 3: Volume integrated proton spectra in terms of particle kinetic energy at (red dotted), (green short dashed), (blue long dashed), (magenta dot-dashed), and (black solid). (A colour version of this figure is available in the online journal.)

Fig. 3 shows the volume integrated CR proton energy spectrum, , where , is the proton kinetic energy, and . We note that in the models with enhanced (i.e., S1 and S2) particles of charge Z can be accelerated to eV by the early ST stage. This is consistent with results of other analogous SNR simulations (e.g.,  Berezhko & Völk, 2007). In the early ST stage the shock structure is significantly modified, with , so is concave upwards. This well-known behavior is a consequence of momentum dependent diffusion across the precursor; CRs near the injection momentum experience only the subshock velocity jump, while higher momentum CRs see a larger velocity jump. Thus, the CR spectrum is softer at lower momentum and harder just below the high energy cutoff than the strong shock test particle result, . Moreover, the Alfvénic drift in the precursor further softens the spectrum of newly accelerated CRs as the shock slows down and the Alfvénic Mach number decreases.

It will be useful below to estimate the maximum CR proton momentum expected at a given time . This we do by integrating the standard expression for the average momentum gained per unit time for a particle (Lagage & Cesarsky, 1983),

where, once again, subscripts and refer to upstream and downstream conditions in the shock frame, respectively. The useful parameter,


is the DSA test particle spectral index, including the upstream Alfvén wave drift term, , for consistency with equation (1). The right-most expression for neglects this term, since it is a relatively minor correction in our model SNRs.

Neglecting the Alfv́en wave term and assuming a constant compression ratio, , but including the diminishing shock speed in the ST solution, equation (2.4) gives for protons accelerated between and ,


For our simulations, which start from , equation (9) asymptotes to eV/c at large using for convenience. At this gives, for example, eV/c for the S1 model and eV/c for the S4 model. These estimates are consistent with the numerically determined cutoff energies in shown in Fig. 3.

Figure 4: Volume integrated electron spectra in terms of particle kinetic energy at (red dotted), (green short dashed), (blue long dashed), (magenta dot-dashed), and (black solid). (A colour version of this figure is available in the online journal.)

The CR electron spectra develop substantial differences from the proton spectra, despite their common DSA interactions (e.g.,  Webb, et al., 1984; Berezhko, Ksenofontov & Völk, 2009; Blasi, 2010; Zirakashvili & Aharonian, 2010). Fig. 4 shows the volume integrated CR electron energy spectra, , at several times for each SNR model. Note that the electron spectra cutoff at much lower energies than the analogous proton spectra shown in Fig. 3. The structures of the integrated electron spectra are complex around and below the cutoff. The low energy break in each spectrum identifies the energy at which the synchrotron/Compton loss time for electrons downstream of the shock equals the SNR age. Below this energy radiative energy losses are negligible. The integrated spectrum between this break and the high energy cutoff is controlled by the downward movement of the electron energy cutoff as recently accelerated electrons move into the SNR interior away from the shock. The high energy cutoff itself develops at the shock and corresponds to the energy where synchrotron losses balance DSA gains.

At the shock the electron cutoff momentum, , can be evaluated by balancing DSA and radiative losses. In particular, corresponds to the momentum at which the average momentum gain by an electron in one pair of shock crossings, , equals the momentum loss from radiation during the same period of time, , where is the electron velocity, with and defined by equation (2.1) in terms of the upstream and downstream effective magnetic fields, respectively (e.g.,  Webb, et al., 1984). This leads to the (test particle) result,


where is the instantaneous shock speed. Here we assume and the second expression corresponds to the limit, , which is applicable in our simulations. This translates into a cutoff energy for CR electrons at the shock,


where the subscript on B indicates field strength in microGauss.

Equation (11) evaluated with and predicts TeV for these models at time , which is consistent with the simulation results shown in Fig. 4. This is also consistent with the upper limits on electron spectra derived from observation for a number of observed SNRs (Reynolds & Keohane, 1999) and with SNR simulation results, for example, of Berezhko, Ksenofontov & Völk (2002). During the early ST stage, the maximum acceleration momentum gradually saturates as given in equation (9), but then it eventually asymptotes to determined by the instantaneous shock speed as given in equation (2.4). As a result, for larger values of , the electron cutoff energy is smaller, but the maximum proton energy is larger.

Electrons advected into the SNR interior will continue to lose radiative energy with the combined synchrotron/IC radiative cooling time, , given by


Setting the cooling time behind the shock to the SNR age, , provides a conservative, rough estimate for the minimum electron momentum influenced by radiative cooling,


Recall that . Below the integrated electron spectrum should be very similar in form to the proton spectrum at the same momentum, consistent with results shown in Fig. 3 and Fig. 4.

To help interpret the integrated electron spectra, , between the break momentum and the cutoff momentum , we note in the thin shell approximation that , and that the local electron spectrum, , below . In this case it is easy to show that the integrated electron energy spectrum below would be steeper by about one in the spectral index than the integrated proton spectrum in the same energy range. The steepened slope would extend down to the energy GeV. On the other hand, in spherical flows following roughly the ST behavior, the postshock flow accelerates away from the high magnetic field region just inside the shock while reducing the magnetic field strength. These effects reduce the rate of cooling from that implied in equation 13, thus shortening the interval between and . Then the spectral steepening of the integrated electron population compared to the proton spectrum will be greater than unity, as observed in our results.

3 Modeling of continuum emissions from the SNR models

3.1 The Cosmicp code for radiative processes

The nonthermal radio to -ray emissions expected from CR electrons and from proton secondary products in the simulated SNRs were computed through the post processing of model data using Cosmicp, a general purpose nonthermal continuum emissions code developed in-house around published radiative process and inelastic particle scattering formulations. Cosmicp computes direct emissions from input electron/positron populations and calculates the inelastic collision products for photon and proton interactions with matter, including pion and lepton secondaries. Lepton secondaries are incorporated into Cosmicp, but they do not make significant contributions to SNR emissions, so are not included in our analysis here. Cosmicp was designed to be very flexible, so makes as few assumptions as possible. For example, the input energy distributions of energetic particles are arbitrary. Using inputs for CR electrons and protons along with the relevant environmental information, such as magnetic field, radiation field, and gas density, Cosmicp calculates the volume emissivity per unit frequency, , for each radiative process, where is the emitted photon frequency. Bremsstrahlung from the thermal electron population is also included in our calculations with Cosmicp, but electronic emissions involving discrete atomic transitions are excluded. All calculations are in cgs units, which allows easy linking to actual observables.

Specific CR electronic emissions incorporated into our analysis here include synchrotron, inverse Compton (IC), and bremsstrahlung processes. The synchrotron, IC, and relativistic bremsstrahlung formalisms followed presentations in Schlickeiser (2002). The nonrelativistic bremsstrahlung formalism used to compute the spectrum due to the thermal population followed Jackson (1999). The IC formalism uses an arbitrary, user supplied spectral form for the ambient photon field, so Cosmicp can model interactions with multiple blackbody or non-blackbody sources. It does, however, currently assume an isotropic incident radiation field.

Hadronic interactions include, as noted, inelastic proton-proton and photon-proton interactions. Photopion production in the latter case can be calculated from an arbitrary photon field using the formalism in Kelner & Aharonian (2008). Pion generation from proton-proton interactions uses relations presented in Kelner, Aharonian, & Bugayov (2006). The principal radiative process that results from hadronic interactions in our SNR context is secondary decay, which is dominated by the channel. Cosmicp can include helium in the hadron-hadron production of photons. However, that correction was not applied in the present calculations. For a CR composition similar to that incident on the earth and an ISM with normal metalicity this correction would increase the pion production rate by % (Schlickeiser, 2002).

3.2 Nonthermal emissions: Local interaction conditions

Rates for nonthermal electron bremsstrahlung and p-p collisions leading to pion decay ’s depend on the product of the local plasma density and the CR electron and proton density, respectively. These are taken in our emission calculations directly from the simulations. Synchrotron and IC emissions depend on the spatial magnetic field distribution and the ambient radiation field respectively. As discussed in §2.2, we assume the local magnetic field strength scales with thermal plasma density; i.e., . The ambient photon field for our calculations is uniform and models the Galactic radiation field in the solar neighborhood as given by Schlickeiser (2002). It consists of a sum of four components: the Cosmic Microwave Background (CMB) blackbody at K with an energy density of , dust at K with an energy density of , old yellow stars at K with an energy density of , and young blue stars at K with an energy density of . The total energy density for the ambient photon field is . In evaluating equation (12) for electron energy losses this translates to an effective magnetic field of . We point out in §3.3.2, however, that IC scattering by the highest energy electrons of all but the CMB and dust emission radiation fields is limited by electron recoil, and reduced by the Klein-Nishina cutoff, so this value overestimates the effective . In our models electron losses are dominated by the synchrotron process, so this is a relatively small correction, in any case.

3.3 Volume integrated radiation spectrum from the model remnants

Figure 5: The volume integrated SED for the S1 model at : total emission (solid black line), synchrotron emission (red dotted line), IC emission (green short dashed line) bremsstrahlung (blue long dashed line), s from decays (magenta dot-dashed line). (A colour version of this figure is available in the online journal.)

A general impression of the radiative emissions produced in the simulated SN Ia remnants can be obtained from Fig. 5, which shows for the S1 simulation at the volume integrated Spectral Energy Distribution (SED), from radio frequencies to PeV -rays. Synchrotron emission dominates the spectrum for Hz222Useful translation factors are: 1 keV Hz, 0.41 keV.. There is a narrow -ray window, dominated by bremsstrahlung from CR electrons interacting with thermal plasma. Above (keV) the spectrum is dominated by photons from decays for the ISM density assumed in the simulation.

Figure 6: The total volume integrated SEDs for models S1-S4 over time. The dotted red line is , the short dashed green line is , the long dashed blue line is , the dot-dashed magenta line is , and the solid black line is . (A colour version of this figure is available in the online journal.)

The spatially integrated radiative emissions at selected times for all four models are shown in Fig. 6. Independent of time, the radiation is mostly dominated by the electronic CR emissions at lower energies, but, for the most part, proton decay at high energies. The additional bump in the UV and X-ray seen at late times especially in model S1 is due to thermal bremsstrahlung. This feature appears as the high energy end of the electron distribution is depleted by radiative losses, revealing thermal emission from hot, postshock plasma. The lower gas density in the S4 model reduces this thermal contribution, and also the nonthermal bremsstrahlung -rays.

In the S4 model the -ray portion of the spectrum is a composite, dominated by decays between Hz ( MeV - MeV) and again above about Hz ( TeV). However, the lower ambient gas density of the S4 model also reduces the production enough that IC emissions are dominant between the synchrotron cutoff and . IC emission is also predominant between Hz ( GeV - TeV).

We outline below some useful relations for interpretation of the calculated radiative emissions in the context of these models, including some simple scaling relations that can be helpful in extending our results to different model parameter choices than we have used. We begin with discussions of the electronic emissions and follow with the essential elements to understand the calculated emissions produced by inelastic p-p collisions.

3.3.1 Electron Synchrotron Emission

The electron synchrotron emission spans a frequency range from radio to roughly the emission peak for electrons near the Lorentz factor . Above this energy the electron population is heavily modified by losses and begins to drop off severely. The associated synchrotron peak frequency is , where


and is the nonrelativistic electron cyclotron frequency (Blumenthal & Gould, 1970). For the calculations presented here we set the angle of the magnetic field with respect to the observer, .

Using equations (2.4) for the upper electron energy cutoff and (14) we can derive the following estimate for the cutoff frequency in the synchrotron spectrum in the limit that ; namely,


For and the factor , so , consistent with the results in Fig. 6. As noted by previous authors (e.g.,  Bamba et al., 2003; Berezhko & Völk, 2004) the cutoff frequency is determined mainly the instantaneous shock speed and independent of magnetic field strength except through its implicit dependence on the factor . This explains the similarity between the values for the synchrotron cutoffs for S1 and S3 even though they have very different magnetic fields.

The break frequency that corresponds to the estimated break momentum depends only on magnetic field strength and the shock age as


The values of and generated from equation (15) and (16), respectively, are consistent with the SEDs shown in Figs. 5 and 6. The photon energies corresponding to the synchrotron cutoff frequency fall in the range keV, consistent with observed Type Ia remnants (Reynolds & Keohane, 1999).

The synchrotron luminosities of all four simulations below the X-ray cutoffs can be approximately related to each other by a simple scaling relation that, in addition, provides a means to extend these results approximately to different parameter choices. For the volume integrated electron energy spectrum shown in Fig. 4 can each be represented by


where represents the suprathermal injection energy of CR electrons, is the volume integrated number of CR electrons and is the energy distribution mean power law slope between and . For reference, recall that momentum and energy distribution power law indices are related as, . The synchrotron luminosity can then be expressed approximately as (Blumenthal & Gould, 1970)


Since and the injection faction is nearly a constant for in these simulations, . The injection energy, scales roughly with the shock speed, so also as . In §2.2, we established the approximate scaling, , so . Then in these simulations. For example, for or for . Since the index varies slowly within the range in these simulations, the radio luminosity should depend very weakly on or , consistent with the results in Fig. 6.

Note that equation (18) alone is not sufficient to estimate synchrotron emission behaviors within the X-ray band. The synchrotron emission in this band is generated predominantly by electrons near the spectral break and cutoff . Consequently, X-ray synchrotron emissions are heavily influenced by the positions of these frequencies.

3.3.2 Electron Bremsstrahlung and IC Emission

The bremsstrahlung and IC contributions are mostly subdominant in our models. At early times synchrotron emission dominates thermal bremsstrahlung over the bands where they both contribute. The bremsstrahlung peak is roughly , while the synchrotron cutoff, , is given by equation (15). Both cutoffs depend primarily on the shock speed squared, so evolve more or less together. The decreasing role of synchrotron emission in the X-ray band results, on the other hand, from the fact that CR electron energy losses downstream of the shock cause that population to become further depleted, dropping the relative emitted power between and , revealing the thermal bremsstrahlung emission. Consequently for in S1 thermal bremsstrahlung is an important X-ray contributor (see Fig. 6).

In another exception to the general emission behaviors, we note that IC emissions near TeV in model S4 exceed the otherwise dominant pion decay gamma-rays. The S4 decay emissions, on the other hand, are smaller by two orders of magnitude, reflecting the similarly lower S4 ambient density (see equation [24]).

Some comments on the contributions and form of the IC SED may be useful. As noted in §3.2 the IC emissions result from a combination of incident photon fields with blackbody forms at different characteristic temperatures. The effective incident field is just their sum. While the CMB is a blackbody in both spectral form and energy density, the radiation from dust and stars has greatly diluted energy density compared to the Planck functions for their respective temperatures. In each case we can assign a frequency-independent dilution factor, , given by the ratio of the local energy density to the appropriate Planck function. The dilution factors for the radiation field properties listed in §3.2 are: for the CMB, , for thermal dust emission, , yellow stars, , and for blue stars, .

The IC scattering source function in the limit of Thomson scattering of incident photons with a Planck spectrum with color temperature, , and dilution factor, , is


where is the classical electron radius, is a dimensionless scaling function as defined in Blumenthal & Gould (1970); for . We have used the form of the electron spectrum given in equation (17).

For a representative CR electron slope, , the different background radiation components contribute in equation (19) according to factors , which leads to a ranked list represented in the ratios . This suggests that the CMB and the dust IR radiation are the predominant contributors to the IC SED, despite the roughly comparable photon energy densities of the stellar components (see §3.2). In this regard it is also worth noticing that the IC spectral peak at a few TeV involves photons scattered to energies close to , the maximum possible, whereas the Thomson scattering limit would take photons into the PeV range when scattered from the incident stellar photon fields. In fact, photons from those higher temperature incident fields are scattered in the Klein-Nishina limit by the most energetic electrons, as given by the condition on the incident photon energy, . This condition can be expressed practically as


For only the CMB and, marginally, the dust radiation escape this limit. The scattering cross section in the Klein-Nishina limit is reduced from the Thomson cross section, so that, in fact, gamma ray IC emissions in these models are dominated by the CMB and dust emissions.

We can use equation (19) to construct an approximate, simple scaling relation for the IC luminosity, , in the Thomson regime assuming an electron distribution of the form in equation (17),


where the radiation energy density, , would be approximated by the sum of the CMB and dust radiation densities. The IC luminosity, , scales as , similar to the synchrotron luminosity . In fact, the ratio of the two for emissions at synchrotron frequency , and IC frequency, in, say, the radio and gamma ray bands, is just


This indicates that the ratio between the two contributions to the SED for each of our SNR models remains roughly constant over time. Since is model independent, the ratio of IC gamma ray luminosity to radio synchrotron luminosity varies in these models roughly as .

3.3.3 Decay Gamma-Ray Radiation

Gamma ray decay products from inelastic p-p collisions, and , in particular, are commonly viewed as the “smoking gun” for hadronic DSA in SNRs. As noted above these emissions mostly dominate the gamma ray SEDs in the SNR models being discussed here. Although the detailed modeling of decay gamma rays is fairly complicated, a simple analytic approximation offered by Pfrommer & Enßlin (2003) provides a good understanding of the gamma ray SED in these models. In particular, if the local CR proton energy spectrum above the threshold for pion production () is a power law, , then the -ray source function, , (photons/time/volume/energy) is


where is an effective p-p cross-section for pion production, , while is the density of the thermal protons, and is the mass.

This source function peaks for MeV. At low energies it asymptotes to a power law , while it asymptotes at high energies to a power-law . The source function cuts off above roughly 10% the proton high energy cutoff (Kelner et al. 2006). The radiated power per unit volume (the SED) will scale with . Using equation (23) for the SED would peak around . There is no SED maximum for , except at the high energy cutoff.

Assuming a proton spectrum of the form in equation (17) one obtains a simple high energy scaling relation for the volume integrated pion decay luminosity


This scaling is consistent with the results shown in Fig. 6, showing in particular why the pion decay emissions are similar in all the models except for S4 because of its lower .

The proton distributions in these SNR models are not really true power laws, of course, but have concave upwards forms; that is, decreases with energy. Protons near the pion production thresholds just above a GeV are mostly responsible for gamma rays near the low energy gamma ray peak. In that range typically (see Figs. 3 and 4), The expected gamma ray SED peak is near 1 GeV (Hz), consistent with Figs. 5 and 6. The model SEDs in Figs. 5 and 6 also exhibit a second, higher energy peak. That results from the concave upwards character of the proton spectra and their eventual cutoffs around or above a PeV for . Except for model S4 the proton spectra have below their cutoffs, making their upper gamma ray peaks quite prominent.

Since the the upper gamma ray peak is so prominent it can be used to determine the peak proton energy in the SNR (e.g.,  Aharonian et al., 2007). From the gamma-ray production relations in Kelner et al. (2006) one can see that the peak in gamma-ray production for a monoenergetic proton of occurs at a photon energy . Looking at Fig. 6 for S1-S3 the gamma-ray peak energy is around about TeV ( Hz) which corresponds to a proton energy of TeV. This is consistent with the peaks in the proton spectra in Fig. 3.

Both IC and source functions depend on the number of CRs accelerated, while the pion-generated gamma rays also depend on the density of the thermal plasma. The ratio of found from equations (23) and (19),


provides a useful comparison of the two processes for estimation and scaling purposes.

Again, taking as representative, then including CMB and dust contributions in the radiation field as discussed earlier, and setting along with , appropriate to postshock conditions for model S1 (see Table 1 and Fig. 2), we can use equation (25) to estimate


which is valid roughly for MeV (Hz). This number is in reasonable agreement with the numerical results displayed in Fig. 5. Equation (26) suggests that decay emissions would need to be reduced by at least two orders of magnitude from those in model S1 before IC emissions would become dominant in the GeV band. Note for scaling purposes from equation (25) that . Consequently, GeV IC emissions would dominate when , consistent with the model S4 results as displayed in Fig. 6.

3.4 X-ray and gamma-ray brightness distributions