CORE-COLLAPSE MODEL OF BROADBAND EMISSION FROM SNR RX J1713.7−3946 WITH THERMAL X-RAYS AND GAMMA RAYS FROM ESCAPING COSMIC RAYS

The direct observation of radio, non-thermal X-ray, and GeV–TeV emission from several young supernova remnants (SNRs) provides unambiguous proof that these objects produce relativistic particles. The morphology of this emission, which is often seen in thin, rim-like structures associated with the supernova (SN) blast wave, provides convincing evidence that diffusive shock acceleration (DSA) is the mechanism producing these particles. This and other evidence, such as the direct measure of efficiency at the Earth bow shock (e.g., Ellison et al. 1990), and the likely existence of amplified magnetic fields in some SNRs, suggests that the efficiency for producing superthermal ions, i.e., cosmic rays (CRs), is high in DSA. A critical issue for DSA, and for the origin of CRs, however, concerns the radiation mechanism responsible for the GeV–TeV emission. Is it dominated by pion-decay emission from ions or inverse-Compton (IC) emission from leptons?

Given relativistic electrons and ions, three mechanisms—non-thermal bremsstrahlung, IC scattering, and pion decay—can produce GeV–TeV photons. Sorting out this mix in a particular remnant not only provides information on how CRs are produced, it constrains the basic physics of DSA, in particular the electron to proton ratio, K_ep, of accelerated particles. While most modelers of young remnants claim that shock-accelerated protons are mainly responsible for the GeV–TeV emission (e.g., Morlino et al. 2009; Berezhko & Völk 2010), some researchers (e.g., Porter et al. 2006; Katz & Waxman 2008; Zirakashvili & Aharonian 2010) have suggested that IC from relativistic electrons dominates the GeV–TeV emission from SNR RX J1713.7−3946 (henceforth SNR J1713). Recently we have definitively shown, with a model that includes the production of thermal X-ray lines self-consistently with the broadband continuum emission (i.e., Ellison et al. 2010), that IC emission from relativistic electrons clearly dominates pion decay if SNR J1713 is in a homogeneous environment typical of a Type Ia SN.

The assumption of a homogeneous environment, although typically assumed in self-consistent models of SNR J1713, is an important caveat because SNR J1713 is more likely a core-collapse SN that may have exploded in a complex environment (e.g., Fukui et al. 2003; Moriguchi et al. 2005; Inoue et al. 2009; Sano et al. 2010). Because of this, we have generalized our model in two important ways. First, we can now model a spherically symmetric, non-homogeneous environment. We consider a SN that explodes in a pre-SN wind surrounded by a dense, swept-up shell of wind material. We further allow for external dense material that might be associated with molecular clouds out of which massive stars form. Second, in addition to γ-rays produced by CRs trapped within the forward shock (FS), we now include γ-ray production from the CR protons that escape from the FS and diffuse rapidly into the dense shell and/or the molecular cloud material beyond the FS.

These generalizations allow a fairly realistic model of a core-collapse SNR J1713 interacting with dense material, at the expense of adding model parameters. However, despite the added parameters and the flexibility in fitting they allow, we find it is impossible to produce a good fit to the broadband data, including the X-ray emission lines, with pion decay dominating IC at GeV–TeV energies unless the SNR is interacting with a dense mass concentration (≳ 10⁴ M_☉) such as a molecular cloud (e.g., Gabici & Aharonian 2007; Gabici et al. 2009). An isolated SNR, regardless of its type, even allowing for escaping CRs interacting with a dense shell of material with M_shell ≲ 100 M_☉ from a pre-SN wind, is excluded.

The conclusion of Ellison et al. (2010), that IC from electrons dominated the GeV–TeV emission, hinged on the assumption that both non-thermal pion decay and thermal X-ray line emission scale in approximately the same way with the ambient density. This is expected to be the case if the shock-accelerated particles are drawn from the target material, i.e., if the material swept up by the FS provides seed particles for DSA and simultaneously provides a dense target for proton–proton collisions. If the environment is homogeneous and has a uniform density everywhere outside of the FS, as assumed in Ellison et al. (2010), photons produced outside of the FS, where the density is considerably less than behind the shock, can be ignored. In this case, the parameters chosen for the injection of electrons and protons into DSA, along with the ambient density and magnetic field, largely determine the IC/pion-decay ratio and the thermal X-ray emission/pion-decay ratio. As shown by Ellison et al. (2010), reasonable parameters for SNR J1713 exclude pion decay as the dominant radiation mechanism for the GeV–TeV emission in this scenario.

There are three main populations of shock-accelerated CRs that are important for producing γ-rays: relativistic electrons producing γ-rays through IC and non-thermal bremsstrahlung, CR ions that remain trapped within the FS, and CR ions that are accelerated by the FS but escape upstream.⁴ These three populations are produced simultaneously by DSA in the blast wave and, while they have different properties, their radiation signatures, particularly IC and pion decay, are not easily distinguished at GeV–TeV energies with current observations. The recent Fermi Large Area Telescope observations of SNR J1713, however, now clearly support a leptonic model for the GeV–TeV emission (Abdo et al. 2011).

An important distinction between trapped and escaping CR ions is that, in a non-homogeneous environment, where regions outside of the FS can have an elevated density, the escaping CR ions can diffuse into these dense regions and produce γ-rays without producing any corresponding thermal X-ray or IC emission.⁵ The production of these external γ-rays depends mainly on the amount of external mass but the diffusion properties of the CRs in the circumstellar medium (CSM) also come into play. In Ellison & Bykov (2011) we presented a simple Monte Carlo model for escaping CR diffusion and investigated how the γ-ray flux depends on the CR diffusion without application to a particular SNR. For reasonable values of the diffusion, it is clear that, if the escaping CRs have enough external mass to interact with, the γ-ray flux from pion decay can overwhelm that from IC (e.g., Gabici & Aharonian 2007; Lee et al. 2008). Other observational considerations, such as the existence of nearby molecular clouds and the spatial coincidence of radio, X-ray synchrotron, and γ-ray emission must be used to discriminate between a hadronic or leptonic origin in this case.

Here, we model the broadband emission from SNR J1713 for four non-homogenous environments. The first is for an isolated core-collapse SNR expanding in a slow pre-SN wind with an external, dense shell of swept-up wind material. For an isolated SNR, the mass of the swept-up wind must be less than the mass of the main-sequence star and we use M_shell = 100 M_☉ as an extreme case. In the second example, we increase the external mass to M_shell = 10⁴ M_☉ to mimic a SNR which is near a mass concentration but has not yet impacted it. In the third example, we allow the SNR to impact a dense shell of 10⁴ M_☉ and follow the evolution of the remnant and the photon production as the FS moves into the shell. Finally, we consider the case of a fast pre-SN wind.

We find, largely independent of diffusion parameters, that escaping CRs interacting with an external mass shell from a pre-SN wind with mass ≲  100 M_☉ will not produce enough pion-decay emission to compete with the IC emission from trapped CR electrons. If larger external masses ≳  10⁴ M_☉ are allowed, as would be the case for a SNR near or interacting with a molecular cloud, the pion decay can dominate the IC at TeV energies. Now, however, other considerations come into play such as the shape of the escaping CR distribution and the details of how CRs can stream into a molecular cloud and produce self-generated turbulence. We expect the escaping CR distribution to be far too narrow to allow a good fit to the GeV–TeV observations without a substantial contribution from IC photons. Nevertheless, the shape is uncertain since escaping CRs will be highly anisotropic and the plasma physics of these distributions is an active area of research.⁶

The broadband fit we obtain for our IC core-collapse model is, in fact, significantly better than that obtained in Ellison et al. (2010), particularly for the highest energy High Energy Stereoscopic System (HESS) points. The reason, as we show below, stems from the fact that the pre-SN wind magnetic field, in which the SN explodes, can be considerably lower than 3 μG, allowing electrons to be accelerated to higher energies before radiation losses dominate.

It is important to emphasize that our model and all self-consistent shock acceleration models of SNR J1713 we are aware of accelerate CR ions efficiently. The models we show below place 25%–50% of the FS ram kinetic energy flux into relativistic ions at any instant. Only 0.25% or less of the instantaneous ram kinetic energy flux goes into relativistic electrons. Leptons dominate the observed emission simply because leptons radiate far more efficiently than ions, not because ions are missing. Furthermore, our best-fit parameters for this remnant result in maximum proton energies of ∼10¹⁴ eV. Iron nuclei would be accelerated to ∼26 × 10¹⁴ eV, well into the CR "knee." We make no claim that the spectral shape produced in this particular SNR matches the CR spectrum observed at Earth. As has been known for decades (e.g., Ginzburg & Syrovatskii 1964), CRs observed at Earth come from a mix of many SNRs, with some admixture of other sources (see Meyer et al. 1997), after a long and tortuous propagation in the interstellar medium (ISM). The spectrum produced by any one source will be modified substantially by propagation (see Figure 1 in Ellison & Eichler 1985), and different elements, e.g., hydrogen and helium, may propagate differently (e.g., Blasi & Amato 2011). Despite some recent claims (Adriani et al. 2011), we believe our detailed fit to SNR J1713 is fully consistent with SNRs being the primary source of galactic CR ions at least to energies into the CR knee.

The CR, hydrodynamic, non-equilibrium ionization (CR-hydro-NEI) model of an evolving SNR we use has been described in detail in a number of previous papers. It is essentially the same as that described in Ellison et al. (2010) and references therein except that we have added the escape of high-energy CRs from the shock precursor and the propagation and γ-ray production of these CRs as they diffuse in the CSM. Our method describing the escaping CRs is given in Ellison & Bykov (2011).

The important features of the CR-hydro-NEI model are (1) the hydrodynamics and evolution of the SNR are coupled to the efficient production of CRs via nonlinear (NL) DSA; (2) the NL CR acceleration is calculated using the semi-analytic technique given in Blasi et al. (2005); (3) the CRs that remained trapped within the FS have a normalization and maximum energy cutoff consistent with the normalization and shape of the high-energy CRs that escape upstream; (4) magnetic field amplification (MFA), presumably produced by NL DSA at the FS, is simply parameterized; (5) as described in detail in Ellison et al. (2007), the NL acceleration model determines the heating of the shocked plasma; (6) the electron temperature and NEI state of the shocked plasma are calculated and followed self-consistently with the SNR dynamics in the interaction region between the contact discontinuity and the FS as the SNR evolves; (7) within the SNR, thermal X-ray emission lines and continuum emission from synchrotron, IC, bremsstrahlung, and pion decay from trapped CR electrons and protons are calculated taking into account the SNR evolution and adiabatic and radiation losses;⁷ and (8) the pion-decay emission from escaping CR protons and trapped CRs is absolutely normalized relative to each other and relative to the continuum and line emission in all other bands.

Since, as in our previous applications of the CR-Hydro-NEI model, we assume spherical symmetry, all of our escaping CRs interact with the external shell. In three-dimensional geometry, where the external mass is concentrated in clumps, only some fraction of the escaping CRs would encounter the clumps before diffusing away (see, for example, Lee et al. 2008). Our γ-ray fluxes from escaping CRs interacting with a given external mass are therefore upper limits.

Any realistic description of a SNR requires a large number of parameters. Parameters are required to model the SN explosion and SNR, the DSA mechanism, the CSM, and the diffusion of escaping CRs in the CSM. In Tables 1–3 we list these input parameters, and in Table 4 we list some output values.

Table 2. Input Parameters for DSA and Line Emission

Modela	${\cal E_\mathrm{DSA}}$	Bamp	Kep	αcut	feq	fnorm
	(%)
A	25	8.5	0.01	0.75	1	0.75
B	25	8.5	0.01	0.75	1	0.65
C	50	10	10−3	0.75	0 or 1	0.9
D	25	...	...	...	...	...
E	25	...	...	...	...	...

Note. ^aAll models use f_sk = 0.1, N_cut = 30, and a solar elemental composition.

Download table as:  ASCIITypeset image

Table 4. Output Values at the End of the Simulation

Model	rtot	npa	B0a	B2b	RFS	VFS	Mswept	CR	esc	χinj
		(cm−3)	(μG)	(μG)	(pc)	(km s−1)	(M☉)	(ESN)	(ESN)
A	4.6	0.013	0.22	10	8.9	4200	4.2	0.15	0.013	3.87
B	4.6	0.06	0.23	10	8.9	4200	4.2	0.15	0.013	3.87
C	5.6	2.0	0.7	75	9.8	850	4.3	0.17	0.03	4.0
D	4.6	0.01	...	...	9.7	2500	0.025	4 × 10−3	3 × 10−4	3.9
E	4.6	1.5 × 10−5	...	...	19	9000	0.04	6 × 10−3	5 × 10−4	3.7

Notes. ^aValue just upstream of the FS at the end of the simulation. ^bMagnetic field just downstream from the FS at the end of the simulation.

Download table as:  ASCIITypeset image

Besides the assumed age and distance to SNR J1713, Table 1 gives values for the SN explosion energy E_SN, and ejecta mass M_ej, as well as the power-law index of the initial ejecta mass distribution, n (see Ellison et al. 2004, for a full discussion of n). All of our examples here assume a pre-SN wind with a speed V_wind, mass-loss rate dM/dt, and temperature T_wind, all of which are assumed constant. The parameter σ_wind shown in Table 1 determines the wind magnetic field at a radius R from the explosion according to

$$\begin{equation} B_0(R) = \frac{(\sigma _\mathrm{wind}V_\mathrm{wind}dM/dt)^{1/2}}{R}, \end{equation} \tag{ 1 }$$

as in Chevalier & Luo (1994) and Ellison & Cassam-Chenaï (2005). The constant parameter σ_wind is the ratio of magnetic field energy density to kinetic energy density in the wind and can be related to properties of the star by (e.g., Walder et al. 2011)

$$\begin{equation} \sigma _{\rm wind}\propto \frac{B_{\ast }^2 R_{\ast }^2 }{(dM/dt)v_w}\times \left(\frac{v_r}{v_w}\right)^2. \end{equation} \tag{ 2 }$$

Here, B_* denotes the surface magnetic field of the star, R_* is the stellar radius, v_r is the star rotation velocity, and v_w = V_wind is the terminal speed of the wind. Obtaining σ_wind ∼ 0.03 by fitting the spectrum of SNR J1713 as we do constrains the progenitor star parameters. As noted by Walder et al. (2011), values of σ_wind ≪ 1 indicate that the stellar wind dominates the magnetic field, producing a roughly radial field far from the star.

All of the parameters in Table 2 have been described previously (see Ellison et al. 2010, and references therein). Briefly, ${\cal E_\mathrm{DSA}}$ is the fraction of ram kinetic energy flux that is placed in superthermal particles at any instant,⁸ B_amp is an amplification factor for the shocked magnetic field, K_ep is the ratio of electrons to protons at relativistic energies, f_sk is the fraction of the FS radius used to determine the maximum CR momentum, α_cut is a factor that determines the shape of the particle cutoff at high momentum for trapped CRs, and N_cut (along with α_cut) determines the shape of the escaping CR distribution (see Ellison & Bykov 2011, for a detailed discussion of α_cut and N_cut). The parameter f_eq = 1 indicates that Coulomb equilibration is used for electron heating, while f_eq = 0 indicates that electrons equilibrate instantly with protons downstream from the shock (see Ellison et al. 2007), and f_norm is the overall normalization applied to the photon emission to match the observations. In all cases, a solar elemental composition has been assumed for the CSM.

The values given in Table 3 determine the characteristics of the CSM and the diffusion of the escaping CRs. As described in Ellison & Bykov (2011), the diffusive mean free path for the escaping CRs in the CSM is parameterized by

$$\begin{equation} \lambda _\mathrm{CSM}= \lambda _\mathrm{CSM,0}(r_g/r_{g,0})^{\alpha _\mathrm{rg}}[n_\mathrm{CSM}(R)/n_0]^{-\beta _n}. \end{equation} \tag{ 3 }$$

Here, r_g = pc/(eB) is the gyroradius, n_CSM(R) is the CSM proton number density,⁹ and α_rg and β_n are parameters. For scaling, we use n₀ = 1 cm⁻³, r_{g, 0} = 10 GeV/(eB_{CSM, 0}), and B_{CSM, 0} = 3 μG. We note that instead of the n_CSM term, we could have scaled λ_CSM with $(B/B_\mathrm{CSM,0})^{-\beta _n}$, or even a combination of density and magnetic field terms. These are essentially equivalent parameterizations unless the connection between background field, ambient density, and wave generation by streaming CRs is specified. The normalization of the CSM diffusion coefficient, D_{CSM, 0} = λ_{CSM, 0} c/3, can be estimated from CR propagation studies (see, for example, Ptuskin et al. 2006; Gabici et al. 2009). For example, with D_{CSM, 0} = 10²⁷ cm² s⁻¹, n_CSM = 0.01 cm⁻³, α_rg = 0.5, and β_n = 1, λ_CSM ∼ 1 pc at 1 GeV, consistent with the fits of Ptuskin et al. (2006). In general, the stronger the diffusion (i.e., the smaller λ_CSM) the greater the γ-ray emission will be in the external material.

The values n_uni and n_shell are the proton number densities for the uniform CSM beyond the dense shell and for the dense shell, respectively. The values M_shell and R_shell are the mass of the dense shell and its inner radius, respectively, and B_shell is the magnetic field in the shell. As shown in Figures 1 and 3, we smooth the transition from the pre-SN wind to the dense shell.

Zoom In Zoom Out Reset image size

We note that within the FS we assume Bohm diffusion for the CRs with a mean free path λ ∼ r_g, which is very much smaller than λ_CSM. This is reflected in the sharp drop in λ within the FS as shown in the bottom panels of Figures 1 and 3.

3.1. Pre-SN Wind Interaction

For our core-collapse model A, we take the SN explosion energy to be E_SN = 10⁵¹ erg, the ejecta mass M_ej = 3 M_☉, and assume a slow, dense, pre-SN wind with a mass-loss rate dM/dt = 10⁻⁵ M_☉ yr⁻¹and wind speed V_wind = 20 km s⁻¹. Our model is a simplified description that might resemble what happens after an early-type star with a fast wind creates a large, low-density bubble before evolving into a red supergiant with a much slower wind (see, for example, Chevalier 1999). As we show below, the critical conditions that result in a good leptonic fit are that the density is relatively low and the B-field in the wind is lower than the normal ISM field. Other than this, none of our conclusions depend critically on particular wind parameters.

To determine the unshocked magnetic field as a function of radius, R, in the pre-SN wind, we take σ_wind = 0.03 in Equation (1). At the assumed age of SNR J1713 (i.e., t_SNR ≃ 1630 yr), the FS has not yet reached the dense material of the swept-up wind. The situation is shown in Figure 1 where, in the top two panels, the proton number density, n_p, and the magnetic field, B, are plotted as functions of radius, R, from the center of the SNR. The dashed curves in the top two panels are the density and magnetic field profiles at the start of the simulation. The solid curves are these profiles at t_SNR = 1630 yr. Parameters have been chosen so the SNR radius is ∼9 pc at t_SNR = 1630 yr, consistent with a distance to SNR J1713 of ∼1 kpc and an angular size of ∼60 arcmin. As seen in the second panel, the pre-SN wind magnetic field just upstream of the FS, as determined by σ_wind, is ∼0.2 μG at t_SNR = 1630 yr and this is increased by compression and amplification to ∼10 μG immediately downstream.

At a radius beyond the FS, we have placed a dense shell with a total mass ∼100 M_☉ and the third panel in Figure 1 gives the mass within R. A shell mass of 100 M_☉ is an extreme case and any shell from swept-up, pre-SN wind material would be considerably less. The fourth panel in Figure 1 shows the number density of escaping CRs that have diffused beyond the FS and the bottom panel shows the scattering mean free path beyond the FS, λ_CSM, for 10 GeV (solid black curve) and 1 TeV (dashed red curve) CRs. Note that we have set B = 3 μG in the shell and beyond. Since we do not calculate the synchrotron emission from electrons beyond the FS, the value of B beyond the FS is unimportant.

In Figure 2, we show the radiation produced by the trapped and escaping CRs, along with the observations of SNR J1713 at radio (Acero et al. 2009), X-ray (Suzaku; Tanaka et al. 2008), GeV (Fermi-LAT; Abdo et al. 2011), and TeV energies (HESS; Aharonian et al. 2007). As described in Ellison et al. (2010), this is a "best-fit" result where we have varied parameters to obtain a good match to the broadband data. The essential features of the comparison are (1) IC emission from relativistic electrons interacting with the cosmic microwave background (CMB) strongly dominates the GeV–TeV emission; (2) the quality of this broadband fit with a pre-SN wind, including weak X-ray emission lines consistent with the Suzaku observations, is superior to that obtained with a uniform CSM model by Ellison et al. (2010); (3) the escaping CRs (black dotted curve), even though they diffuse through a region with 100 M_☉ of material, give a small contribution to the TeV pion-decay emission; and (4) even though there are a lot of relatively free parameters for this fit, the critical result that IC dominates over pion decay at GeV–TeV energies is robust.

Zoom In Zoom Out Reset image size

We do not show a fit where parameters have been optimized for pion-decay dominance at GeV–TeV energies because the result is essentially the same as shown in Ellison et al. (2010) for the homogeneous case. As is clear from Figure 2, if K_ep is decreased to allow pion decay to dominate IC, the overall density will have to be increased substantially to produce enough flux at GeV–TeV energies. Any increase in density (here, for the pre-SN wind through varying dM/dt and/or V_wind) will increase the X-ray line emission along with the pion-decay emission, producing a conflict with the Suzaku observations. In fact, as long as the contribution from the escaping CRs remains well below the IC, all fitting parameters, such as the efficiency for DSA, ${\cal E_\mathrm{DSA}}$, are constrained here just as in the fits shown in Ellison et al. (2010).

We note that the parameter f_sk determines the precursor length in our model but that we do not explicitly model the spatial properties of the precursor in our implementation of the Blasi et al. DSA model (see Caprioli et al. 2009, and references therein, for an updated version of the model that does include the precursor explicitly). We do, however, account for all CRs that are accelerated and, except for the escaping CRs, assume they are all trapped behind the FS. We, therefore, overestimate the pion-decay flux since the trapped CRs interact with the dense shocked plasma rather than the thinner precursor material. The only possible case where the precursor CRs might enhance the pion-decay emission is if the precursor is interacting with external material denser than the shocked plasma. We consider this case unlikely because it requires fine tuning. If the external material is fully outside of the precursor, it is contained in the case we discuss in Section 3.2 below. If the FS is also interacting with the external material, it is the case we discuss in Section 3.3 below. Only if the precursor, but not the FS, is impacting dense material is there a possibility that the pion-decay emission will be greater than we estimate. With a precursor length determined by f_sk = 0.1, we consider this to be unlikely. While larger values of f_sk are possible, they imply strong self-generated turbulence far upstream where the density of accelerated CRs has dropped significantly from that at the shock (this is, in effect, CR "dilution" as discussed by Berezhko et al. 1996a, 1996b). It is noteworthy that neglecting the spatial properties of the precursor is less likely to overestimate the electron contribution to the GeV–TeV emission since radiation losses will prevent the highest energy electrons from streaming far upstream.

There are two important improvements over the broadband IC fit given in Figure 4 of Ellison et al. (2010). One is that the newer Fermi-LAT data for SNR J1713 now clearly favor an IC model, whereas the preliminary Fermi-LAT data available for Ellison et al. (2010) were less clear. The second improvement is in our match to the highest energy HESS points. In the constant ISM model used in Ellison et al. (2010), the IC fit fell below the highest energy HESS points. Now, with our core-collapse model, we are able to fit the highest energy points successfully with only the CMB photon field. This improvement comes about because the magnetic field at the FS is lower in the core-collapse case and electrons can obtain a higher energy before synchrotron losses dominate. Magnetic field amplification is still important (we fit the data with B_amp = 8.5 for this case) but starting with a lower ambient field is advantageous.

3.2. External Molecular Cloud Interaction

The contribution from the escaping CRs to the TeV emission will increase as the external target material increases, as would be the case if the escaping CRs from the SNR diffused into a nearby molecular cloud (e.g., Aharonian & Atoyan 1996; Ptuskin & Zirakashvili 2005; Gabici & Aharonian 2007; Gabici et al. 2009; Caprioli et al. 2010; Ohira et al. 2011). In Figures 3 and 4, we show an example (Model B) where the mass external to the FS is 10⁴ M_☉ with a density of ∼100 cm⁻³. All other parameters are the same as for Model A and again we note that our model is spherically symmetric so that all escaping CRs interact with the outer shell. Now, the pion-decay emission from the escaping CRs at TeV energies is well above the pion decay from the trapped CR protons and is comparable at the highest energies to the IC from the trapped CR electrons.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While a larger external mass would clearly result in pion decay dominating the IC, it remains to be seen if a satisfactory fit for SNR J1713 can be found with just pion decay. The first problem concerns the shape of the pion-decay emission from the escaping CR distribution. The distributions shown in Figures 2 and 4 are much too narrow to produce a good fit to both the Fermi-LAT and HESS fluxes and a substantial contribution from IC is required to produce a good fit. We note, for example, that Zirakashvili & Ptuskin (2008) and Caprioli et al. (2010) find similar narrow escaping CR distributions, as does Vladimirov et al. (2006) with Monte Carlo simulations that directly determine the escaping flux, but that Gabici et al. (2009) and Ohira et al. (2011) assume broader distributions when averaged over the age of the remnant.¹⁰ As was discussed in Ellison & Bykov (2011), the shape of the escaping CR particle distribution depends on the details of how the highest energy particles, trapped and escaping, generate turbulence, particularly via long-wavelength effects. The long-wavelength, wave–particle interactions that are important for these high-energy, escaping particles involve highly anisotropic distributions that are only beginning to be modeled self-consistently (e.g., Caprioli et al. 2010; Bykov et al. 2011; Schure & Bell 2011). It is certainly possible that broader distributions of escaping CRs than we show may come from these calculations or others and that a better fit to the GeV–TeV flux without the IC component could be obtained. Nevertheless, we feel it is unlikely that escaping CRs streaming away from a relatively young SNR will have a distribution much broader then we show here.

Regardless of the details of the wave generation, however, the shapes of the trapped CR ions and the escaping ones are related since the turnover in the trapped ion spectrum comes about as CRs escape from the FS. In the case of a high enough B-field, the shape of the electron turnover is determined mainly by radiation losses. With any attempt to broaden the spectral cutoff to allow a match to the GeV–TeV emission, care must be taken so that the X-ray synchrotron is not broadened beyond acceptable limits. This is particularly important for a core-collapse model where the effective magnetic field is less than in a homogeneous model.

Another problem that must be addressed if γ-ray production in external material is important concerns the observed spatial coincidence of γ-ray and X-ray emission. If γ-rays from pion production in external clouds were to dominate the γ-ray emission, one would not expect the observed good agreement between the X-ray and γ-ray morphologies (i.e., Aharonian et al. 2006).

3.3. SNR Impacting External Material

A possibility that has been discussed for SNR J1713 (e.g., Berezhko & Völk 2006, 2010) is that the FS is currently impacting a pre-SN shell or molecular cloud at a radius of ∼10 pc. In this case, timescales considerably less than t_SNR become important as the FS moves into the steep density ramp of the shell or cloud edge. Since there is a time delay between the shock heating of the plasma and the non-equilibrium production of X-ray emission lines, the question naturally arises: is it possible for GeV–TeV emission to be produced before strong X-ray emission lines are generated?

In Figure 5, we show our model C as a time sequence of the SNR density profile and the photon emission as the FS runs into a density gradient. The parameters are similar to Model B except we have moved the inner edge of the dense shell to ∼9.6 pc so that the SNR radius is about 10 pc at 1600 yr to be consistent with observations. We have also lowered K_ep to 10⁻³ and increased the wind speed to 30 km s⁻¹ from 20 km s⁻¹ and made other changes (see Tables 1–3) in the input parameters in the attempt to obtain a good fit to the X-ray and GeV–TeV observations at 1600 yr with pion decay, from a combination of trapped and escaping CR protons, dominating the GeV–TeV emission at t_SNR = 1600 yr. As is clear from Figures 2 and 4, the mass concentrated in the dense shell must be ≫100 M_☉ to make the pion-decay contribution from the escaping CRs comparable to IC and we use M_shell = 10⁴ M_☉ for Figure 5.

Zoom In Zoom Out Reset image size

In Figure 5, the black dashed curves in the left-hand plots are the initial density profiles and the red solid curves are the density profiles at the ages shown. The right-hand plots show the broadband emission at the indicated ages. The important properties of Figure 5 are the following.

• 1.  
  In order to not overproduce the thermal X-rays, the density upstream of the FS at ∼1600 yr must be ≲ 1 cm⁻³, forcing the FS to be ascending the density ramp and not in the high-density plateau.
• 2.  
  The pion-decay emission from the escaping CRs varies weakly with t_SNR since it depends largely on the constant amount of external mass in the dense shell.
• 3.  
  The pion-decay emission from the trapped CRs increases rapidly as the FS ascends the dense shell because it scales ∼n²_p.
• 4.  
  As the FS moves into the "molecular cloud," the X-ray synchrotron and TeV IC increase less rapidly than the pion decay from the trapped CRs because radiation losses increase due to the increasing magnetic field in the cloud. This is particularly noticeable in the synchrotron between 1400 and 1600 yr.
• 5.  
  Between 1600 and 1800 yr, the relative intensity of thermal X-rays increases by more than a factor of 10 relative to the non-thermal continuum emission.
• 6.  
  The distance between the forward and reverse shocks (RSs) drops substantially as the FS runs up the density gradient.

In Figure 6, we show the thermal X-rays along with the synchrotron emission (solid blue curves) in the X-ray band for Model C. The solid black histograms are calculated assuming that electrons immediately behind the shock are cold and then heated via Coulomb collisions with the shock-heated protons. The dashed green histograms are calculated assuming the electrons are heated instantly with the protons. While the method of electron heating produces important differences in the details of the emission lines, all major conclusions for broadband models are independent of this heating (see Ellison et al. 2007, 2010, for a full discussion of our model for electron heating). The calculation of the emission lines includes the non-equilibrium effects as the plasma is shock heated and compressed and then expands and cools behind the FS. For this set of parameters, until ∼1600 yr, the emission lines lie mainly at or below the synchrotron continuum while at 1800 yr, the line flux is more than 10 times the synchrotron flux at ∼0.1 keV. The variation in line intensity between 0.2 and 1 keV over the 400 yr period from 1400 to 1800 yr is greater than an order of magnitude, while the variation in synchrotron intensity in this energy range is less than a factor of three.

Zoom In Zoom Out Reset image size

Also evident in Figure 6 is that the high-energy synchrotron emission decreases with increasing age. We note that Patnaude et al. (2011) have observed a decline in the non-thermal emission from Cassiopeia A on timescales ≲ 10 yr with a 1.5%–2% yr⁻¹ decline in the 4.2–6.0 keV range. Such rapid changes in synchrotron emission could result from a number of causes as the FS enters a steep density gradient. These include the slowing of the FS, an increase in the ambient magnetic field strength, as is the case in Figure 5, and/or the damping of magnetic turbulence in the weakly ionized dense material.

3.4. High-velocity Pre-SN Wind

Another possibility is that SNR J1713 exploded in a fast wind with V_wind ∼ 1000 km s⁻¹, and dM/dt = 10⁻⁷ M_☉ yr⁻¹, typical of an early B-type or late O-type star (e.g., de Jager et al. 1988). If we accept the constraints that the remnant currently has an age of ∼1600 yr, a radius of ≲ 10 pc, and an explosion energy of E_SN ∼ 10⁵¹ erg, it is clear that the only possibility is that the FS is now impacting dense material that stands ∼10 pc from the explosion center.

To illustrate this, we show in Figure 7, examples (Models D and E), where V_wind = 500 km s⁻¹, dM/dt = 10⁻⁶ M_☉ yr⁻¹, and M_ej ∼ 10 M_☉. The position of the FS at 1600 yr will be well beyond 10 pc unless the FS is slowed by dense material. We note that we have chosen extreme values for Figure 7; a faster V_wind, a lower dM/dt, and/or a smaller M_ej will result in an FS beyond 20 pc at 1600 yr. The slower speed we use also encompasses the case where the CSM contains regions of both high- and low-speed pre-SN winds.

Zoom In Zoom Out Reset image size

The black curves in Figure 7 show the case where the shell edge is at R_shell ∼ 9.5 pc, while the red curves show the results for the same input parameters except that the dense shell is beyond the radius the FS obtains at t_SNR = 1630 yr. The same arguments used to exclude the FS penetrating the dense shell in Figure 5 apply for a fast wind. The constraints of age, distance, and radius are not consistent with the lack of X-ray line emission if the FS is interacting with a dense shell.

We have generalized our CR-hydro-NEI model of an evolving SNR to include the production of radiation in a non-homogeneous CSM such as would exist with a pre-SN wind or for a SN exploding near a molecular cloud. The model includes the production and escape of high-energy CR protons, as well as CR electrons and protons that remain trapped within the FS. As the escaping CR protons diffuse in the CSM, they produce γ-rays by pion decay and this production is added to the synchrotron, IC, and pion-decay emission from the trapped CRs. As in our previous work (e.g., Ellison et al. 2007; Patnaude et al. 2009; Ellison et al. 2010), we simultaneously calculate the thermal X-ray line emission with the broadband continuum allowing strong constraints to be placed on some of the many parameters required to model the SNR J1713.

CR protons that escape the SNR and interact with external material produce γ-rays without producing thermal X-rays. If escaping CRs are not considered, homogeneous models clearly favor IC from trapped CR electrons as the mechanism responsible for the GeV–TeV emission in SNR J1713. Including escaping CRs opens the possibility that the constraint on ambient density imposed by the lack of X-ray line emission in SNR J1713 can be satisfied and still have the GeV–TeV emission produced predominately by CR protons in the external material.

As in our previous work, we do not calculate the thermal X-ray emission at the remnant RS, nor do we explicitly include the acceleration of elements heavier than protons at the FS. We do include a 10% by number contribution of helium in the target nuclei. The calculations of Caprioli et al. (2011) include heavy nuclei in the shock modification and show an increase in the pion-decay emission by as much as a factor of 2.5 along with changes in the shape of the emitted spectrum. However, the enhancement of pion-decay emission from heavy elements is not enough to modify our main conclusion that leptons dominate the GeV–TeV emission, even without considering the thermal emission from the RS. Including the RS thermal emission would produce an even greater enhancement of IC over pion decay at GeV–TeV energies.

We have considered four cases. The first is when the external material is a shell from a slow pre-SN wind, in which case the external mass is <100 M_☉. The second is when a larger external shell of M_shell = 10⁴ M_☉ lies near but outside of the FS at t_SNR = 1630 yr. The third case examines what happens when the FS is interacting with the dense 10⁴ M_☉ shell. The fourth case considers a SNR in a fast pre-SN wind interacting with a 100 M_☉ shell.

Our major conclusion is that it is not possible to obtain a satisfactory broadband fit to SNR J1713 dominated by pion decay if the remnant is isolated. An excellent fit can be obtained (Figure 2) with IC dominating the GeV–TeV emission, confirming that an isolated SNR J1713, whether from a core-collapse or thermonuclear SN, produces most of its GeV–TeV emission by IC radiation. It is possible, but unlikely, for pion decay from escaping CRs to dominate the GeV–TeV emission if this remnant is near or interacting with ≳ 10⁴ M_☉ of external material.

If the SNR is near a massive region ∼10⁴ M_☉, then pion decay from escaping CRs can, in principle, outshine the IC from trapped CR electrons at TeV energies. However, there are fundamental uncertainties in the properties of escaping CRs having to do with the shape of the escaping CR distribution and their diffusion in the partially ionized CSM. The model we use produces a distribution that is too narrow to match the combined Fermi-LAT and HESS observations with pion decay alone and a substantial contribution from IC is required for a satisfactory fit, see Figure 4. The highest energy CRs, whether trapped or escaping, will have highly anisotropic distributions and substantial work is needed to better understand wave generation in such distributions before a firm prediction for the shape can be made.

Our conclusion that the GeV–TeV emission from an isolated SNR J1713 must be dominated by IC is essentially independent of the CR diffusion parameters in the CSM since escaping CRs contribute an insignificant fraction of the emission for M_shell ⩽ 100 M_☉. However, if the escaping CRs are interacting with enough mass (≳ 10⁴ M_☉), the diffusion properties will be important. In that case, the energy dependence of D_CSM, as well as the normalization, will influence the shape and intensity of the pion-decay emission. As the CRs diffuse through the molecular cloud, D_CSM will be determined by wave–particle interactions in a dense, partially ionized medium (e.g., Reville et al. 2007, and references therein). A detailed analysis of this propagation is beyond the scope of this paper.

A new aspect of our work is the modeling of the remnant as the FS impacts a density gradient as shown in Figures 5 and 6. These figures show that the emission, and particularly the relative thermal and non-thermal fluxes, will vary substantially on short timescales ∼100 yr depending on how long the shock has been interacting with the cloud. We note that an important approximation of our CR-hydro-NEI model will modify this conclusion somewhat. As the FS runs into a steep density gradient, the mass flux crossing the shock will increase and, depending on the change in shock parameters, more mass will be injected into the DSA mechanism. In reality, it will take some time for this additional mass to be accelerated to TeV energies but we assume the acceleration occurs essentially instantaneously. The X-ray line emission is calculated consistently with time so the relative increase in thermal over TeV emission shown in Figure 5 would actually be greater than what is shown.

Significantly, we have found that the core-collapse scenario presented here affords a better fit to the broadband SNR J1713 observations than the homogeneous Type Ia SN model presented in Ellison et al. (2010). The lower pre-SN wind magnetic field present in the core-collapse model yields a lower post-shock field, even with strong MFA, allowing electrons to be accelerated to higher energies than in the homogeneous ISM. These higher energy electrons produce IC emission against the CMB consistent with the highest energy HESS points.

The low magnetic field value we find is not necessarily inconsistent with the rapid time variations observed for SNR J1713. Uchiyama et al. (2007) observed ∼1 yr variations in X-ray synchrotron emission and interpreted this as the radiation loss timescale for TeV electrons in ∼1 mG fields. However, other estimates yield lower values (see references in Ellison & Vladimirov 2008), and there is an alternative explanation for rapid time variations that does not require such large fields. Bykov et al. (2008) showed that a steeply falling electron distribution in a turbulent magnetic field can produce intermittent synchrotron emission consistent with the Uchiyama et al. (2007) observations. In the Bykov et al. (2008) model, root-mean-square fields of 10's of μG are adequate to explain this time variation and are consistent with our findings here.

It is important to note that while we feel the case is all but closed for SNR J1713, particularly considering the recent Fermi-LAT observations, we make no claim that all SNRs will show IC dominance at GeV–TeV energies. The relative brightness of pion decay versus IC depends largely on environmental parameters and pion decay may well dominate in some SNRs, particularly if they show strong X-ray line emission. It is also possible that, in a particular remnant, high-energy emission in some areas will be dominated by IC and in others by pion decay.

It also must be noted that in this paper we only consider a smooth, spherically symmetric ambient density distribution. Clumpiness of the matter, which may originate from the shock interaction with either a clumpy stellar wind or molecular cloud clumps, may affect the ratio of the X-ray thermal line luminosity to the GeV–TeV γ-ray luminosity. Indeed, the core of a shocked, dense clump may have a temperature below the X-ray emitting regime if the density contrast between the clump and the interclump matter is large enough. Hadrons accelerated at the FS could then produce γ-rays in these clumps without the accompanying X-ray line emission. To calculate the ratio of the γ-ray emission due to pion decay to the thermal X-ray emission one must model the clump evaporation that results in a hot X-ray emitting corona surrounding the shocked clump. This corona is sensitive to the uncertain magnetic field structure of the clump and of the shocked ambient matter. In the case of older, evolved SNRs in molecular clouds like IC 443, it has been shown that clumps can produce significant γ-ray emission from both trapped electrons (e.g., Bykov et al. 2000) and protons (e.g., Castro & Slane 2010). The constraints that come into play when clumpiness is considered for young SNRs require a thorough investigation of the clump evaporation effect and is beyond the scope of this paper.

A final important point is that even though IC from relativistic CR electrons dominates pion decay from relativistic ions, all of our models, and all nonlinear models we are aware of, show the efficient acceleration of protons. Far more energy is put into relativistic protons than into electrons by our DSA models. Our work does not in any way call into question the widely held belief that SNRs are the primary sources of CRs, at least up to the so-called knee at 10^{15 − 16} eV. The evidence that some individual SNRs produce CR ions with high efficiency, while indirect, is compelling (see, for example, Reynolds & Ellison 1992; Hughes et al. 2000; Warren et al. 2005; Helder et al. 2009; Morlino & Caprioli 2011).

D.C.E. acknowledges support from NASA grants ATP02-0042-0006, NNH04Zss001N-LTSA, 06-ATP06-21, and NNX11AE03G. D.P. and P.S. acknowledge partial support from NASA Contract NAS8-03060 and Grant NNX09AT68G. A.M.B. was supported in part by the Laboratory of Astrophysics with Extreme Energy Release at the Politechnical University, St. Petersburg, RBRF grant OFIm 11-02-12082, and by the RAS Presidium Program. D.C.E. and P.S. are grateful to Aspen Center for Physics where part of this work was done when the authors were participating in a summer program.